基于Matlab实现离散系统分岔图的绘制基于Matlab实现离散系统分岔图

1.一维离散分岔图
2.二维离散分岔图
3.封面图绘制

1.一维离散分岔图一维那非常简单哈，就循环着画呗，以下举两个简单的例子：

文章图片

% x(n+1)=1-r*x(n)^2% (r∈(0,2),x∈[-1,1])的分支混沌图。hold onf=@(x,r)1-r.*x.^2; r=0:.01:2; x=0; % x初值for n=1:1000 x=f(x,r); if n>100 % 稳定后开始绘图 plot(r,x,'k.','MarkerSize',1); drawnowendend

文章图片

% Logistic系统% x(n+1)=r*x(n)-r*x(n)^2% (r∈(2.6,4),x∈(0,1])的分支混沌图。hold onf=@(x,r)r.*x-r.*x.^2; r=2.6:.01:4; x=0.6; % x初值for n=1:1000 x=f(x,r); if n>100 % 稳定后开始绘图 plot(r,x,'k.','MarkerSize',1); drawnowendend

文章图片

横坐标代表参数的数值，纵坐标表示该参数数值下序列可能的取值，n>100再开始画图是为了让序列通过迭代稳定下来，事实上我么可以不设置n>100，同时将颜色设置为随着n变化的渐变色，可以发现几乎看不出渐变来，该序列稳定的很快(以下是绘图部分代码的微调)：

c1=[0 0.4470 0.7410]; c2=[0.6350 0.0780 0.1840]; N=1000; for n=1:Nx=f(x,r); plot(r,x,'.','Color',(n.*c1+(N-n).*c2)./N,'MarkerSize',2); drawnowend

文章图片

当然我们可以设置n为奇数和偶数时绘制不同颜色，下图所示，对于该系统而言，其序列的数值是反复横跳的(以下是绘图部分代码的微调)：

文章图片

当然可以设置更多颜色：

for n=1:1000x=f(x,r); switch mod(n,4)case 3,plot(r,x,'.','Color',[0.4660 0.6740 0.1880],'MarkerSize',2); case 2,plot(r,x,'.','Color',[0.8500 0.3250 0.0980],'MarkerSize',2); case 1,plot(r,x,'.','Color',[0 0.4470 0.7410],'MarkerSize',2); case 0,plot(r,x,'.','Color',[0.6350 0.0780 0.1840],'MarkerSize',2); enddrawnowend

文章图片

2.二维离散分岔图绘制Henon系统的分岔图：

文章图片

定住b值不变，改变a值，观察y序列，不同b值时绘制效果不同：

% x(n+1)=1+y(n)-a*x(n)^2% y(n+1)=b*x(n)% Henon系统hold onfx=@(x,y,a)1+y-a.*x.^2; fy=@(x,b)b.*x; a=0:.002:1.4; b=0.2; x=0; y=0; for n=1:800 lx=x; ly=y; x=fx(lx,ly,a); y=fy(lx,b); if n>100 % 稳定后开始绘图 plot(a,y,'k.','MarkerSize',1); drawnowendend

b=0.2时绘制效果

文章图片

b=0.3时绘制效果

文章图片

3.封面图绘制经典体现理科生工科生艺术情怀环节，我们怎么能够将分岔图的美忽视？
感觉大家很多也是因为看封面图点进来的，虽然不短，但还是把代码放一下叭，原理很简单，构造一个矩阵统计各个位置点数量，然后依据点数量映射到颜色：
图一

% x(n+1)=1+y(n)-a*x(n)^2% y(n+1)=b*x(n)% Henon系统fx=@(x,y,a)1+y-a.*x.^2; fy=@(x,b)b.*x; a=0:.002:1.4; b=0.3; x=0; y=0; % 填充矩阵pntMat=zeros(451,701); for n=1:12000lx=x; ly=y; x=fx(lx,ly,a); y=fy(lx,b); disp(['进度:[',num2str(n),'/12000]']); ty=round((y+0.4)*500); ta=a*500; index=round((ta).*451+ty); pntMat(index)=pntMat(index)+1; end% 矩阵上下翻转(坐标y轴方向与图片序数相反)pntMat=flipud(pntMat); % 绘图imagesc(pntMat); caxis([0,50])ax=gca; hold on; ax.XTick=[]; ax.YTick=[]; % 颜色映射map=[0.1294 0.0549 0.1725; 0.2196 0.1608 0.2902; 0.3882 0.1804 0.4941; 0.4392 0.1922 0.4706; 0.5333 0.2235 0.4392; 0.6471 0.2588 0.3686; 0.7137 0.2745 0.3294; 0.7725 0.3059 0.2902; 0.8510 0.3725 0.2275; 0.9137 0.4196 0.1804; 0.9608 0.5020 0.2000; 0.9765 0.5529 0.2078; 0.9804 0.6431 0.2549; 0.9843 0.6627 0.2706; 0.9765 0.7176 0.3412; 0.9765 0.7686 0.4000; 0.9765 0.8118 0.4902; 0.9725 0.8510 0.5961; 0.9882 0.9020 0.6667; 1.0000 0.9451 0.8431; 1.0000 0.9961 0.9804; 1.0000 1.0000 1.0000]; Xi=1:size(map,1); Xq=linspace(1,size(map,1),800); map=[interp1(Xi,map(:,1),Xq,'linear')',...interp1(Xi,map(:,2),Xq,'linear')',...interp1(Xi,map(:,3),Xq,'linear')']; colormap(map)

文章图片

图二

% x(n+1)=1-r*x(n)^2% (r∈(0,2),x∈[-1,1])的分支混沌图。f=@(x,r)1-r.*x.^2; r=0:.0025:2; x=0; % x初值pntMat=zeros(801,801); for n=1:20000 x=f(x,r); disp(['进度:[',num2str(n),'/20000]']); if n>1 % 稳定后开始绘图tx=round((x+1)*400); tr=r*400; index=round((tr).*801+tx); pntMat(index)=pntMat(index)+1; endend% 为了减少锯齿化高斯模糊一下pntMat=imgaussfilt(pntMat,0.3); % 矩阵上下翻转(坐标y轴方向与图片序数相反)pntMat=flipud(pntMat); % 绘图imagesc(pntMat); caxis([0,80])ax=gca; hold on; ax.XTick=[]; ax.YTick=[]; % 颜色映射map=[0.14000.11000.15000.28000.09000.41000.27000.21000.51000.23000.32000.55000.19000.42000.56000.15000.51000.56000.12000.58000.55000.14000.67000.51000.24000.73000.46000.39000.80000.37000.59000.85000.2500]; Xi=1:size(map,1); Xq=linspace(1,size(map,1),800); map=[interp1(Xi,map(:,1),Xq,'linear')',...interp1(Xi,map(:,2),Xq,'linear')',...interp1(Xi,map(:,3),Xq,'linear')']; colormap(map)