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

第基于Matlab繪制洛倫茲吸引子相圖目錄1.公式及Lorenz函數(shù)2.混沌吸引子圖像3.混沌吸引子圖像4.封面圖繪制洛倫茲吸引子(Lorenzattractor)是由MIT大學(xué)的氣象學(xué)家EdwardLorenz在1963年給出的，他給出第一個(gè)混沌現(xiàn)象蝴蝶效應(yīng)。。。。。。。。廢話不多說(shuō)。

反正咱就是，好看且有用咱就寫代碼，第零部分給出公式。第一部分給出混沌吸引子圖像，第二部分給出龐加萊截面法分岔圖繪制。

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

Lorenz微分方程組定義如下：

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

functiondL=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))

gridon

修飾動(dòng)態(tài)圖：

%ode45求解

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

%修飾及屬性設(shè)置

ax=gca;

holdon;gridon

plhdl=plot3(0,0,0,'Color',[0.98430.85880.53330.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.GridLine;

ax.FontName='cambria';

ax.Color=[000];

ax.DataAspectRatio=[1,1,1];

view([-159,18]);

%循環(huán)繪圖

fori=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=[];

forr=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);

fork2=2:size(L,1)

k1=k2-1;

ifD(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

plot(Z,'.','markersize',1)

title('Lorenz映射分岔圖')

xlabel('r'),ylabel('|y|wherex=y')

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

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

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

其次代碼中用了D(k1)*D(k2)=0來(lái)判斷是否采點(diǎn)，

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

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

另外說(shuō)明一下：

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

是啥。

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

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

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

L(:,1)=L(:,2)-L(:,1);

fork2=2:size(L,1)

k1=k2-1;

ifL(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=[];

forr=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,:));

%找到穿過(guò)直線y=x的前后兩個(gè)點(diǎn)

D=L(:,2)-L(:,1);

logInd=D(2:end).*D(1:end-1)

k1=[logInd;false];k2=[false;logInd];

%對(duì)找到的兩個(gè)點(diǎn)進(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];

plot(Z,'.','markersize',1)

title('Lorenz映射分岔圖')

xlabel('r'),ylabel('|y|wherex=y')

4.封面圖繪制

fig=gcf;

%左圖

ax1=axes('Parent',fig);

ax1.Position=[1/121/121/2-1/61-1/6];

holdon;gridon

[~,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',[00.25100.45100.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.GridLine;

ax1.FontName='cambria';

ax1.DataAspectRatio=[1,1,1];

view([-159,18]);

ax2=axes('Parent',fig);

ax2.Position=[1/21/121/2-1/181-1/6];

holdon;gridon

Z=[];

forr=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,:));

%找到穿過(guò)直線y=x的前后兩個(gè)點(diǎn)

D=L(:,2)-L(:,1);

logInd=D(2:end).*D(1:end-1)

k1=[logInd;false];k2=[false;logInd];

%對(duì)找到的兩個(gè)點(diǎn)進(jìn)行插值

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

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

plot(Z,'.','markersize',1,'Color',[00.25100.45100.5])

ax2.YLabel.String='|y|wherex=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.GridLine;

ax2.FontName='cambria';

%Lorenz函數(shù)

functiondL=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)