分形几何第四章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