一维那非常简单哈，就循环着画呗，以下举两个简单的例子 ：
% x(n+1)=1-r*x(n)^2
% (r∈(0,2),x∈[-1,1])的分支混沌图。
hold on
f=@(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);
drawnow
end
end
% Logistic系统
% x(n+1)=r*x(n)-r*x(n)^2
% (r∈(2.6,4),x∈(0,1])的分支混沌图。
hold on
f=@(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);
drawnow
end
end
横坐标代表参数的数值，纵坐标表示该参数数值下序列可能的取值，n>100再开始画图是为了让序列通过迭代稳定下来，事实上我么可以不设置n>100，同时将颜色设置为随着n变化的渐变色，可以发现几乎看不出渐变来，该序列稳定的很快(以下是绘图部分代码的微调)：
c1=[0 0.4470 0.7410];
c2=[0.6350 0.0780 0.1840];
N=1000;
for n=1:N
x=f(x,r);
plot(r,x,'.','Color',(n.*c1+(N-n).*c2)./N,'MarkerSize',2);
drawnow
end
当然我们可以设置n为奇数和偶数时绘制不同颜色，下图所示，对于该系统而言，其序列的数值是反复横跳的(以下是绘图部分代码的微调)：
当然可以设置更多颜色：
for n=1:1000
x=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);
end
drawnow
end
绘制Henon系统的分岔图：
定住b值不变，改变a值，观察y序列，不同b值时绘制效果不同：
% x(n+1)=1+y(n)-a*x(n)^2
% y(n+1)=b*x(n)
% Henon系统
hold on
fx=@(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);
drawnow
end
endb=0.2时绘制效果
b=0.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:12000
lx=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;
end
end
% 为了减少锯齿化高斯模糊一下
pntMat=imgaussfilt(pntMat,0.3);
% 矩阵上下翻转(坐标y轴方向与图片序数相反)
pntMat=flipud(pntMat);
% 绘图
imagesc(pntMat);
caxis([0,80])
ax=gca;
hold on;
ax.XTick=[];
ax.YTick=[];
% 颜色映射
map=[0.1400 0.1100 0.1500
0.2800 0.0900 0.4100
0.2700 0.2100 0.5100
0.2300 0.3200 0.5500
0.1900 0.4200 0.5600
0.1500 0.5100 0.5600
0.1200 0.5800 0.5500
0.1400 0.6700 0.5100
0.2400 0.7300 0.4600
0.3900 0.8000 0.3700
0.5900 0.8500 0.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)
