基于Matlab繪制洛倫茲吸引子相圖

更新時間：2022年04月12日 10:41:58 作者：slandarer

洛倫茲吸引子(Lorenz attractor)是由MIT大學(xué)的氣象學(xué)家Edward Lorenz在1963年給出的。本文將利用Matlab實現(xiàn)洛倫茲吸引子相圖的繪制，感興趣的可以了解一下

1.公式及Lorenz函數(shù)

Lorenz微分方程組定義如下：

非常容易能寫出該微分方程組函數(shù)：

function dL=Lorenz(t,L)
% L=[x;y;z;a;r;b];
% dL=[dx/dt;dy/dt;dz/dt;0,0,0];
% dz/dt=-a*(x-y)
% dy/dt=x*(r-z)-y
% dz/dt=x*y-b*z
dL=zeros([6,1]);
dL(1)=-L(4)*(L(1)-L(2));
dL(2)=L(1)*(L(5)-L(3))-L(2);
dL(3)=L(1)*L(2)-L(6)*L(3);
dL(4:6)=0;
end

2.混沌吸引子圖像

基礎(chǔ)繪圖：

[~,L]=ode45(@(t,L)Lorenz(t,L),0:.01:100,[1;1;1;10;28;8/3]); 
plot3(L(:,1),L(:,2),L(:,3))
grid on

修飾動態(tài)圖：

% ode45求解
[~,L]=ode45(@(t,L)Lorenz(t,L),0:.01:100,[1;1;1;10;28;8/3]); 
% 修飾及屬性設(shè)置
ax=gca;
hold on;grid on
plhdl=plot3(0,0,0,'Color',[0.9843 0.8588 0.5333 0.5],'LineWidth',1.3);
ax.XColor=[1,1,1].*.6;ax.XLim=[-20,20];
ax.YColor=[1,1,1].*.6;ax.YLim=[-30,30];
ax.ZColor=[1,1,1].*.6;ax.ZLim=[0,50];
ax.LineWidth=1.5;
ax.GridAlpha=.09;
ax.GridLineStyle='-.';
ax.FontName='cambria';
ax.Color=[0 0 0];
ax.DataAspectRatio=[1,1,1];
view([-159,18]);
% 循環(huán)繪圖
for i=1:size(L,1)
    plhdl.XData=L(1:i,1);
    plhdl.YData=L(1:i,2);
    plhdl.ZData=L(1:i,3);
    drawnow
end

3.混沌吸引子圖像

基本代碼：

這里使用龐加萊截面法，即繪制y=x平面上|y|的圖像，基本代碼如下：

Z=[]; 
for r=1:500
    % 舍棄前面迭帶的結(jié)果，用后面的結(jié)果畫圖 
    [~,L]=ode45(@(t,L)Lorenz(t,L),[0,1],[1;1;1;10;r;8/3]);  
    [T,L]=ode45(@(t,L)Lorenz(t,L),[0,50],L(end,:)); 
    D=L(:,2)-L(:,1); 
    
    for k2=2:size(L,1)
        k1=k2-1;
        if D(k1)*D(k2)<=0
            y=(L(k2,1).*L(k1,2)-L(k1,1).*L(k2,2))./(D(k2)-D(k1));
            Z=[Z,r+abs(y').*1i];
        end  
    end
    
end
plot(Z,'.','markersize',1)
title('Lorenz映射分岔圖')
xlabel('r'),ylabel('|y| where x=y')

代碼有一些地方詳細(xì)講解一下，首先說明為什么要用

Z=[Z,r+abs(y').*1i];

的格式進(jìn)行存儲，這樣存儲可以少構(gòu)造一個數(shù)組，一般情況下我們需要分別存儲γ和|y|到兩個矩陣，存儲為復(fù)數(shù)形式就可以復(fù)平面繪圖減少初始化矩陣數(shù)量。

其次代碼中用了D(k1)*D(k2)<=0來判斷是否采點，

D(k1)=x1-y1,D(k2)=x2-y2

當(dāng)D(k1)*D(k2)<=0時說明(x1,y1),(x2,y2)兩點分別在 y=x直線兩側(cè)。

另外說明一下：

y=(L(k2,1).*L(k1,2)-L(k1,1).*L(k2,2))./(D(k2)-D(k1));

是啥。

其實就是構(gòu)造的兩點連線與直線y=x的交點：

PS:為了進(jìn)一步減少空間復(fù)雜度，我們可以將上述函數(shù)更改為完全由x,y差值以及y代替，這樣就可以直接將中間變量D存儲到原來x的位置，減少中間變量的數(shù)量：

因此代碼可以改寫為（當(dāng)然為了可讀性最后并沒有采取這個策略hiahiahia）：

L(:,1)=L(:,2)-L(:,1); 
    
    for k2=2:size(L,1)
        k1=k2-1;
        if L(k1,1)*L(k2,1)<=0
            y=L(k2,2)+(L(k1,2)-L(k2,2)).*L(k2,1)./(L(k2,1)-L(k1,1));
            Z=[Z,r+abs(y').*1i];
        end  
    end

最后，這部分代碼依賴循環(huán)我們完全可以將其向量化，即修改為：

Z=[]; 
for r=1:500
    % 舍棄前面迭帶的結(jié)果，用后面的結(jié)果畫圖 
    [~,L]=ode45(@(t,L)Lorenz(t,L),[0,1],[1;1;1;10;r;8/3]);  
    [T,L]=ode45(@(t,L)Lorenz(t,L),[0,50],L(end,:)); 
    
    % 找到穿過直線y=x的前后兩個點
    D=L(:,2)-L(:,1); 
    logInd=D(2:end).*D(1:end-1)<=0;
    k1=[logInd;false];k2=[false;logInd];
    
    % 對找到的兩個點進(jìn)行插值
    y=(L(k2,1).*L(k1,2)-L(k1,1).*L(k2,2))./(D(k2,:)-D(k1,:));
    Z=[Z,r+abs(y').*1i];
end
plot(Z,'.','markersize',1)
title('Lorenz映射分岔圖')
xlabel('r'),ylabel('|y| where x=y')

4.封面圖繪制

fig=gcf;

% 左圖 
ax1=axes('Parent',fig);
ax1.Position=[1/12 1/12 1/2-1/6 1-1/6];
hold on;grid on
[~,L]=ode45(@(t,L)Lorenz(t,L),0:.01:100,[1;1;1;10;28;8/3]); 
plot3(L(:,1),L(:,2),L(:,3),'Color',[0 0.2510 0.4510 0.5],'LineWidth',1.2)
ax1.XColor=[1,1,1].*.6;
ax1.YColor=[1,1,1].*.6;
ax1.ZColor=[1,1,1].*.6;
ax1.LineWidth=1.5;
ax1.GridAlpha=.09;
ax1.GridLineStyle='-.';
ax1.FontName='cambria';
ax1.DataAspectRatio=[1,1,1];
view([-159,18]);

% 右圖
ax2=axes('Parent',fig);
ax2.Position=[1/2 1/12 1/2-1/18 1-1/6];
hold on;grid on
Z=[]; 
for r=1:500
    % 舍棄前面迭帶的結(jié)果，用后面的結(jié)果畫圖 
    [~,L]=ode45(@(t,L)Lorenz(t,L),[0,1],[1;1;1;10;r;8/3]);  
    [T,L]=ode45(@(t,L)Lorenz(t,L),[0,50],L(end,:)); 
    
    % 找到穿過直線y=x的前后兩個點
    D=L(:,2)-L(:,1); 
    logInd=D(2:end).*D(1:end-1)<=0;
    k1=[logInd;false];k2=[false;logInd];
    
    % 對找到的兩個點進(jìn)行插值
    y=L(k2,2)+(L(k1,2)-L(k2,2)).*D(k2,:)./(D(k2,:)-D(k1,:));
    Z=[Z,r+abs(y').*1i];
end
plot(Z,'.','markersize',1,'Color',[0 0.2510 0.4510 0.5])
ax2.YLabel.String='|y| where x=y';
ax2.YLabel.FontSize=14;
ax2.XColor=[1,1,1].*.4;
ax2.YColor=[1,1,1].*.4;
ax2.ZColor=[1,1,1].*.4;
ax2.LineWidth=1.5;
ax2.GridAlpha=.09;
ax2.GridLineStyle='-.';
ax2.FontName='cambria';

% Lorenz函數(shù)
function dL=Lorenz(t,L)
% L=[x;y;z;a;r;b];
% dL=[dx/dt;dy/dt;dz/dt;0,0,0];
% dz/dt=-a*(x-y)
% dy/dt=x*(r-z)-y
% dz/dt=x*y-b*z
dL=zeros([6,1]);
dL(1)=-L(4)*(L(1)-L(2));
dL(2)=L(1)*(L(5)-L(3))-L(2);
dL(3)=L(1)*L(2)-L(6)*L(3);
dL(4:6)=0;
end