快捷導(dǎo)航

基于Matlab實(shí)現(xiàn)離散系統(tǒng)分岔圖的繪制

更新時(shí)間：2022年04月27日 15:40:17 作者：slandarer

這篇文章主要介紹了如何利用Matlab實(shí)現(xiàn)離散分岔圖的繪制，文中的示例代碼講解詳細(xì)，對(duì)我們學(xué)習(xí)Matlab有一定的幫助，需要的可以參考一下

1.一維離散分岔圖

一維那非常簡(jiǎn)單哈，就循環(huán)著畫唄，以下舉兩個(gè)簡(jiǎn)單的例子：

% 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 % 穩(wěn)定后開始繪圖 
        plot(r,x,'k.','MarkerSize',1);
        drawnow
    end
end

% Logistic系統(tǒng)
% 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 % 穩(wěn)定后開始繪圖 
        plot(r,x,'k.','MarkerSize',1);
        drawnow
    end
end

橫坐標(biāo)代表參數(shù)的數(shù)值，縱坐標(biāo)表示該參數(shù)數(shù)值下序列可能的取值，n>100再開始畫圖是為了讓序列通過迭代穩(wěn)定下來，事實(shí)上我么可以不設(shè)置n>100，同時(shí)將顏色設(shè)置為隨著n變化的漸變色，可以發(fā)現(xiàn)幾乎看不出漸變來，該序列穩(wěn)定的很快(以下是繪圖部分代碼的微調(diào))：

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

當(dāng)然我們可以設(shè)置n為奇數(shù)和偶數(shù)時(shí)繪制不同顏色，下圖所示，對(duì)于該系統(tǒng)而言，其序列的數(shù)值是反復(fù)橫跳的(以下是繪圖部分代碼的微調(diào))：

當(dāng)然可以設(shè)置更多顏色：

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

2.二維離散分岔圖

繪制Henon系統(tǒng)的分岔圖：

定住b值不變，改變a值，觀察y序列，不同b值時(shí)繪制效果不同：

% x(n+1)=1+y(n)-a*x(n)^2
% y(n+1)=b*x(n)
% Henon系統(tǒng)
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 % 穩(wěn)定后開始繪圖 
        plot(a,y,'k.','MarkerSize',1);
        drawnow
    end
end

b=0.2時(shí)繪制效果

b=0.3時(shí)繪制效果

3.封面圖繪制

經(jīng)典體現(xiàn)理科生工科生藝術(shù)情懷環(huán)節(jié)，我們?cè)趺茨軌驅(qū)⒎植韴D的美忽視？

感覺大家很多也是因?yàn)榭捶饷鎴D點(diǎn)進(jìn)來的，雖然不短，但還是把代碼放一下叭，原理很簡(jiǎn)單，構(gòu)造一個(gè)矩陣統(tǒng)計(jì)各個(gè)位置點(diǎn)數(shù)量，然后依據(jù)點(diǎn)數(shù)量映射到顏色：

圖一

% x(n+1)=1+y(n)-a*x(n)^2
% y(n+1)=b*x(n)
% Henon系統(tǒng)

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(['進(jìn)度:[',num2str(n),'/12000]']);
    ty=round((y+0.4)*500);
    ta=a*500;
    index=round((ta).*451+ty);
    pntMat(index)=pntMat(index)+1;
end

% 矩陣上下翻轉(zhuǎn)(坐標(biāo)y軸方向與圖片序數(shù)相反)
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(['進(jìn)度:[',num2str(n),'/20000]']);
    if n>1 % 穩(wěn)定后開始繪圖     
        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);

% 矩陣上下翻轉(zhuǎn)(坐標(biāo)y軸方向與圖片序數(shù)相反)
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)