声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

楼主: mechanic05

[稳定性与分岔] 请教关于分岔程序

[复制链接]
发表于 2007-7-17 16:08 | 显示全部楼层
太乱,感觉有点不是很好这个分叉图
回复 支持 反对
分享到:

使用道具 举报

发表于 2007-7-17 16:25 | 显示全部楼层
你这个图和我在95楼的差不多吗!

我上面的那两个图,一个是a,一个是-a,效果还不错,呵呵!
程序和你的不大一样!
function xdot=dddd(t,x)
eta=0.1 ; a= -48.704; b= 24.35e4;gamma= 1;
global omega;
xdot=[x(2);-a*x(1)-b*x(1).^3-eta*x(2)+gamma*cos(omega*t)];


clear all
global omega;
t0=[0 200];%积分时间
y0=[0,0];

%bifurcation
for omega=0.01:0.01:4; %omega的变化精度
    [t,y]=ode45('dddd',t0,y0);
    [Xmax]=getmax(y(:,2));

    plot(omega,Xmax,'b','markersize',1)
    hold on
    clear Xmax
end


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
 楼主| 发表于 2008-1-7 10:17 | 显示全部楼层

回复 #107 octopussheng 的帖子

新年好!:handshake
不好意思,一直没回复是因为忙别的工作了.现在闲下来整理时又拿出这个工作,决定再次研究这个分岔问题.发现其实无水前辈给的分岔程序很适合我的这个工作.我也研读了彭老师的书.基本弄懂了分岔及分岔程序的意义.下面是我贴出来的分岔图,看起来还不错.是关于参数a的.抽空准备再用前辈给的getmax方法做一下,比较一下两个程序的差别.
rezaihe.jpg
发表于 2008-1-7 13:27 | 显示全部楼层
这个图还可以啊,我建议你在5~26这个区间做个放大的图看看!
 楼主| 发表于 2008-1-7 21:05 | 显示全部楼层

回复 #109 octopussheng 的帖子

:handshake 多谢您的建议!!
我正在算呢.本来还想贴一个关于gamma的分岔图,结果由于图太大(大约247K)就没有贴.图有办法变小吗?
发表于 2008-1-7 21:24 | 显示全部楼层
用acdsee之类的工具,可以将图像的尺寸设置的小一些,这样就可以上传了!
 楼主| 发表于 2008-1-10 18:05 | 显示全部楼层
各位新年好!!
1)我把上面分岔图特殊部分放大了.现在传上来,麻烦前辈给看看.从图上看,当参数取8.05以前时应该是1周期运动;当参数取8.13时应该是多周期运动;当参数取8.5时系统应该是混沌运动,对吧?
2)我把参数取4时的各图形及程序贴上来,麻烦前辈高人给看看,程序有无不妥之处.为什么当我取参数为8.13时(程序其他部分不动)并不是想象的多周期运动,而是和参数取4时结果相似?问题是否是画图点数取的范围有技巧?麻烦指点.非常感谢!!

clc;clear;
nt=4;
eta=0.1; b=24.35*10^4; omega=2; gamma=1;T=2*pi/omega;
str{1}='庞加莱截面-周期1吸引子';
    [t,x]=ode23('nnn',[0:T/1000:200*T],[0.02,0.02],[],nt,eta,b,omega,gamma);
    figure
    subplot(2,1,1)
    plot(t,x(:,1));
    title('位移曲线')
    xlabel('t');ylabel('x');
   
    subplot(2,2,3)
    plot(x(50000:end,1),x(50000:end,2));
    axis([-0.02 0.02 -0.2 0.2])
    xlabel('x');ylabel('dx/dt');
    title('相图')
   
   subplot(2,2,4)
   axis([-0.06 0.06 -0.05 0.05])
   hold on
   for i=50000:1000:100000
       plot(x(i,1),x(i,2),'k.');
   end
   title(str{1});
end
rezaihe5.jpg
rezaihe_zhouqi1.jpg
 楼主| 发表于 2008-1-10 18:14 | 显示全部楼层

回复 109楼 的帖子

麻烦前辈给看看!!谢谢!!
发表于 2008-1-10 19:44 | 显示全部楼层
1)我把上面分岔图特殊部分放大了.现在传上来,麻烦前辈给看看.从图上看,当参数取8.05以前时应该是1周期运动;当参数取8.13时应该是多周期运动;当参数取8.5时系统应该是混沌运动,对吧?
2)我把参数取4时的各图形及程序贴上来,麻烦前辈高人给看看,程序有无不妥之处.为什么当我取参数为8.13时(程序其他部分不动)并不是想象的多周期运动,而是和参数取4时结果相似?问题是否是画图点数取的范围有技巧?麻烦指点.非常感谢!!
第一个问题你说的对,但是第二个就不知道为什么了,系统的状态不能按照想象来,不过你的程序中为什么取中间的那部分来画图,for i=50000:1000:100000,,[0:T/1000:200*T]不是200000个点吗?不过从相图看倒是比较有意思
发表于 2008-1-11 07:50 | 显示全部楼层
相图看起来是一个多周期的运动哦,怎么Poincare截面只有一个点啊?不是很理解!


建议新开一个帖子吧,这个帖子已经够长了!
 楼主| 发表于 2008-1-11 09:14 | 显示全部楼层
谢谢各位前辈关注!!非常感谢!!:handshake
那我就另外开一个帖子.其实咕噜噜前辈的疑问正是我的疑问.为什么要从中间取值呢?这个程序的写法,我是仿照前辈给的那本彭老师的书编写的.我也对这个有疑问.麻烦大家看看.互相交流、提高!!:@) 谢谢了!!
我这里就我一个人做这个,所以有问题也没人可以讨论。导师也没做。所以我非常感谢群里的前辈及同仁。
祝贺大家新年快乐!!
发表于 2008-6-9 10:02 | 显示全部楼层

回复115楼

同意你的观点,这是个3周期运动吧
发表于 2008-6-25 16:16 | 显示全部楼层

请教一下我的分岔程序错在哪?

function dy=dufen(t, x, E)
miu=0.168;
wi2=0.5;
bita=0.5;
p=1;
dy=[x(2);E*sin(p*t)-miu*x(2)+wi2*x(1)-bita*x(1)^3]


close all, clear all, clc
T=2*pi;
f=0.1:0.001:0.2;
for i=1:length(f)
[t,x]=ode45('dufen',[0:T/100:50*T],[0 0],[],f(i));
plot(f(i),x(1:100:end,1),'--');hold on
end
 我是参照上面的程序写的
怎么就运行不了啊?
谢谢前辈指点
 楼主| 发表于 2009-3-22 01:31 | 显示全部楼层

回复 118楼 jieshigood 的帖子

你的E没有定义啊,好像。
发表于 2010-8-31 08:37 | 显示全部楼层
回复 octopussheng 的帖子
请问这个程序中的r是什么意思呢,一般研究的不是取系统中的某个参数作为分岔参数,进而看它取何值时产生hopf分岔吗,就像前面无水兄发的那个程序似的,这个是否也能改变成那种形式呢,谢谢,我现在正在困惑与此,请指点

   
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-5-19 09:07 , Processed in 0.057430 second(s), 20 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表