MATLAB非线性可视化（引2）牛顿迭代分形

接着上次的Mandelbrot集，这次再介绍一个牛顿迭代中的非线性现象。

牛顿迭代法是一种非常简单的求解根的方法，利用该点处导数的信息，通过每一次的迭代，使得点逐渐向解靠近。

牛顿迭代法的公式为：

我们以复数平面中，简单的二次方程为例：

在[-2,2]区间内，绘制出每一个点牛顿迭代过程的轨迹，如下图：

可以看到，方程的根只有x=1和x=-1两个，在短短几步之内，整个平面的点都可以快速收敛到这两个根。

其中左半边(实部小于零)的点全部收敛到了x=-1这个根；右半边的点，全部收敛到x=1这个根。这看上去很简单。

然而当阶数上升，变为p阶多项式的时候，非线性就开始出现了。我们以3次方程为例，f(x)=x^3-1。三个根分别为：

x=1

x=-0.5+0.866i

x=-0.5-0.866i

我们继续观察平面内点的收敛过程，如下图所示：

可以看到这些点并不是直接向这三个根靠拢，有很多点被这三个根来回吸引，很久之后才逐渐收敛。

我们把整个平面内，初始点x0对应的最终根xi进行绘制，可以得到下面这个图：

蓝色对应的第一个根(x=1)，黄色对应着第二个根(x=-0.5+0.866i)，绿色对应着第三个根(x=-0.5-0.866i)。在根的附近处，所有点都可以直接收敛到对应的根。然而在每个根的交接处，受不同根的互相吸引的影响，交界面处出现了分形现象。

此后，当继续增加多项式的阶次，会诞生出越来越复杂的分形图案。比如下图为p=6时的图案

混沌中，3是个很神奇的数字，很多系统当数量大于3时，就会出现混沌。正所谓一生二，二生三，三生万物，这也是科学界所说的more is different的体现。

对于其他大于三阶的多项式，大家可以自己在后面的程序中简单更改查看。而且不仅仅限于多项式，只要有比较多的解，对于其它形式的方程，也可以观察到分形现象。

下图是x^2*sin(x)-1=0方程的牛顿迭代根的分布图：

后面附上程序：

%牛顿迭代法程序
clear;
close all
res=1024;xc=0;yc=0;

xoom=1;
x0=xc-2/xoom;x1=xc+2/xoom;
y0=yc-2/xoom;y1=yc+2/xoom;
x=linspace(x0,x1,res);
y=linspace(y0,y1,res);
[xx,yy]=meshgrid(x,y);
z=xx+yy*1i;

p=3;%3次多项式

R = ones(res,res);
for k=1:numel(R)
    zk=z(k);
    dz=1;
    while abs(dz)>1e-3
         dz=(zk^p-1)/(p*zk^(p-1));%f(x)/f'(x)的形式
%          dz=(zk^2*sin(zk)-1)/(zk^2*cos(zk)+2*zk*sin(zk));%x*x*sinx的方程.改p=24
        zk = zk-dz;
    end
    R(k)=zk;
end

Z_real=real(R);
Z_imag=imag(R);
idx=kmeans([Z_real(:),Z_imag(:)],p);

RootNum=zeros(res,res);
RootNum(:)=idx(:);

figure()
pcolor(x,y,RootNum);shading flat
% colormap(gallery('uniformdata',[p,3],13))%改p=24%x*x*sinx的方程演示用
colormap(gallery('uniformdata',[p,3],7))

参考资料：

[1]俞碧川. 基于牛顿迭代图形的丝绸提花织物纹理设计方法[D]. 浙江理工大学, 2011.