自然界中线性总是特例,非线性才是自然界的常态。线性系统往往可以优雅的化整为零,一步一步从最简单的项叠加为最终结果。但是非线性系统往往无法直观的用某个函数去解析,而各种分岔、分形、混沌等行为,导致非线性系统更加难以被认知。
因此,随着人们这些年对非线性研究的发展,诞生出了很多非线性可视化方法,从繁琐的数学方程中解放出来,帮助人们直观的理解认知非线性系统的特性。在介绍常见的非线性动力系统中用的可视化方法前,先利用几个小引子,来直观的认识非线性的特征。
首先介绍一个研究迭代分形中,最经典的Julia集。
设置一个复数域上的函数f(z)=z^2+C。在初始值z0确定的情况下,可以通过迭代生成一些列的z值:
z0
z1=f(z0)
z2=f(z1)
z3=f(z2)...
对于不同的初始值z0,有的序列收敛,有的序列发散。
因此,我们根据迭代收敛的特性,把二维复数平面内每个点都代入,把收敛越快的点赋予越大的值,发散的点赋予最小值。这样就构成了Julia集。
下图展示了C=0.279的Julia集的可视化。颜色图提取了参考资料[2]中颜色。具体分为两步,第一步是计算出Julia集,第二步是为了可视化进行插值。
其中第一张图为没有插值的效果图。由于每个点的值都是整数,所以即使采用了shading interp,画出来的图仍然是台阶样式的。
因此,采用等高线类似的方法,提取等高线边缘的点,对数据再一次插值计算,得到下面的光滑连续的图:
Mandelbrot集的求解方法与Julia集方法类似,只是里面的C需要替换为每一个点的坐标z0,也就是:
f(z)=z^2+z0
Mandelbrot集的图像就像一个横躺着的大葫芦。绘制出的图形如下图:
我们可以看到,无论是Julia集还是Mandelbrot集,都存在非常多的微小结构。我们以Mandelbrot集为例,将(-0.912,-0.2611)这一点进行放大,可以得到:
可以看到将细节放大,还有更多的细节在等着我们。同时每一个细节似乎与周围的样子相似,存在规律,但又不同,无法精确的用数学语言描述。它的复杂度不会因为尺度的减小而消失。
从这个简单的迭代中,我们可以管中窥豹,见识非线性的复杂之处:能够感受到规律却无法用数学描述出来;处处存在细节无法获取全部信息;微小的差异也会导致最终结果的与众不同;然而背后的方程并不复杂。
最后放一个放大的动图在文章后面。微信上传有大小限制,所以画质动图质量可能不是太高。
代码附在后面:
%Mandelbrot
%?res是目标分辨率,iter是循环次数,(xc,yc)是图像中心,xoom是放大倍数
%res分辨率越多,在一张图上能够计算的点越多。iter越多,细节越丰富。
clear;
clc
close all
% res=1024;iter=100;xc=-0.748766710846959;yc=-0.123640847970064;xoom=6.503637847981780e+08;
% res=2048;iter=10000;xc=-0.912;yc=-0.2611;xoom=150;%为了生成图形好看,计算时间非常长
res=1024;iter=100;xc=-0.5;yc=0;xoom=1;
% res=2500;iter=200;xc=0;yc=0;xoom=1.7;
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;
C=z;
% C=0.279;
tic
%matlab官方推荐
N = ones(res,res);
%N=uint32(N);
for k = 0:iter
z = z.*z + C;
inside = abs( z )<=2;
%N = N + uint32(inside);
N = N + inside;
end
toc
%蓝色-橘色颜色条
mycolorpoint=[
[235 255 255];...
[243 253 242];...
[248 247 193];...
[255 211 49];...
[236 133 0];...
[134 36 21];...
[75 6 34];...
[15 5 64];...
[3 8 96];...
[7 35 154];...
[24 84 205];...
[87 169 229];...
[165 223 247];...
[235 255 255];...
];
mycolorposition=[1 3 7 14 27 45 57 66 74 86 96 109 119 128];
mycolormap_r=interp1(mycolorposition,mycolorpoint(:,1),1:128,'pchip','extrap');
mycolormap_g=interp1(mycolorposition,mycolorpoint(:,2),1:128,'pchip','extrap');
mycolormap_b=interp1(mycolorposition,mycolorpoint(:,3),1:128,'pchip','extrap');
mycolor=[mycolormap_r',mycolormap_g',mycolormap_b']/255;
figure()
pcolor(x,y,log(N));
mycolor=interp1(linspace(1,256,128),mycolor,1:256);
colormap([mycolor;mycolor;mycolor])
shading interp
Nq=LowPrecision2High(N);
figure()
pcolor(x,y,log(Nq));
colormap([mycolor;mycolor;mycolor])
shading interp
function Zq=LowPrecision2High(Z)
%对阶梯数据进行平滑
[Nx,Ny]=size(Z);
[X,Y]=meshgrid(1:Nx,1:Ny);
XP=[];YP=[];ZP=[];
for kx=1:(Nx-1)
for ky=1:(Ny-1)
if Z(kx,ky)~=Z(kx,ky+1)
x_t=0.5*(X(kx,ky)+X(kx,ky+1));
y_t=0.5*(Y(kx,ky)+Y(kx,ky+1));
z_t=0.5*(Z(kx,ky)+Z(kx,ky+1));
XP=[XP;x_t];YP=[YP;y_t];ZP=[ZP;z_t];
end
if Z(kx,ky)~=Z(kx+1,ky)
x_t=0.5*(X(kx,ky)+X(kx+1,ky));
y_t=0.5*(Y(kx,ky)+Y(kx+1,ky));
z_t=0.5*(Z(kx,ky)+Z(kx+1,ky));
XP=[XP;x_t];YP=[YP;y_t];ZP=[ZP;z_t];
%pause(0.1)
end
end
end
%边界值
for kx=1:Nx
XP=[XP;X(kx,1)]; YP=[YP;Y(kx,1)]; ZP=[ZP;Z(kx,1)];
XP=[XP;X(kx,end)];YP=[YP;Y(kx,end)];ZP=[ZP;Z(kx,end)];
end
for ky=2:Ny-1
XP=[XP;X(1,ky)]; YP=[YP;Y(1,ky)]; ZP=[ZP;Z(1,ky)];
XP=[XP;X(end,ky)];YP=[YP;Y(end,ky)];ZP=[ZP;Z(end,ky)];
end
Zq=griddata(XP,YP,ZP,X,Y,'natural');%'v4'最好,实在不行用linear
end
参考资料:
[1] ww2.mathworks.cn/help/gpucoder/gs/gpu-code-generation-the-mandelbrot-set.html#responsive_offcanvas
[2] www.matrix67.com/blog/archives/292