主题：【原创】千奇百怪话分形 -- 安德的游戏

分形，混沌与非线性动力学是经常混在一起说的概念。其实从严格意义上来说，非线性动力学描述了一个耗散系统，比如说天气的变化，而一个耗散系统有可能收敛，有可能发散，也有可能既不收敛也不发散，而是在一定的范围内呈现一种不规则的震荡。这时候就出现了混沌。而当把出现混沌状态的系统图形化了以后，往往会出现分形图形的特征。所以，这三者既有联系，又有区别。不过既然讲到了分形，那么就不得不涉及到一些关于混沌的东西。

说到混沌，就一定要提到洛仑兹吸引子。洛仑兹生于1917年（嗯，老先生已经90高龄了），是美国的数学家和气象学家，目前在MIT。他在上个世纪五十年代末六十年代初研究长期气象预报的时候，用计算机进行模拟，结果发现，初始状态的微小差异会在一定时间以后对状态产生非常大的影响。这个发现后来被表达成著名的“蝴蝶效应”：南美的一只蝴蝶扇一下翅膀会在美国引起一场飓风。其实关于蝴蝶效应的表达方式有很多种，大部分都是差不多的意思，我这里只是找了一种比较常见的说法。2004年还拍过一个以此命名的影片。

基本上，这个效应揭示了我们的世界是呈发散性的。以天气预报为例，对于初始状态，如温度，湿度，风速等测量的微小差异，随着系统的演化，预测值和实际值的差异会逐渐变大。即使测量得足够精确，测量的点也不可能做到无限密集。这样局部的不确定性很快就会影响到全局。所以，从理论上来说，长期的准确天气预报是不可能做到的。

一般来说，天气系统的模型是包含几十到上百个参数的偏微分方程，要用大型计算机经过大量的计算才能推导出系统的演化来。而洛仑兹把这个系统简化成只包含三个参数，依然保持混沌的特性。这三个方程是这样的：

dx/dt=σ(y-x)

dy/dt=x(ρ-z)-y

dz/dt=xy-βz

这里面σ，ρ和β可以取任意大于0的数，是决定系统特性的常数。常用的组合是σ=10，β=8/3，而ρ的值就决定了系统是稳定收敛的还是混沌的。当ρ取小于28的数值的时候，系统会收敛，下面的图是ρ取14的时候系统随时间演变的曲线。因为状态量x，y，z有三个，所以状态空间是三维的。这里画的是在x-z平面上的投影。

很明显系统会逐渐收敛到一个稳定的点上。因为曲线的样子很像是被一个点所逐渐吸引过去的，所以叫做吸引子。这张图反映的是一个一般的吸引子。

当ρ的值取大于28的时候，系统会呈现混沌的状态，这时候会出现所谓的奇异吸引子。系统的演变不再是趋于稳定，而是围绕着两个中心点不停地变化。下面这张图显示了当ρ取32的时候系统的变化。仿照三视图的做法，这里给出了在x-z，x-y和y-z平面上的投影。其中左下图是在x-y平面上的投影，一般画洛仑兹吸引子的时候，大多给的是这张图。而左上方的是在x-z平面上的投影，形状看起来很象是一只展翅的蝴蝶。也有说法说“蝴蝶效应”的命名是起源于这张图。

下面这张图，显示了初值的选择对于系统变化的影响。蓝色和红色的曲线唯一的差异是x的初值一个取10，一个取10.1。红色箭头指出的位置是初始点。经过一段时间的演化以后，绿色箭头表示出系统的状态。可以看到状态已经差很多了。如果系统继续演化下去，二者的差异会变得越来越大，以至于完全无关。

其实这些图画起来很简单，比曼德尔布诺特集画起来要容易得多。我把Matlab的程序给出来，有兴趣的人可以自己改变参数和初始状态试试看。

close all

clear all

x0=10;

y0=15;

z0=20;

N=5000;

step=0.01;

sigma=10;

beta=8/3;

rou=32;

state(1,1)=x0;

state(1,2)=y0;

state(1,3)=z0;

for i=2:N

x=state(i-1,1);

y=state(i-1,2);

z=state(i-1,3);

deltax=sigma*(y-x);

deltay=x*(rou-z)-y;

deltaz=x*y-beta*z;

state(i,1)=x+deltax*step;

state(i,2)=y+deltay*step;

state(i,3)=z+deltaz*step;

end

subplot(2,2,1);

plot(state(:,1),state(:,3));

subplot(2,2,3);

plot(state(:,1),state(:,2));

subplot(2,2,4);

plot(state(:,2),state(:,3));

关于混沌与非线性系统的东西讲下去还有很多。不过我们这里说的是分形，所以就到这里为止吧。

现在，我们要来讲一讲曼德尔布诺特集是怎么生成的了。

首先，我们建立一个复平面。如果不知道什么叫复平面的，那个，就不用继续看说明了，后面还有几张图看看，然后就算看完了。

对于复平面上面的任意一点C=x+y*i，和初值Z(0)=0，用以下的迭代公式进行迭代：Z(i+1)=Z(i)^2+C。所有使得序列Z收敛的点C，就构成了曼德尔布诺特集。

这个生成方法够简单的吧？这实际上就是分形的第五个特征，可以用简单的迭代方法生成。不过在实际操作的时候，还是有一些其他的问题需要考虑进去的。

在画出曼德尔布诺特集的时候，用人工的方法进行计算是不现实的，一般都是用计算机来进行迭代的运算。首先要确定一个范围，超出这个范围的点一定不收敛。我们知道在第一次迭代以后，可以得到Z(1)=C。那么如果对于序列中的任意一个Z(i)，可以有以下不等式：||Z(i+1)||=||Z(i)^2+C||>=||Z(i)^2||-||C||=||Z(i)||^2-||C||。对于序列中的任意一个Z(i)，如果||Z(i)||>=||C||（实际上对于i=1，上式一定成立），那么||Z(i)||^2-||C||>=||Z(i)||^2-||Z(i)||=||Z(i)||*(||Z(i)||-1)>= ||Z(i)||*(||C||-1)。那么对于任意||C||>2，||Z(i+1)||>= ||Z(i)||*(||C||-1)>||Z(i)||*(2-1)=||Z(i)||。所以如果||C||>2，那么序列一定发散。也就是说只要计算所有模小于等于2的点就够了。

接下来，因为复平面上面的点是无限稠密的，所以不可能计算所有的点，而且在一定的尺度上面，只要点足够稠密，对于更小的细节反正也画不出来，所以就不需要计算了。计算机做图都是以像素为单位的，所以如果要画1000x1000像素的图，给定实部和虚部各为正负2的范围的话，那需要计算的点就是C=((x-500)+(y-500)*i)/250，其中x和y分别取1到1000作为图像的横纵坐标。给定不同的偏移量，加在x和y上，就画出不同位置的图像。除以不同的除数，就给出不同的放大率（也就是尺度）。x和y取值的范围，就给出了图像的大小。

因为计算机判断收敛不可能计算无限序列，所以通常给定一个迭代次数，如果迭代相应的次数以后，结果的模还是不大于一个某个域值（比如说2），那么就认为是收敛的，否则就认为是发散的。习惯上曼德尔布诺特集内部的点画成黑色。如果发散的点画成白色，那就有了一个黑白两色的曼德尔布诺特集图像了。需要注意的是，图像放大的倍数越多，也就是越要看更精细的部分，相应的迭代次数就要越大，这样才能画出正确的图像。否则有些区域本来是发散的，也被画成在曼德尔布诺特集内部了。

如果想要一幅彩色的图像，只需要把发散的区域根据迭代到超过域值需要的迭代次数映射成相应的颜色就可以了。颜色映射的算法决定了画出来的图像是否好看，这个就需要自己摸索了。

另外，还有一种曼德尔布诺特集的推广，叫做“多布诺特集”（miltibrot set）。迭代的公式稍有不同，是Z(i+1)=Z(i)^d+C。

如果d=3，给出的就是下面的这个图形。

如果d=4，给出的就是下面的这个图形。

到这里，曼德尔布诺特集的部分就算是基本讲完了。跟曼德尔布诺特集相关联的还有一种朱丽亚集（Julia set），不过那个比较复杂，就不放在这里讲了。