利用matlab实现复数域空间牛顿迭代法的分形图案展示（newton法）

本文首发于 matlab爱好者 微信公众号，欢迎关注。

惯例声明：本人没有相关的工程应用经验，只是纯粹对相关算法感兴趣才写此博客。所以如果有错误，欢迎在评论区指正，不胜感激。本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。

本文为学习非线性动力学时，搜索资料时发现的一个小例子，感觉很有趣，试着自己尝试一下。相关实现案例非常多，所以这篇文章就简单水一下了。

1 一维函数的牛顿迭代法

牛顿迭代法是一种经典的求解函数根的方法，它是一种迭代方法。其思想为沿着当地的切线方向不断移动计算点，使其不断逼近真实根值。

这种方法在导数比较简单，或者初始点距离根足够近的时候，计算结果还是比较稳定的。

假设初始值为xk，函数为f(x)
$$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)}$$
其下一步预测的根就是$x_{k+1}$。然后逐次迭代，x值会越来越接近真实的根值。

下面以y=x^2-2为例。可以看到每次迭代，各个点逐渐向根处靠近的过程。由于该函数的特点，所有点第一次迭代后都跑到了根的外侧，然后再逐渐向根靠近。

clear
clc
close all

x=-2.5:0.01:2.5;
y=Fx(x);%绘制函数用

xp=[-2.5:0.1:2.5];%绘制点用

figure(1);
clf
hold on
plot(x,y,'color','k','linewidth',1);
plot(x,0*x,'color','k','linewidth',0.5);
yp2=Fx(xp);
scatter(xp,yp2,30,1:length(xp),'filled');%绘制点
hold off
xlim([min(x),max(x)])
F=getframe(gcf);
I=frame2im(F);
[I,map]=rgb2ind(I,256);
%保存动图
imwrite(I,map,'test.gif','gif','Loopcount',inf,'DelayTime',0.8);

for k=2:6
    fxp=Fx(xp);
    dfxp=dFx(xp);%导数
    xp2=xp-fxp./dfxp;%经过牛顿迭代一次
    figure(1);
    clf
    hold on
    plot(x,y,'color','k','linewidth',1);
    plot(x,0*x,'color','k','linewidth',0.5);
    yp2=Fx(xp2);
    scatter(xp2,yp2,30,1:length(xp),'filled');%绘制点
    hold off
    pause(0.1)
    xp=xp2;%更新点的位置
    xlim([min(x),max(x)])
    %保存动图
    F=getframe(gcf);
    I=frame2im(F);
    [I,map]=rgb2ind(I,256);
    imwrite(I,map,'test.gif','gif','WriteMode','append','DelayTime',0.8);
end
function y=Fx(x)
y=x.^2-2;
end
function y=dFx(x)
%导数
y=2*x;
end

2 复平面的牛顿迭代法

复平面的牛顿迭代方法和普通的一维迭代方法相同，公式都一样，就是把原本实数x替换为复数x。其余不变。
$$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)}$$

2.1 简单方程结果

我们以复平面中简单的二次方程为例：
$$f(x)=x^2-1$$
它的牛顿迭代求根的迭代方程为：
$$x_{k+1}=x_k-\frac{x_k^2-1}{2x_k}$$
当然，上面的数都是复数。

我们绘制出其在实部和虚部都在[-2,2]区间内，每一个点牛顿迭代过程的轨迹，如下图：

可以看到，方程的根只有 x=1+0i 和 x=-1+0i 这两个。而且在比较短的迭代次数里，整个平面的点都可以快速收敛到这两个根。

如果仔细观察，还可以看到左半边（实部小于0）的点都收敛到了x=-1这个根上，而右半边的点都收敛到了x=1这个根上。

对于$f=x^2-1$，还没有非线性的特性。现在保持$f=x^n-1$形式不变，当n>3时，再进行牛顿迭代，就会产生非线性的收敛结果。

下面以$f=x^3-1$为例，
它的牛顿迭代求根的迭代方程为：
$$x_{k+1}=x_k-\frac{x_k^3-1}{3x_k}$$
该方程有3个根，分别为x=1,x=-0.5+0.866i,x=-0.5-0.866i其在实部和虚部都在[-2,2]区间内，每一个点牛顿迭代过程的轨迹，如下图：
可以看到相比较k=2的结果，这些点收敛速度比较慢，而且并不是直接向这3个根靠拢。·有些点被三个根来回吸引，很久之后才能收敛。

%绘制牛顿迭代法的详解
clear;
clc
close all
res=100;
x=linspace(-2,2,res);
y=linspace(-2,2,res);
[xx,yy]=meshgrid(x,y);
zz=xx+yy*1i;
p=3;%3次多项式
N=50;
R=zeros(res,res,N);
zk=zz;
for k=1:N
    R(:,:,k)=zk;
    dz=(zk.^p-1)./(p*zk.^(p-1));%f(x)/f'(x)的形式
    zk = zk-dz;
end
%指定每个点的颜色
Color_R=interp2([-2,2;-2,2],[-2,-2;2,2],[0,0;1,1],xx,yy);
Color_G=interp2([-2,2;-2,2],[-2,-2;2,2],[1,0.2;0.3,0],xx,yy);
Color_B=interp2([-2,2;-2,2],[-2,-2;2,2],[0,1;0,1],xx,yy);
%动画中，每个帧之间插值补充m个帧
figure(1)
m=9;
count=0;
for k=1:9
    R_k_real_1=real(R(:,:,k));
    R_k_imag_1=imag(R(:,:,k));
    R_k_real_2=real(R(:,:,k+1));
    R_k_imag_2=imag(R(:,:,k+1));
    dt=1/(m+1);
    for j=0:1:m
        t=dt*j;
        R_k_real=R_k_real_1*(1-t)+R_k_real_2*t;
        R_k_imag=R_k_imag_1*(1-t)+R_k_imag_2*t;
        clf
        scatter( R_k_real(:) , R_k_imag(:) ,6,[Color_R(:),Color_G(:),Color_B(:)],'filled')
        xlim([-2,2])
        ylim([-2,2])
        box on
        pause(0.1)
        count=count+1;
        F=getframe(gca);
        I=frame2im(F);
        [I,map]=rgb2ind(I,192);
        if count == 1
            imwrite(I,map,'test.gif','gif','Loopcount',inf,'DelayTime',0.9);
        else
            imwrite(I,map,'test.gif','gif','WriteMode','append','DelayTime',0.1);
        end
    end
end

还是这个n=3函数的牛顿迭代。

如果我们把最终收敛到x=1的根的所有点标记为蓝色，把所有最终收敛到x=-0.5+0.866i的点标记为黄色，所有收敛到x=-0.5-0.866i的点标记为绿色，那么在空间中的分布为下面的图案：

同理，更改n的阶数，当n=6的时候，会形成下面这个更为复杂的图案。

正所谓二生三，三生万物，混沌现象与数字3有着不解之缘。

另外注：当n=4时，由于根在像素的对角线上，所以在对角线上的像素点甚至会出现不收敛的情况。建议这时候把方程稍微加一点扰动，变为$f=x^4-1-10^{-3}\ast i$。其它涉及像素刚好在根连线处的（比如水平线垂直线或对角线），处理方式相同。
代码见下一小节。

2.2 其它非线性方程结果

由于我上面方程挑的方程为x^n-1的形式，所以根的分布都呈现出中心对称的图案。

如果变化方程，还会产生很多新的图案。下面就举几个简单的例子：

我们取$f=0.5\ast x^3-1.2\ast x^2+0.5\ast x+0.9$为例，其牛顿迭代公式如下：
$$x_{k+1}=x_k-\frac{0.5x_k^3-1.2x_k^2+0.5x_k+0.9}{1.5x_k^2-2.4x_k+x_k}$$
生成的图片如下：

可以看到其基本框架与标准的三阶相似，因为都是由3个根主导。但是生成图案的细节结果却大不相同。对于更高阶图案，如果更改参数，可以产生更为丰富的结果。

当然不限于多项式，其它函数也可以，只要能够用复数表达即可。比如：
$$f=x^2\sin x$$
其牛顿迭代公式如下：
$$x_{k+1}=x_k-\frac{x^2\sin x}{x^2\cos x+2x\sin x}$$
由于其复数根比较多，因此只取其一部分根进行染色，效果如下：

代码如下：

%Mandelbrot
%?res是目标分辨率，（xc,yc）是图像中心，xoom是放大倍数
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
        %在这里更改牛顿迭代的f(x)/f'(x)形式
        %dz=(zk^p-1)/(p*zk^(p-1));%f(x)/f'(x)的形式
        %dz=(zk^p-1-1e-3*1i)/(p*zk^(p-1));%f(x)/f'(x)的形式，这里是n=4以及其它情况，加一个1e-3的小扰动
        dz=(0.5*zk^3-1.2*zk^2+0.5*zk+0.9)/(1.5*zk^2-2.4*zk+0.5);%x^3+x^2+x+1的方程,对应p=3
        %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用，这里只要好看就行，没有特别要求
% colormap(gallery('uniformdata',[p,3],13))%改p=24%x*x*sinx的方程演示用
colormap(gallery('uniformdata',[p,3],7))