基于Matlab繪制洛倫茲吸引子相圖
更新時(shí)間：2022年04月12日 10:41:58 作者：slandarer
洛倫茲吸引子(Lorenz attractor)是由MIT大學(xué)的氣象學(xué)家Edward Lorenz在1963年給出的。本文將利用Matlab實(shí)現(xiàn)洛倫茲吸引子相圖的繪制，感興趣的可以了解一下
洛倫茲吸引子(Lorenz attractor)是由MIT大學(xué)的氣象學(xué)家Edward Lorenz在1963年給出的，他給出第一個(gè)混沌現(xiàn)象——蝴蝶效應(yīng)。。。。。。。。廢話不多說(shuō)。
反正咱就是，好看且有用咱就寫(xiě)代碼，第零部分給出公式。第一部分給出 混沌吸引子 圖像，第二部分給出龐加萊截面法 分岔圖 繪制。
1.公式及Lorenz函數(shù)
Lorenz微分方程組定義如下：
非常容易能寫(xiě)出該微分方程組函數(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
修飾動(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;
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é)果畫(huà)圖
[~,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ì)講解一下，首先說(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ù)量：
因此代碼可以改寫(xiě)為（當(dāng)然為了可讀性最后并沒(méi)有采取這個(gè)策略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é)果畫(huà)圖
[~,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)<=0;
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];
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é)果畫(huà)圖
[~,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)<=0;
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];
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
以上就是基于Matlab繪制洛倫茲吸引子相圖的詳細(xì)內(nèi)容，更多關(guān)于Matlab洛倫茲吸引子相圖的資料請(qǐng)關(guān)注腳本之家其它相關(guān)文章！
相關(guān)文章
c++ 子類構(gòu)造函數(shù)初始化及父類構(gòu)造初始化的使用
這篇文章主要介紹了c++ 子類構(gòu)造函數(shù)初始化及父類構(gòu)造初始化的使用，文中通過(guò)示例代碼介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來(lái)一起學(xué)習(xí)學(xué)習(xí)吧2020-07-07
并查集 是一種樹(shù)型的數(shù)據(jù)結(jié)構(gòu),用于處理一些不相加集合的合并和查詢問(wèn)題。本文給大家分享C++并查集常用操作及算法實(shí)現(xiàn)，感興趣的朋友跟隨小編一起看看吧2021-07-07
Cocos2d-x 3.x入門教程（一）：基礎(chǔ)概念
這篇文章主要介紹了Cocos2d-x 3.x入門教程（一）：基礎(chǔ)概念,本文講解了Director、Scene、Layer、Sprite等內(nèi)容,需要的朋友可以參考下2014-11-11
這篇文章主要介紹了C++編程中的匿名類類型和位域,是C++入門學(xué)習(xí)中的基礎(chǔ)知識(shí),需要的朋友可以參考下2016-01-01
C++中const char*、char const*、char * const三者的區(qū)別
這篇文章主要介紹了C++中const char*、char const*、char * const三者的區(qū)別，文中通過(guò)示例代碼介紹的非常詳細(xì)，具有一定的參考價(jià)值，感興趣的小伙伴們可以參考一下2021-09-09
C++/STL實(shí)現(xiàn)判斷平面內(nèi)兩條線段的位置關(guān)系代碼示例
這篇文章主要介紹了C++/STL實(shí)現(xiàn)判斷平面內(nèi)兩條線段的位置關(guān)系代碼示例，具有一定參考價(jià)值，需要的朋友可以了解下。2017-11-11
以下是對(duì)C/C++中的SIZEOF進(jìn)行了詳細(xì)的分析介紹，需要的朋友參考下2013-07-07
C++實(shí)現(xiàn)LeetCode(122.買股票的最佳時(shí)間之二)
這篇文章主要介紹了C++實(shí)現(xiàn)LeetCode(122.買股票的最佳時(shí)間之二),本篇文章通過(guò)簡(jiǎn)要的案例,講解了該項(xiàng)技術(shù)的了解與使用,以下就是詳細(xì)內(nèi)容,需要的朋友可以參考下2021-07-07
Opencv學(xué)習(xí)教程之漫水填充算法實(shí)例詳解
這篇文章主要給大家介紹了Opencv學(xué)習(xí)教程之漫水填充算法的相關(guān)資料，文中給出了詳細(xì)的示例代碼供大家參考學(xué)習(xí)，對(duì)大家具有一定的參考價(jià)值，需要的朋友們下面跟著小編一起來(lái)學(xué)習(xí)學(xué)習(xí)吧。2017-06-06
今天小編就為大家分享一篇關(guān)于一圖總結(jié)C++中關(guān)于指針的那些事，小編覺(jué)得內(nèi)容挺不錯(cuò)的，現(xiàn)在分享給大家，具有很好的參考價(jià)值，需要的朋友一起跟隨小編來(lái)看看吧2019-02-02