脚本之家 服务器常用软件
快捷导航
您的位置：首页 → 软件编程 → C 语言 → Matlab洛伦兹吸引子相图

基于Matlab绘制洛伦兹吸引子相图

 更新时间：2022年04月12日 10:41:58   作者：slandarer  
洛伦兹吸引子(Lorenz attractor)是由MIT大学的气象学家Edward Lorenz在1963年给出的。本文将利用Matlab实现洛伦兹吸引子相图的绘制，感兴趣的可以了解一下

洛伦兹吸引子(Lorenz attractor)是由MIT大学的气象学家Edward Lorenz在1963年给出的，他给出第一个混沌现象——蝴蝶效应。。。。。。。。废话不多说。

反正咱就是，好看且有用咱就写代码，第零部分给出公式。第一部分给出 混沌吸引子 图像，第二部分给出庞加莱截面法 分岔图 绘制。

1.公式及Lorenz函数

Lorenz微分方程组定义如下：

非常容易能写出该微分方程组函数：

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.混沌吸引子图像

基础绘图：

[~,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

修饰动态图：

% ode45求解
[~,L]=ode45(@(t,L)Lorenz(t,L),0:.01:100,[1;1;1;10;28;8/3]); 
% 修饰及属性设置
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]);
% 循环绘图
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
    % 舍弃前面迭带的结果，用后面的结果画图 
    [~,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')

代码有一些地方详细讲解一下，首先说明为什么要用

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

的格式进行存储，这样存储可以少构造一个数组，一般情况下我们需要分别存储γ和|y|到两个矩阵，存储为复数形式就可以复平面绘图减少初始化矩阵数量。

其次代码中用了D(k1)*D(k2)<=0来判断是否采点，

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

当D(k1)*D(k2)<=0时说明(x1,y1),(x2,y2)两点分别在 y=x直线两侧。

另外说明一下：

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

是啥。

其实就是构造的两点连线与直线y=x的交点：

PS:为了进一步减少空间复杂度，我们可以将上述函数更改为完全由x,y差值以及y代替，这样就可以直接将中间变量D存储到原来x的位置，减少中间变量的数量：

因此代码可以改写为（当然为了可读性最后并没有采取这个策略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

最后，这部分代码依赖循环我们完全可以将其向量化，即修改为：

Z=[]; 
for r=1:500
    % 舍弃前面迭带的结果，用后面的结果画图 
    [~,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];
    
    % 对找到的两个点进行插值
    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
    % 舍弃前面迭带的结果，用后面的结果画图 
    [~,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];
    
    % 对找到的两个点进行插值
    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函数
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绘制洛伦兹吸引子相图的详细内容，更多关于Matlab洛伦兹吸引子相图的资料请关注脚本之家其它相关文章！

您可能感兴趣的文章:

相关文章

  • C++指向函数的指针实例解析

    C++指向函数的指针实例解析

    这篇文章主要介绍了C++指向函数的指针,需要的朋友可以参考下
    2014-07-07
  • Linux下semop等待信号时出现Interrupted System Call错误(EINTR)解决方法

    Linux下semop等待信号时出现Interrupted System Call错误(EINTR)解决方法

    本篇文章是对在Linux下semop等待信号时出现Interrupted System Call错误(EINTR)的解决方法进行了详细的分析介绍，需要的朋友参考下
    2013-05-05
  • C语言单双线性及循环链表与实例

    C语言单双线性及循环链表与实例

    这篇文章主要介绍了C语言的单链表、双链表和循环链表，还有一些相关的实例，感兴趣的同学可以借鉴一下
    2023-03-03
  • C语言算法练习之佩奇借书

    C语言算法练习之佩奇借书

    这篇文章主要该大家分享C语言佩奇借书的练习，文章主要通过描述佩奇借书的问题然后确定程序框架将结果运算出来，下面来看详细内容吧,需要的朋友可以参考一下
    2022-03-03
  • C语言超详细解析函数栈帧

    C语言超详细解析函数栈帧

    在C语言中，每一个正在运行的函数都有一个栈帧与其对应，栈帧中存储的是该函数的返回地址和局部变量。从逻辑上讲，栈帧就是一个函数执行的环境：函数参数、函数的局部变量、函数执行完后返回到哪里等等
    2022-03-03
  • 关于C++中0是十进制还是八进制的问题

    关于C++中0是十进制还是八进制的问题

    本篇文章中，小编将为大家介绍关于C++中0是十进制还是八进制的问题，有需要的朋友可以参考一下
    2013-04-04
  • C++实现LeetCode(27.移除元素)

    C++实现LeetCode(27.移除元素)

    这篇文章主要介绍了C++实现LeetCode(27.移除元素),本篇文章通过简要的案例,讲解了该项技术的了解与使用,以下就是详细内容,需要的朋友可以参考下
    2021-07-07
  • C语言的变量类型及内存大小详解

    C语言的变量类型及内存大小详解

    这篇文章主要介绍了CC和C++变量类型及内存大小,是C++入门学习中的基础知识,需要的朋友可以参考下，希望能够给你带来帮助
    2021-09-09
  • C++ LeeCode题目:比特位计数和买卖股票的最佳时机

    C++ LeeCode题目:比特位计数和买卖股票的最佳时机

    这篇文章主要介绍了基于C语言计算比特位计数和买卖股票的最佳时机,文中通过示例代码介绍的非常详细，对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下
    2021-07-07
  • C语言设计三子棋小游戏

    C语言设计三子棋小游戏

    这篇文章主要为大家详细介绍了C语言设计三子棋小游戏，文中示例代码介绍的非常详细，具有一定的参考价值，感兴趣的小伙伴们可以参考一下
    2021-04-04

最新评论