声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

声振论坛 展示 振动理论 非线性振动 查看内容

【总结帖】分岔图绘制不同方法的总结、比较

2015-10-12 04:40| 发布者: aspen| 查看: 1344| 评论: 89|原作者: octopussheng|来自: 声振论坛

摘要: 经过近期的研究发现,目前对于系统单参数分岔图的计算共有以下的几种方法: 1)最大值法 即对系统微分方程(组)进行求解,对求解的结果用getmax函数进行取点,并绘图。 2)Poincare截面法 对系统参数的每一 ...
经过近期的研究发现,目前对于系统单参数分岔图的计算共有以下的几种方法:
1)最大值法
     即对系统微分方程(组)进行求解,对求解的结果用getmax函数进行取点,并绘图。
2)Poincare截面法
    对系统参数的每一次取值,绘制其Poincare截面,进而得到其分岔图。
    这种方法需要注意的是,自治系统的Poincare截面是选取一超平面,平面上点的分布即构成一Poincare截面,非自治系统的Poincare截面则是根据系统激励的频率进行取点并绘图。
本帖将以Lorenz系统为例,对这两种方法进行比较
首先对第二种方法进行阐述。
编程如下(matlab)
Lorenz系统:
function dy = Lorenz(t,y)
% Lorenz系统
% 系统微分方程:
%        dx/dt = -a(x-y)
%        dy/dt = x(r-z)-y
%        dz/dt = xy-bz
%    a=y(4)
%    r=y(5)
%    b=y(6)
dy=zeros(6,1);
dy(1)=-y(4)*(y(1)-y(2));
dy(2)=y(1)*(y(5)-y(3))-y(2);
dy(3)=y(1)*y(2)-y(6)*y(3);
dy(4)=0;
dy(5)=0;
dy(6)=0;
随r的分岔图求解程序:——按照x=y平面取截面
function Lorenz_bifur_r
Z=[];
for r=linspace(1,500,1000);
    % 舍弃前面迭带的结果,用后面的结果画图
    [T,Y]=ode45('Lorenz',[0,1],[1;1;1;16;r;4]);  
    [T,Y]=ode45('Lorenz',[0,50],Y(length(Y),:));
    Y(:,1)=Y(:,2)-Y(:,1);
    % 对计算结果进行判断,如果点满足x=y,则取点
    for k=2:length(Y)
        f=k-1;
        if Y(k,1)<0   
            if Y(f,1)>0
                y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));
                Z=[Z r+abs(y)*i];
            end  
        else     
            if Y(f,1)<0  
                y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1));   
                Z=[Z r+abs(y)*i];   
            end  
        end
    end
end
plot(Z,'.','markersize',1)
title('Lorenz映射分岔图')
xlabel('r'),ylabel('|y| where x=y')
结果如图1所示。
[/url]
getmax法——取最大值法
function [Xmax] = getmax(y)
a=length(y);
j=1;
for i=(a-1)/2:a
    b=(y(i,1)-y(i-2,1))/2;
    c=(y(i,1)+y(i-2,1))/2-y(i-1,1);
    if y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)&c==0
        Xmax(j)=y(i-1,1);
        j=j+1;
    elseif y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)
        Xmax(j)=y(i-1,1)-b^2/(4*c);
        j=j+1;
    end
end
function Lorenz_bifur_r_getmax
% 最大值法求解分岔图
clear all
t0=[0 100];%积分时间
%bifurcation
for r=linspace(1,500,1000);  %r的变化精度
    [t,y]=ode45('Lorenz',t0,[1;1;1;16;r;4]);
    [Xmax]=getmax(y(:,1));
    plot(r,Xmax,'b','markersize',1)
    hold on
    clear Xmax
end
计算结果如图2所示!
[/url]
最后上传一下参考计算机仿真第22卷第12期上一篇文章“李雅普诺夫指数的研究与仿真”中Lorenz系统的分岔图计算结果,大家比较一下即可看出孰优孰劣了!
[/url]
附件:
发表评论

最新评论

引用 octopussheng 2008-3-14 09:59
下面也请无水、小咕把你们有的关于非自治系统的频闪绘制分岔图的例子也贴一贴吧!让大家都能学习学习!呵呵!
引用 16443 2008-3-14 15:00
支持!
引用 咕噜噜 2008-3-14 15:52
回头一定贴出来,不过论坛里面这类程序其实很多了,总结一下就足够大家用了
我这几天先做试验,oct你先帮我看看实验
引用 xiaoqiu810818 2008-3-14 20:45
:handshake 各位师兄辛苦了。请允许我表达对你们真挚的感谢。
引用 yina_111 2008-3-18 09:50
楼主,那个图3的Lorenz方程的Lyapunov指数图的程序能给我吗?我想学习一下,非常感谢!
引用 liliangbiao 2008-3-18 15:56
这些程序我觉得都不算太好!
引用 octopussheng 2008-3-19 08:00
不知道你这个是怎么画的?能否分享一下经验!
引用 xiaoqiu810818 2008-3-20 21:35
等待liliangbiao 大哥的分享
引用 咕噜噜 2008-3-21 09:05
是啊,不要光是把图贴出来让大家眼馋,分享一下程序啊
引用 无水1324 2008-3-21 10:56
大家自己捉摸一下就好了,等别人的程序很难,虽然我也期待,:@L
引用 octopussheng 2008-3-21 13:36
个人看法——他这个图应是步长取的很细算得的!
引用 无水1324 2008-3-21 14:23
oct,给出系统的详细参数,我给你算一个吧,一定比较好的分岔图
引用 21172485 2008-3-21 15:12
看的我手痒痒啊,现在没时间,等过一阵我也画一个,看看效果:lol
引用 无水1324 2008-3-21 20:22
就以这个系统为例,oct给出参数,会做得都做一个出来,我就不相信这个问题不能彻底解决!
引用 octopussheng 2008-3-22 09:50
呵呵,在上面就有的哦,a=16;b=4
现在正在想怎么画到我写的参考文献中的最大Lyapunov指数谱图呢!有结果了马上贴出来!呵呵!
引用 无水1324 2008-3-22 14:21
哎,没有办法了,用你得参数,我始终算得跟你得不同。大哥你搞清楚没有,那篇论文上的参数是什么?a=10,b=8/3。你得参数都不一样,怎么会得到一样的图呢,使我看错了还是?:@L
引用 octopussheng 2008-3-22 14:51
汗,参数搞错了哈,正在算,马上把结果贴出来,看看差别多大!
引用 octopussheng 2008-3-22 16:55
纠正前面的错误,参考资料上面的计算条件是a=10,b=8/3,我前面的程序计算条件不一样,算了一下午,结果如下:
引用 无水1324 2008-3-22 17:46
是不是初值的影响,导致图有差异?

查看全部评论(89)

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-5-8 10:34 , Processed in 0.046963 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

返回顶部