分形几何第四章new文档格式.docx
《分形几何第四章new文档格式.docx》由会员分享，可在线阅读，更多相关《分形几何第四章new文档格式.docx（26页珍藏版）》请在冰豆网上搜索。
双曲余弦：
coshx=[e^x+e^（-x）]/2,
双曲正弦：
sinhx=[e^x-e^（-x）]/2.
在复迭代程序中，这些基本运算经常用到，最好将它们都编制成函数（过程），存到一个单元中去，使用时可以把它们当成“标准函数（过程）”调用。
设单元为COMPLEX.TPU，COMPLEX.PAS的内容如下：
unitcomplex;
｛PASCAL单元：
复数四则运算｝
interface
procedureadd（x1,y1,x2,y2:
real;
varx3,y3:
real）;
{复数加法}
proceduresub（x1,y1,x2,y2:
{复数减法}
proceduremult（x1,y1,x2,y2:
{复数乘法}
procedurecdiv（x1,y1,x2,y2:
{复数除法}
functioncosh（x:
real）:
{双曲余弦函数}
functionsinh（x:
{双曲正弦函数}
procedurecsin（x,y:
varx1,y1:
real）;
{复数正弦}
procedureccos（x,y:
{复数余弦}
procedurecexp（x,y:
{复数指数}
implementation
｛复数加法｝
{计算z_3=z_1+z_2，
其中:
z_1=x_1+iy_1,z_2=x_2+iy_2,z_3=x_3+iy_3}
begin
x3:
=x1+x2;
y3:
=y1+y2
end;
｛复数减法｝
{计算z_3=z_1-z_2}
=x1-x2;
=y1-y2
｛复数乘法｝
{计算z_3=z_1*z_2}
=x1*x2-y1*y2;
=y1*x2+x1*y2
｛复数除法｝
{计算z_3=z_1/z_2}
var
denom:
{分母}
denom:
=x2*x2+y2*y2;
=（x1*x2+y1*y2）/denom;
{实部}
y3:
=（x2*y1-x1*y2）/denom{虚部}
functioncosh（x:
{计算双曲余弦cosh（x）}
cosh:
=（exp（x）+exp（-x））/2.0
functionsinh（x:
{计算双曲正弦sinh（x）}
sinh:
=（exp（x）-exp（-x））/2.0
procedureccos（x,y:
{计算复数余弦z=cos（x+iy），其中:
x1是z的实部，y1是z的虚部}
x1:
=cos（x）*cosh（y）;
y1:
=-sin（x）*sinh（y）
procedurecsin（x,y:
{计算复数正弦z=sin（x+iy），其中:
x_1是z的实部，y_1是z的虚部}
=sin（x）*cosh（y）;
=cos（x）*sinh（y）
procedurecexp（x,y:
{计算复数指数z=e^（x+iy），其中：
=exp（x）*cos（y）;
=exp（x）*sin（y）
end.{复数四则运算单元结束}
4.2曼德布罗特集
复解析映射总是将复平面分解为两个互不相交的子集合。
一个是稳定的集合，其上的动力学行为是平庸的；
另一个是朱丽亚集J，其上的映射是混沌的。
数学家已经证明朱丽亚集J是完备集（即J与它的导集J′是一回事），并且是不变的。
在复迭代中影响最大的当属迭代z→z^2+c，实际上它只是形式更一般的复解析迭代z_（n+1）=F（z_n）+c的一种，F是一个非线性函数。
显然z→z^2+c也是最简单的一种，它在复迭代中的地位相当于逻辑斯蒂映射x_（n+1）=ax_n（1-x_n）在实迭代中的地位（见第八章）。
考虑一般形式的F，令z=x+iy,c=c_（X）+ic_（Y），其中i表示虚数，i=SQRT（-1）。
分离实部与虚部，具体化迭代关系便有：
x→f（x,y）+c_（X），
y→g（x,y）+c_（Y）.
80年代初曼德布罗特在迭代z→z^2+c时，发现了著名的曼德布罗特集，简称M集。
当时迭代精度和色彩调配均不理想，显现的M集也不好看。
但是过了不久，许多高质量的M集图片纷至沓来，尤以德国布来梅大学动力系统图形室所作的图片最为精美，受到举世赞誉。
随之而来的是，各大学的教师、研究生以及本科生纷纷利用自己的计算机试算复迭代，M集一时泛滥高等院校。
据说有一段时间校方被迫作出规定，不允许利用实验室的公用计算机玩曼德布罗特集。
在当今时代，您自己购买了一台微机，如果不在上面玩一玩M集和J集，实在太遗憾了。
看了本书后，读者一定要亲自试一试。
编写计算M集和J集的程序并不复杂，您可以参照本书给出的PASCAL程序，编制适合您自己使用的更好的程序，您可以用PASCAL，C，C++，Java,甚至BASIC语言。
通常所说的M集是迭代二次函数z→z^2+c产生的，此函数具体化就是
x→x^2-y^2+c_（X），
y→2xy+c_（Y）.
其中z=x+iy,c=c_（X）+ic_（Y）,以横轴x记录实数的实部，以纵轴y记录实数的虚部。
M集合实际上是常数c=（c_（X）,c_（Y））构成的图象。
让c从屏幕左上角开始变化，逐行增加，一直变到屏幕右下角。
如果取的区域是200×
200，则一共要计算40,000个点，把计算的结果用不同的颜色标记下来，就得到一幅图象，这就是M集。
对于不同的c值，如何得到表征迭代性质的不同的结果呢?
容易知道，无穷远处肯定是迭代的一个吸引子，即对于复平面上相当多的初始条件，迭代最终都跑到无穷远处。
但研究发现，在原点附近还存在一个奇特的区域，在迭代过程中此区域永远不会跑掉。
在非严格的意义上，这个不变的集合就是M集，我们的主要任务就是画出这个集合的边界——实际上边界是分形曲线，极其复杂，M集图象的全部魅力就在这里。
在复平面上，我们以原点（0,0）作为参考点，观察迭代过程是否远离原点，以及逃离原点的速度如何。
为此规定一个距离函数
D=x^2+y^2，
其实D有许多不同的取法，以上取法是最普通的。
可以看出，如果D较大，表明迭代点离原点较远，如果D较小，表明迭代点离原点较近。
假设对于任何一个c，迭代都从（x_0，y_0）=（0,0）开始，我们观察迭代点列
（x_1，y_1），（x_2，y_2），（x_3，y_3），…，（x_（100），y_（100）），…
的变化状况。
每一次都计算一下D值，即该点与原点（0，0）的距离（为方便计，这里实际上计算的是距离的平方）。
取一个参考距离R，不妨取得大一些，比如R=40（实际上取20就足够了）。
现在想知道迭代多少次后实际的距离D大于R。
在迭代过程中如果D小于R，则继续让计算机迭代，要规定一个上限，比如300次。
如果迭代了300次后结果仍然不跑掉（即D仍然小于R），则可以近似认为此点属于M集合。
对于迭代次数小于300次的情况，如果迭代10次D就大于R，则标记c点为白色；
如果迭代35次D开始大于R，则标记c点为红色；
如果迭代110次D开始大于R，则标记c点为蓝色，等等。
标色有很多技巧，表面看来好像属于计算机技术，但实际上这属于传统的美术。
懂得传统美术色彩理论的人，在此大有用武之地。
一幅分形图形标上不同的颜色，就有完全不同的视觉效果，虽然本质上具有同样的数学结构。
因此，从这种意义上说，分形图形艺术是传统美术与计算机的结合。
如果掌握了计算机这个工具，甚至完全可以说，这与传统的美术没有本质区别。
4.3朱丽亚集
曼德布罗特集M记录的是整个区域上的c值情况，而朱丽亚集是取一固定的c值后，观察复平面上每一点（x，y）在迭代中的表现，并把结果记录下来。
朱丽亚集简称J集。
因此对于同一个迭代函数，既有M集又有J集。
而且这两种集合的确有密切联系。
一个M集可以对应无数种J集，实际上M集就是所有J集的浓缩。
M集不同部位的形状，反映了对应于该处的J集的形状，于是用M集可以对J集进行分类，至少在计算机图形学的层次上可以这样说。
计算J集时，初始迭代点就不能总取（0，0）了，而是要根据实际位置取实际的（x，y）坐标值。
仍以迭代z→z^2+c为例说明。
先取定一个c值，比如c=（1.0221,0.2433），迭代关系化为
x→x^2-y^2+1.0221，
y→2xy+0.2433.
从屏幕左上角算起，逐行计算，一直算到屏幕右下角。
当然，也可以不取整个屏幕那么大，只选一个200×
200的小区域做。
标色的原理与上面讲M集合时完全类似，此略。
改变常数c的取值，可以得到各式各样的J集。
在源程序中，M集与J集的计算方法十分相似，只需改变两处语句就可以互换为对方，具体说明见表6.1。
表6.1从计算机程序上看M集与J集的区别与联系
初始迭代点迭代关系说明
M集x0:
=0;
y0:
=0；
x1:
=f（x,y）+p0；
y1:
=g（x,y）+q0；
p_0,q_0不断变化
J集x0:
=p0;
=q0；
=f（x,y）+c_X；
=g（x,y）+c_Y；
<
/TD>
c_（X），c_（Y）固定不变
下面给出通用M类集合与J类集合逐级放大程序，此程序没有考虑优化问题，读者可仿照它们编写出更好的程序。
{MandelbrotSetandJuliaSet,Huajie,Fart1995.PAS}
usesgraph,crt,vgafont,shelib;
label10;
label20;
const
MaxIteNum=300;
{实际上可以取得小些，如120}
Distance=40;
scalewidthheight=1；
bps,i,PixelNumber,gd,gm:
integer;
pal:
array［0..256*3-1］ofchar;
f:
file;
p:
pointer;
b:
char;
IteNumber,np,nq,kcolor,scale:
integer;
r,p0,q0,x0,y0,x1,y1,temp,eps,lins1,lins2:
real;
pmin,qmin,pmax,qmax,deltap,deltaq:
ppx,qqy,s:
string;
ProcedurePrepareFile（varFilename:
string）;
BEGIN
writeln（'
Pleaseinputimagefilename'
）;
write（'
Filename='
readln（s）;
IFs='
'
THENREPEAT
YoumustinputatrueDOSfilename!
readln（s）UNTILs<
>
;
Writeln（'
Pressanykeytwotimes!
（IFyouwantTOabortpresskey:
ESC）'
assign（f,s）;
rewrite（f,1）;
blockwrite（f,win.width,sizeof（win.width））;
blockwrite（f,win.height,sizeof（win.height））;
blockwrite（f,bps,sizeof（bps））;
blockwrite（f,pal,sizeof（pal））;
END;
procedureAbort（Msg:
string）;
Writeln（Msg,'
:
'
GraphErrorMsg（GraphResult））;
Halt
（1）;
Clrscr;
IfRegisterBGIDriver（@VGADriver）<
0thenAbort（'
EGA/VGA'
IfRegisterBGIfont（@SansFont）<
0thenAbort（'
SansSerif'
{3}
IfRegisterBGIfont（@TripFont）<
Triplex'
{1}
Dr.LiuHuajie'
ResearchCenterforScience&
Society'
PekingUniversity,100871'
Beijing,CHINA'
Email:
liuhj@'
writeln;
Nowwestart...'
Pleasepressanykeytocontinue'
readln;
TextColor（RED）;
Hint:
-3.0<
pmin<
+2.0'
Pleaseinputpmin='
readln（pmin）;
qmin<
Write（'
Pleaseinputqmin='
readln（qmin）;
qmin:
=-qmin;
{将屏幕坐标转化成正常坐标}
50<
number<
5000'
Pleaseinputthenumberofpixelsononedirection='
Readln（PixelNumber）;
win.width:
=PixelNumber;
win.height:
=round（PixelNumber/scalewidthheight）;
bps:
=8;
scale:
=PixelNumberdiv400+1;
win.xstart:
win.ystart:
pmax:
=pmin+5.8;
qmax:
=qmin+5.8*（win.height-1）/（win.width-1）;
{SetPalette}
FORi:
=0TO20doBEGIN
pal［3*i］:
=chr（255-10*i）;
pal［3*i+1］:
pal［3*i+2］:
=chr（55+10*i）END;
=21TO30doBEGIN
=chr（11）;
=chr（255-21*（i-21））;
=chr（11）END;
=31TO43doBEGIN
=chr（255-11*（i-31））;
=chr（20）;
=chr（20）END;
=44TO53doBEGIN
=chr（60）;
=chr（255-21*（i-44））;
=chr（60）END;
…
=201TO255doBEGIN
pal［i*3］:
=chr（255-20*（201-i））;
pal［i*3+1］:
=chr（0）;
pal［i+3+2］:
=chr（i）END;
{callproceduretowritetheinitialpartofimagefile}
PrepareFile（s）;
gd:
=VGA;
gm:
=VGAHI;
initgraph（gd,gm,'
SetColor（6）;
SetTextStyle（1,0,2）;
OutTextXY（400,5,'
FractalImage!
SetColor（4）;
OutTextXY（400,50,'
Pleasewait...'
REPEAT
deltap:
=（pmax-pmin）/（win.width-1）;
deltaq:
=（qmax-qmin）/（win.height-1）;
FORnq:
=0TOwin.height-1do
FORnp:
=0TOwin.width-1doBEGIN
p0:
=pmin+np*deltap;
q0:
=qmin+nq*deltaq;
IteNumber:
x0:
{forJuliaset，x0:
=q0;
}
10:
x1:
=x0*x0-y0*y0+p0;
=2*x0*y0+q0;
{forJuliaset，useothervaluesintheplaceofp0,q0}
{x1:
=x0*x0*x0-3*x0*y0*y0-0.431183125;
=3*x0*x0*y0-y0*y0*y0+0.01712875;
=sin（x0*x0）-y0*y0/3+0.0136;
=3*x0*y0+0.21;
}
{lins1:
=（1+eps）*（2*x0*cos（pi/10）-y0*sin（pi/10））;
lins2:
=（1+eps）*（y0*cos（pi/10）-x0*sin（pi/10））;
=x0*x0-y0*y0+lins1;
=2*x0*y0+lins2+0.7;
inc（IteNumber）;
r:
=x1*x1+y1*y1;
IFr>
DistanceTHENBEGIN
=1TOIteNumberdoBEGIN
IF（IteNumber>
i）and（IteNumber<
=i+1）THEN
kcolor:
=iEND;
putpixel（npdivscale,nqdivscale,kcolor）;
=chr（kcolorand255）;
blockwrite（f,b,sizeof（b））END
ELSEIF（r<
=Distance）and（IteNumber
grOKTHENHalt
（1）;
SetTextStyle（1,0,4）;
str（pmin,ppx）;
str（-qmin,qqy）;
OutTextXY（10,380,'
px='
Moveto（80,380）;
OutText（ppx）;
OutTextXY（10,420,'
qy='
Moveto（80,420）;
OutText（qqy）;
settextstyle（1,0,2）;
OutTextXY（10,200,'
PressESCtoEXIT!
!
SetViewPort（win.xstart,win.ystart,
win.width,win.height,ClipOn）;
ch:
=ReadKey;
ifch=escthenexit;
UNTIL（ch=esc）;
Sound（500）;
Delay（200）;
nosound;
ClearViewPort;
CloseGraph；
END
升级会员 