6.9 分形(三)

逃逸时间法

本课有几个作图与“逃逸时间法”有关,我们先介绍一下这种方法。

设 f(z) 是复数集C上的一个变换,定义它的迭代为 

设正实数R为阈半径,正整数N为最大迭代次数,z为集合A(常用矩形区域,作为扫描框)内的任一点,R为半径的圆(也可以是某个区域)包含了A。规定V={z∈C||z|>R},即圆的外部。当时间n(当前迭代次数)<N时,若f"(z)不属于V,则进入下一次迭代,若f"(z)∈V,则终止迭代。当迭代次数n=N时终止迭代。定义终止迭代时的时间n为逃逸时间,记为et,n=et时逃逸临界点f"(z)的模记为em,并对点z赋以关联参数et,em的颜色。这种计算机画图法叫做逃逸时间算法,阈半径R也叫逃逸半径。

图霸中对扫描区域A内的每一个点z用for()函数可以快速求出逃逸时间,计算出颜色值,使用分段着色的曲线或轨迹,经过线扫描成为图象,或用分片着色的曲面直接绘制。颜色函数控制色调,如果要使颜色有红绿蓝等三个参数控制,可以用参数颜色,对点迭代后扫描。
例6.9.1 :用逃逸时间法绘制Sierpinski三角形

如图,将外接圆半径为12的正三角形及其内部作为扫描区域,ΔABC外即为陷阱。以三角形中心为原点建系,把ΔABC向内压缩后得到ΔAFE、ΔFBD、ΔEDC的三个变换构成迭代函数系统,可以绘制IFS谢氏三角(回看分形二的练习2),由于是向内压缩变换,三角形越来越小,且总在大三角形内。设点P为ΔABC内一点,现在反过来求这个点在第几次迭代的小三角形内,即逃逸时间。为此,先求出逆变换:

根据点P在三个小三角形ΔAFE、ΔFBD、ΔEDC的哪一个内部,选择对应的逆变换。若点P在ΔDEF中(落入陷阱),则无论选用哪一个逆变换,下次必然到ΔABC外。上图中n=4时发生逃逸,迭代步骤为:

开始时,P在ΔEDC内,它必是由ΔABC内某点按照变换W3得到的,用逆变换求出是点P1,此点在ΔAFE内,用W1逆变换求出点P2,...,P3不在三个角上的小三角形内,可以认为已经落入陷阱,下次变换后必然逃逸到P4。逆变换是伸放变换,点离原点越来越远,总会发生逃逸。

那么如何判断点是否在区域内呢?图霸中的逻辑运算不难做到。点P(x,y)在ΔABC内可以用不等式组表示为

由于在大三角形内,再去判断在哪个小三角形内就方便了。y>3,就在ΔAFE内,x<0在左边,x>0在右边。现在我们给出作图步骤:

法一:扫描曲线

1. 新建参数n,整数,从2到10,参数y,从-6到12。选取y与n,添加方程曲线,输入参数方程:x=t,y=No1,z=0。No1即选取的参数y,是当前点P的纵坐标。三角形下边宽上边尖,根据y可以求出x范围,输入参数t范围从(No1-12)/sqrt(3)到(12-No1)/sqrt(3)。由于页默认单位长20个象素,所以每个象素0.05个单位,x=t,所以精度设定0.05,可以逐点设定颜色。勾选分段着色,点击编辑着色函数(可以在记事本中输入,复制后再粘贴):

for(set(0,0)+set(1,t)+set(2,No1),get(2)>=-6&&get(2)<=sqrt(3)*get(1)+12&&get(2)<=-sqrt(3)*get(1)+12&&get(0)<No2,(if(get(2)>3,(set(1,2*get(1)),set(2,2*get(2)-12)),get(1)<0,(set(1,2*get(1)+6*sqrt(3)),set(2,2*get(2)+6)),(set(1,2*get(1)-6*sqrt(3)),set(2,2*get(2)+6))),accum(0)),get(0)/(No2+2)+0.56)

这个函数用一个for循环计算逃逸时间,是制作的关键。参考下图去理解。由于for()的使用,大大提高了制作及运算效率。

2. 选取曲线,【显示】中设置线宽为1,【追踪】。选取y,添加【动画】按钮,从-6到12,参数增速0.05,点击执行扫描。拖动n=6,y=-6后重新扫描,如下图。

法二:只用一个曲面

添加参数曲面,编辑方程:x=u,y=v,z=for(set(0,0)+set(1,u)+set(2,v),get(2)>-6&&get(2)<sqrt(3)*get(1)+12&&get(2)<-sqrt(3)*get(1)+12&&get(0)<6,(if(get(2)>3,(set(1,2*get(1)),set(2,2*get(2)-12)),get(1)>0,(set(1,2*get(1)-6*sqrt(3)),set(2,2*get(2)+6)),(set(1,2*get(1)+6*sqrt(3)),set(2,2*get(2)+6))),accum(0)),4-get(0)/2)

u从-6*sqrt(3)到6*sqrt(3),v从-5.95到12,步长均为0.05,u、v可扩展,光照模式。自定义颜色:0.56-z/4,裁剪条件:v>=12+sqrt(3)*u||v>=12-sqrt(3)*u,确定,得到带“厚度”的谢氏三角。

M集和J集

例6.9.2 :曼德勃罗集合(Mandelbrot set)

曼德勃罗集合是由曼德勃罗(分形学之父)在1980年提出的一类分形集,简称M集。有人认为 Mandelbrot 集合是“人类有史以来做出的最奇异、最瑰丽的几何图形”,曾被称为“上帝的指纹”。

考虑复函数变换f(z)=z^2+c,从z=0开始迭代,值构成数列 c,c^2+c,(c^2+c)^2+c,…,它们或者延伸到无穷大,或者只停留在有限半径的圆盘内。 如 c=0 时,各项永远是0,并不发散,因此 c=0属于M集;又如c=3i时,对应的数列为3i,−9+3i,,72−51i,2583-7341i,…,数字越来越大,因此3i就不属于M集。对f(z)施以逃逸时间法,将所有不属于M集的点按照其逃逸时间赋予不同的颜色,就可以得到M集的经典彩色图像。一个点属于M集当且仅当它对应的数列中的任何项的模都不大于2,这里的2就是上面提到的圆陷阱的“逃逸半径”。由于数列的的元素有无穷多个,我们只能取有限的迭代次数来模拟,迭代次数达到N时仍未逃逸的区域就是著名的曼德勃罗湖。我更喜欢把它想象为一个海岛,岛边风光无限,值得去仔细欣赏。它远胜无忌小时候生活过的冰火岛,也许小龙女也想到那儿过一过过儿也没有过过的生活。

扫描法制作时与例6.9.1类似,着色函数改为:for(set(0,0)+set(1,t)+set(2,No1),get(1)^2+get(2)^2<=4&&get(0)<64,(set(3,get(1)^2-get(2)^2+t),set(2,2*get(1)*get(2)+No1),set(1,get(3)),accum(0)),get(0)/68+0.55)

扫描范围:x从-2.05到1.05,y从-1.35到1.35,由于矩形域较小,可以在页属性中把单位长设为100(或调整坐标映射关系),参数t增量改为0.01,效果图如上图。

曲面法绘制时参数方程为:x=u,y=v,z=for(set(0,0)+set(1,u)+set(2,v),get(1)^2+get(2)^2<=4&&get(0)<200,(set(3,get(1)^2-get(2)^2+u),set(2,2*get(1)*get(2)+v),set(1,get(3)),accum(0)),(ln(1+get(0))-lg(0.1+ln(1+sqrt(get(1)^2+get(2)^2)))))*0.1,参数范围同扫描法,着色函数为z+0.49。利用z值不仅简化计算,而且曲面是三维的,着色中关联了et和em,更有立体感。我的着色方法学的也不成熟,诸位可以自由修改,各显神通。

M集的结构非常复杂,极其精细,曼德勃罗湖边风光无限,令人流连忘返。利用曲面绘制,并动态改变参数范围,放大局部图形,可以制成动画,发现层层嵌套的自相似结构,体会数学美与艺术美的高度统一。

下面介绍以(x,y)为局部中心,放大一百万倍的制作方法:

1. 二维视图,页属性中单位长改为50。新建参数x,y,表示中心点坐标,代表复数x+yi,再新建参数t,范围从0到7。添加它们的计算:set(7,x)+set(8,y),set(9,10^(-t))。在内置变量中保存坐标和放大倍数的倒数,可以进行高精度计算。图霸度量时使用单精度浮点数,有效数字6个,内部变量使用双精度,有效数字15个。

2.选取t,添加参数曲面,方程为:x=u,y=v,z=for(set(0,0)+set(1,u*get(9)+get(7))+set(2,v*get(9)+get(8))+set(4,get(1))+set(5,get(2)),get(1)^2+get(2)^2<=4&&get(0)<1024,(set(3,get(1)^2-get(2)^2+get(4)),set(2,2*get(1)*get(2)+get(5)),set(1,get(3)),accum(0)),get(0)/5000)

参数范围均为从-2到2,步长0.02(电脑独显性能好的,可以在页属性中设置单位长为100,这时步长用0.01,图形更大)。颜色自定义为4*z+0.6,确定。

3,添加t的计算10^t,表示当前放大的倍数。拖动t,观看放大过程。切换到三维视图,可见它是一个三维曲面,高度正比于逃逸时间。图形生成后,再进行视旋转与放缩会很快。

如果要方便地选取放大区域,可以仿照下图设置8个曲面作为动画的关键帧,每一个比前一个放大10倍。拖动方框中心,坐标存入内置变量,后一个图形动态改变。最后在下一页根据上页的内置变量值绘制放大图,制作动画。

按照下列设置添加曲面,可以把M集绘制到球面上,其中坐标系使用球坐标系,参数方程中r=10+for(set(0,0)+set(1,u)+set(2,v),get(1)^2+get(2)^2<=4&&get(0)<64,(set(3,get(1)^2-get(2)^2+u),set(2,2*get(1)*get(2)+v),set(1,get(3)),accum(0)),get(0)/80)

这样看起来有立体感,更象哪个星球上的绿洲。

使用条件裁剪,按照下列设置,可以制作镂空效果的M集。裁剪条件为:for(set(0,0)+set(1,u)+set(2,v),get(1)^2+get(2)^2<=4&&get(0)<200,(set(3,get(1)^2-get(2)^2+u),set(2,2*get(1)*get(2)+v),set(1,get(3)),accum(0)),get(0)<15||get(0)==200)

使用隐函数绘制等势线或等势面,可以制作M集的边界线或三维的M集。

【构造】隐函数曲面,柱坐标系下方程为:for(set(0,0)+set(1,-z)+set(2,r),get(1)^2+get(2)^2<=4&&get(0)<100,(set(3,get(1)^2-get(2)^2-z),set(2,2*get(1)*get(2)+r),set(1,get(3)),accum(0)),get(0)-80)=0

这些表达式大同小异,理解一个可以自由变化。

例6.9.3 :朱利亚集合(Julia set)

朱利亚集合是以法国数学家加斯顿·朱利亚(Gaston Julia)的名字命名的,它和曼德勃罗集合差不多,也是研究复函数f(z)=z ^2 +c,但这次我们固定c=a+bi,取某一z值(如z=z0),可以得到序列z1、z2、z3…, 这一序列可能发散于无穷大或始终处于某一范围之内并收敛于某一值。我们将使其不扩散的z0=x+yi值的集合称为朱利亚集合。不同的c值下Julia图案的差别很大。


制作步骤:

新建参数a、b,选取后添加曲面,参数方程为:x=4*u,y=4*v,z=for(set(0,0)+set(1,u)+set(2,v),get(1)^2+get(2)^2<=4&&get(0)<400,(set(3,get(1)^2-get(2)^2+No1),set(2,2*get(1)*get(2)+No2),set(1,get(3)),accum(0)),get(0)/500),

u从-1.6到1.6,v从-1.5到1.5,步长0.01(用0.025更好),自定义着色函数z+0.56.

陷阱技术

在函数迭代过程中,我们设置若干个陷阱,当迭代点落入某个陷阱中,则停止迭代,而给对应的点赋予该井所设置的颜色。陷阱其实是复平面上的一些区域,如圆形、环形、多边形、长带形 ,井的形状将决定生成的分形具有独特形状。在复平面设置互不重叠的陷阱,迭代序列落入不同井中,可以设置不同的颜色,也可以设置相同的颜色。但在同一个井中,离井中心越近,颜色越深,离井中心越远,颜色越浅,这样可以产生立体感。三维软件中用点的z值关联此距离,光在上方照射,产生镜面光。单个井可以根据落入井时的迭代深度确定颜色的色调。

例6.9.4 :陷阱法绘制J集

新建参数d,从0.01到0.03,拖动到0.02,添加曲面x=u,y=v,z=for(set(0,0)+set(1,u)+set(2,v)+set(4,1),get(0)<100&&get(4)>No1&&get(4)<4,(set(4,get(1)^2+get(2)^2),set(3,get(1)^2-get(2)^2+0.15),set(2,2*get(1)*get(2)+0.6),set(1,get(3)),accum(0)),if(get(4)>=No1,0.4,get(0)%7/8+2*get(4)+0.8)),u和v均从-1.2到1.2,步长0.01(当前页单位长用100),自定义着色函数h=z.

个性设置中设置灯光正上方射向XOY面。

 

练习:使用映射 f(z)=z^3+c 制作广义的曼德勃罗集,如果熟悉复数的三角形式,可以制作f(z)=z^m+c对应的图形。

 

 

进入下一章

返回帮助目录