迭代与分形

迭代与分形的概念

分形的概念

分形(fractal)本意是不规则的、破碎的。用来描述整体上看极不规则但整体和局部又存在某种相似性的几何对象。曼德勃罗(Mandelbrot)把部分和整体以某种方式相似的形体称为分形。

迭代的概念

计算机进行数值计算的一个基本特点就是要重复成千上万次的计算，这就是所谓的迭代计算，迭代计算是编程最基本的技巧之一，在数值计算以及分形画图中都有广泛应用。

迭代是重复反馈过程的活动，其目的通常是为了逼近所需目标或结果。每一次对过程的重复称为一次“迭代”，而每一次迭代得到的结果会作为下一次迭代的初始值。

迭代计算的方法

迭代计算除了for循环和while循环之外，还可以通过递归调用，或矢量化编程的方法实现。

分形树

如下图所示，将一条线段三等分，在等分点上各画一条长度为原线段长度三分之一的线段，该线段与原线段成固定的夹角(如$\pi /3$)，就会得到一种像树枝一样的图形。

还是老样子，我们先列出几种绘制分形树的算法，读者可自行调试代码，然后感受算法思路。

树程序1

%% 树程序1
clc;clear all;close all
tic
u=[0,0;0,1];    % 线段oc端点的坐标
subplot(3,3,1),
plot(u(:,1),u(:,2)),    % 绘制线段oc
axis([-0.5,0.5,0,1])
theta=pi/6; 
for n=1:6
    uuu=[];
    for I=0:(length(u)/2-1) % 用上次分形树的线段画下一个分形树
        p1=(u(2*I+1,:)*2+u(2*I+2,:))/3;     % 线段的1/3处
        p2=(u(2*I+1,:)+u(2*I+2,:)*2)/3;     % 线段的2/3处
        lp=[cos(theta),-sin(theta);sin(theta),cos(theta)]*...
            [(u(2*I+2,1)-u(2*I+1,1));u(2*I+2,2)-u(2*I+1,2)]/3;
        lp=p1+lp';			% 向左转动的点
        rp=[cos(theta),sin(theta);-sin(theta),cos(theta)]*...
            [(u(2*I+2,1)-u(2*I+1,1));u(2*I+2,2)-u(2*I+1,2)]/3;
        rp=p2+rp';			% 向右转动的点
        uu=[u(2*I+1,:);p1;p1;lp;p1;p2;p2;rp;p2;u(2*I+2,:)];    % 将点排序
        uuu=[uuu;uu];    	% 将第二层for循环所得的点都储存到uuu中
    end
    u=[uuu];    % 将本轮的所有点的集合作为下一轮的初始点
    subplot(3,3,n+1),plot(u(:,1),u(:,2)),   % 绘制新的分形树
    axis([-0.5,0.5,0,1])
end
toc

代码思路

树程序1是比较容易理解的方法。每次迭代时将所有点储存到矩阵中，使用u(:,1)即矩阵u的第一列储存点的横坐标，使用u(:,2)即矩阵u的第二列储存点的纵坐标，再通过 plot(u(:,1), u(:,2)) 将要绘制的线段连接起来。

对于上一轮循环得到的新分支oc，迭代时会产生新的4个点，所有线段点的排列顺序为 (o, a, a, d, a, b, b, c, b, c)，而plot函数绘制时会将相邻两点都连接一次，因此会重复画线，进而导致代码效率下降。

树程序2

树程序2，使用复数表示点并引入了矢量化编程的方法。

%% 树程序2  
clc;clear all;close all
tic 
new=[0,i]; subplot(3,3,1);
plot(new); axis([-0.5,0.5,0,1]);
for k=1:8
    old=new;
    n=length(old)/2-1;      % 上次生成的线段数目
    diff=(old(2:2:end)-old(1:2:end-1))/3;   % 每条线段1/3的长度
    p1=old(1:2:end-1)+diff;
    p2=p1+diff;
    lp=p1+diff*(sqrt(3)/2+1/2*i);		% 通过乘以模为1的复数实现了旋转功能
    rp=p2+diff*(sqrt(3)/2-1/2*i);
    
    % 以下的计算在第一次循环中计算10个端点
    new(1:10:10*n+1)=old(1:2:end-1);    % 各次循环分别计算第1,11,21,...个点
    new(2:10:10*n+2)=p1; 	% 各次循环分别计算第2,12,22,...个点
    new(3:10:10*n+3)=p1;
    new(4:10:10*n+4)=lp;
    new(5:10:10*n+5)=p1;
    new(6:10:10*n+6)=p2;
    new(7:10:10*n+7)=p2;
    new(8:10:10*n+8)=rp;
    new(9:10:10*n+9)=p2;
    new(10:10:10*(n+1))=old(2:2:end);
    subplot(3,3,k+1),
    plot(new)
    axis([-0.5,0.5,0,1]);   
end
toc

使用复数可以同时表示点的横纵坐标，且对复数构成的列表z_list有 plot(z_list)=plot(real(z_list), imag(z_list))。

通过矩阵直接赋值的方法使树程序2相较树程序1减少了一重循环，代码执行效率更高，这也正是矢量化编程的应用。

树程序3

利用分形结构自相似性

%% 树程序3 
clc;clear all;close all
u=[0,i];
subplot(3,3,1);
plot(u)
axis([-0.5,0.5,0,1])
for k=2:8
    m=u/3;      % 将图形缩小到1/3
    uu=[m,...   % 放置缩小后的图形m
        i/3+m*(sqrt(3)*0.5+0.5i),... % 将m在虚轴上向上移动1/3后向左旋转
        m+i/3,...       % 将m在虚轴上向上移动1/3
        2i/3+m*(sqrt(3)*0.5-0.5i),...   % 将m在虚轴上向上移动2/3后向右旋转
        m+2i/3];    % 将m在虚轴上移动2/3
    subplot(3,3,k);
    plot(uu)
    axis([-0.5,0.5,0,1])
    u=uu;
end

colormap winter

该方法妙在利用分形函数自相似性，在每次循环时只将前一次的结果缩小1/3，然后通过加复数平移再乘复数旋转，将前一次结果移到新的位置，重复操作即可得到新的分形结构。

分形的产生

产生分形结构的几种方法

直接编程法 (例如分形树案例1和2)；

L系统 (例如分形树案例3，以及之后的Sierpinski三角形与科赫雪花案例)；

扩散凝聚(蒙特卡洛法)；

复变函数迭代。

L系统

Sierpinski三角形

连接等边三角形的三条边的中点，将三角形四等分，得到四个较小的三角形，舍去中间的三角形，保留周围的三个，然后按此规律继续操作下去，得到的图形就是Sierpinski三角形。

利用分形的特性绘制Sierpinski三角形

clc;clear all;close all
u=[0,1,0.5+0.5i*sqrt(3),0];
subplot(3,3,1)
fill(real(u),imag(u),'g')
for k=1:8
    m=0.5*u;
    u=[m(1:end-k), m+0.5, m+0.25+0.25*i*sqrt(3),0];
    subplot(3,3,k+1)
    fill(real(u),imag(u),'g')
end

科赫(Koch)雪花

将等边三角形的每条边三等分在去掉中间一段，用更小的等边三角形的两边取代它，依次做下去，得到雪花分形图，称为科赫雪花曲线。

科赫雪花的分形不可直接将整体平移+旋转得到，而是需要先画出一条边，然后旋转得到另外两条边。而第一条边则由上一次迭代得到的边变为原来的1/3，然后再三等分在去掉中间一段，用更小的等边三角形的两边取代中间段所得到。

clc;clear all;close all
a=[0,i];
b=a*(-0.5-0.86603i)+i;
c=(a-i)*(-0.5+0.86603i);
subplot(2,2,1)
plot([a,b,c])
axis equal
for k=1:3
    subplot(2,2,k+1)
    a1=a/3;
    a=[a1,a1*(0.5+0.86603i)+i/3,(a1-i/3)*(0.5-0.86603i)+i*2/3,a1+2*i/3];
    b=a*(-0.5-0.86603i)+i;
    c=(a-i)*(-0.5+0.86603i);
plot([a,b,c])
axis equal
end

总结

如果前面三个案例一样，通过某个基本图形(生成元)反复迭代而生成分形的方法称为L系统。L系统的精髓就是给出一个起点(种子)和替换规则，连续利用该替换规则就能生成新一代系统。

扩散凝聚

DLA (Diffusion-limited Aggregation，扩散限制凝聚) 是由Witten和Sander于1981年共同提出来的。其基本思想是：首先置一初始粒子作为种子，在远离种子的任意位置随机产生一个粒子使其做无规行走，直至与种子接触，成为集团的一部分；然后再随机产生一个粒子，重复上述过程，这样就可以得到足够大的DLA团簇。

初始种子可以是一个点，可以是一条线，也可以是其他边界。

Point attractor

Line attractor

Internal circular attractor

复变函数迭代

我们在文章开始就介绍了，迭代是重复反馈过程的活动，复变函数由于其特殊的性质，在其迭代计算 (例如：$Z_{k+1}=Z_k^m+$