首页 > 代码库 > matlab 混沌，分形
matlab 混沌，分形
对于函数f(x)=λsin(πx)，λ∈(0,1]，使用matlab计算随着λ逐渐增大，迭代x=f(x)的值，代码如下：
function y=diedai(f,a,x1)
N=32;
y=zeros(N,1);
for i=1:1e4
x2=f(a,x1);
x1=x2;
y(mod(i,N)+1)=x2;
end
end
%f=@(a,x)a*x*(1-x);
f=@(a,x)a*sin(pi*x);
%x0=0.1;
hold on;
for x0=-1:0.05:1
for a=0:0.01:1
y=diedai(f,a,x0);
for count=1:32
plot(a,y(count),‘k.‘);
hold on;
end
end
end
得到的图像如下：其中横轴为λ，纵轴为x
可以看到随着λ的逐渐增大，出现了倍周期分叉的情况。
由图中可以看出第一个分叉值大约在0.3附近，第二个在0.73到0.75之间，第三个在0.8到0.85之间，混沌大约出现在0.86附近。接下来编写代码计算分叉值，代码如下：
format long;
x0=0.1;
for a=0.3182:0.0000001:0.3183
y=diedai(f,a,x0);
if max(y)>0.001
disp(a);
break;
end
end
得到第一个分叉值大约为0.3182298
format long;
x0=0.1;
for a=0.7199:0.000001:0.72
y=diedai(f,a,x0);
if max(y)-min(y)>0.001
disp(a);
break;
end
end
得到第二个分叉值大约为0.719911
format long;
x0=0.1;
for a=0.8332:0.000001:0.8333
y=diedai(f,a,x0);
if abs(y(32)-y(30))>0.001
disp(a);
break;
end
end
得到第三个分叉值大约为0.833267
利用Feigenbaum常数估计第三个分叉值，得到0.805939
分形图
周常青
画mandelbrot分形图，主要使用了三个函数：iter=mandelbrot1(x0,y0,maxIter)，用来计算迭代后是否收敛，方程z=z2+z0。c=color(iter,maxIter)计算颜色值，返回[r g b]。draw_mandelbrot1用来绘制图像。
function iter=mandelbrot1(x0,y0,maxIter)
x=x0;
y=y0;
for i=1:1:maxIter
if (x*x+y*y)>=4
iter=i;
break;
else
tem=x*x-y*y+x0;
y=x*y*2+y0;
x=tem;
end
iter=i;
end
end
--------------------------------------------------------------------------------------------------------
function c=color(iter,maxIter)
if iter==maxIter
c=[1,0,0];
else
c1=abs(mod((iter*20+255),510)-255);
c2=abs(mod((iter*15+85+255),510)-255);
c3=abs(mod((iter*30+171+255),510)-255);
c=[c1/255,c2/255,c3/255];
end
end
-------------------------------------------------------------------------------------------------------
function draw_mandelbrot1
for y=0:ymax
for x=0:xmax
xt=2*r*x/xmax+x0-r;
yr=r*ymax/xmax;
yt=2*yr*y/ymax+y0-yr;
iter=mandelbrot1(xt,yt,maxIter);
%c=color(iter,maxIter);
c=color(iter,maxIter);
plot(xt,yt,‘.‘,‘color‘,c);
hold on;
end
end
end
输入draw_mandelbrot1(-0.5,300,0,200,2,300)，得到的图像如下：
对收敛部分的颜色也进行绘制，画julia分形图，主要使用了三个函数[xList,yList,iter]=julia1(x0,y0,maxIter,jx0,jy0)，c=color2(iter,maxIter,xList,yList)，draw_julia1(x0,xmax ,y0,ymax, r,maxIter,jx0,jy0)。其中julia1使用方程z’=z2+zm进行迭代，z0=x0+j*y0,zm=jx0+j*jy0;返回iter为迭代次数，color2根据不同的iter，和x，y值进行差值，计算出c=[r g b]，draw_julia1进行绘图。
在matlab里输入draw_julia1(0,300,0,200,1.6,50,-0.78888,0.212325)，得到如下图
代码如下：
function [xList,yList,iter]=julia1(x0,y0,maxIter,jx0,jy0)
xList(1)=x0;
yList(1)=y0;
M=256;
lnln_m=log(abs(log(M)));
xbck=0;
ybck=0;
x=x0;
y=y0;
for i=1:maxIter
if (x*x+y*y>=M)
% iter=i;
break;
end
xbck=x;
ybck=y;
tem=x*x-y*y+jx0;
y=x*y*2+jy0;
x=tem;
xList(i+1)=x;
yList(i+1)=y;
end
if i~=maxIter
lnln_z=log(abs(log(x*x+y*y)));
lnln_zbak=log(abs(log(xbck*xbck+ybck*ybck)));
iter=i-2-(lnln_z-lnln_m)/(lnln_z-lnln_zbak);
else
iter=i;
end
end
----------------------------------------------------------------------------------------------------
function c=color2(iter,maxIter,xList,yList)
% xList=xyList(1:length(xyList)/2);
% yList=xyList(length(xyList)/2:length(xyList));
sincolorf=@(x)(sin(x*2*pi/510-pi*0.5)+1)*0.5*255;
if iter==maxIter
x=xList(maxIter);
y=yList(maxIter);
z=sqrt(x*x+y*y);
zd=z-sqrt(xList(maxIter-1)*xList(maxIter-1)+yList(maxIter-1)*yList(maxIter-1));
r1=sincolorf(z*2000);
g1=sincolorf(y*x*1000);
b1=sincolorf(zd*1000);
% errcolor(1)=errcolor(1)+r1;
% errcolor(2)=errcolor(2)+g1;
% errcolor(3)=errcolor(3)+b1;
else
r1=sincolorf(iter*20);
g1=sincolorf(iter*15+85);
b1=sincolorf(iter*30+171);
% errcolor(1)=errcolor(1)+r1;
% errcolor(2)=errcolor(2)+g1;
% errcolor(3)=errcolor(3)+b1;
end
% result_r=fTcolor(errcolor(1));
% result_b=fTcolor(errcolor(2));
% result_g=fTcolor(errcolor(3));
r=r1/255;
g=g1/255;
b=b1/255;
c=[r g b];
end
-----------------------------------------------------------------------------------------------------------------------------
function draw_julia1(x0,xmax ,y0,ymax, r,maxIter,jx0,jy0)
for y=1:ymax
for x=1:xmax
xt=2*r*x/xmax+x0-r;
yr=r*ymax/xmax;
yt=2*yr*y/ymax+y0-yr;
[xList,yList,iter]=julia1(xt,yt,maxIter,jx0,jy0);
% xyList=[xList yList];
c=color2(iter,maxIter,xList,yList);
plot(xt,yt,‘.‘,‘color‘,c);
hold on;
%drawnow;
end
end
end
可以看到，通过color2对收敛区域内部的颜色进行计算，可以得到色彩更为丰富的分形图，但是总体来说，因为对于相关知识的欠缺，比较难找到合适的颜色方案。
matlab 混沌，分形