声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2003|回复: 10

[稳定性与分岔] 大家帮忙给看看这分岔程序

[复制链接]
发表于 2010-8-13 09:52 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
最近做分岔,解微分方程,总是出问题,各位帮忙给看看,指点一下
原系统是这样的
file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/2AYXH8$IC8PJ3~{UU{RJ(EP.jpg

程序如下:
function v=yuanxitong(t,x)
%参数取值
g1=9;g2=2;d1=4;d2=3;beta1=4;beta2=7;h1=5;h2=-6;
f2=4;c1=7;c2=5;b1=8;b2=9;l1=9;l2=-7;a2=10;
w2=1;u1=6;u2=5;
a1=-22.9126;   %bifurcation point is a1
w1=0;f1=0;
omega1=sqrt(-44.24268);
omega2=sqrt(-3.33333);
omega=2;

v=zeros(4,1);
v(1)=x(2);
v(2)=-omega1^2*x(1)-a1*x(1)*cos(omega*t)-b1*x(1)*sin(omega*t)+g2*x(3)+h2*x(3)*cos(omega*t)...
    +l2*x(3)*sin(omega*t)+u2*x(4)+beta2*x(4)^3-c1*x(2)-d1*x(2)^3+f1*cos(w1*t);
v(3)=x(4);
v(4)=-omega2^2*x(3)-a2*x(3)*cos(omega*t)-b2*x(3)*sin(omega*t)+g1*x(1)+h1*x(1)*cos(omega*t)...
    +l1*x(1)*sin(omega*t)+u1*x(2)+beta1*x(2)^3-c2*x(4)-d2*x(4)^3+f2*cos(w2*t);

clc
clear all

[t,x]=ode45(@yuanxitong,[0:0.001:200],[0.1;0.1;0.1;0.1]);

figure(1)
plot(x(1,:),x(2,:));
figure(2)
plot(t,x(1,:));   
由于刚做这东西,好多还不明白,高手给运行看看,谢谢

[ 本帖最后由 kangarooli 于 2010-8-13 16:35 编辑 ]
回复
分享到:

使用道具 举报

发表于 2010-8-13 10:21 | 显示全部楼层

回复 楼主 kangarooli 的帖子

sol=ode45(@yuanxitong,[1,2.5],[0.1;0.1;0.1;0.1],[0:0.001:20]);   语法是不是有问题啊,那个[1,2.5]表示什么
另外你这个是分岔程序吗?感觉你是想解相图和时间序列
 楼主| 发表于 2010-8-13 10:42 | 显示全部楼层

回复 沙发 hsfy919 的帖子

这应该是一系列时间取值吧,我也不是很明白,这个值是我随便取得,你说得对,我是想解相图
 楼主| 发表于 2010-8-13 15:51 | 显示全部楼层

回复 楼主 kangarooli 的帖子

程序运行老是提示如下
Warning: Failure at t=1.132377e-001.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.220446e-016) at time t.
> In ode45 at 355
  In fenchatu at 3

[ 本帖最后由 kangarooli 于 2010-8-13 16:17 编辑 ]
发表于 2010-8-13 16:23 | 显示全部楼层
建议你把后面那段改成下面的程序试试
clc
clear all

[t,x]=ode45(@yuanxitong,[0,200],[0.1;0.1;0.1;0.1],[]);

figure(1)
plot(x(:,1),x(:,2));
figure(2)
plot(t,x(:,2));


另外请你仔细检查一下你的方程,包括参数设置,为什么你定义了b2,sigma1,sigma2,但在方程里都没出现,至于你所得那个警告我估计是和你方程参数有关

评分

1

查看全部评分

 楼主| 发表于 2010-8-13 16:43 | 显示全部楼层

回复 5楼 hsfy919 的帖子

谢谢,是我写错了,不过改过来还是不行
发表于 2010-8-13 16:48 | 显示全部楼层

回复 6楼 kangarooli 的帖子

如果你想画相图,那我给你改过的程序模式应该没问题,你可以调试积分时间和初值,如果还是不行,我个人觉得你就应该从方程入手了,比如方程参数是否合理,模型是否符合实际等等
 楼主| 发表于 2010-8-13 17:15 | 显示全部楼层

回复 7楼 hsfy919 的帖子

虽然还没得到想要的,但现在有东西了,不过我有点不明白,问什么改成x(:,1)的形式呢,这跟x(1,:)有什么意义

[ 本帖最后由 kangarooli 于 2010-8-13 17:17 编辑 ]
发表于 2010-8-13 22:49 | 显示全部楼层

回复 8楼 kangarooli 的帖子

x(:,1)表示矩阵x的第一列,x(1,:)表示矩阵x的第一行
常微分积分出来的x是列向量
发表于 2010-10-18 09:12 | 显示全部楼层
回复 hsfy919 的帖子

是的,有时候时间取到个别值不满足
发表于 2010-10-19 11:22 | 显示全部楼层
本帖最后由 Vickyvictoria 于 2010-10-19 11:23 编辑

Warning: Failure at t=1.378812e+000.  Unable to meet integration tolerances
without reducing the step size below the smallest value allowed (3.552714e-015)
at time t.

从运行的情况来看,你的程序根本就没有求解成功
根据上面的提示来分析,要嘛是你的方程太特殊,要嘛是你的方程有问题
总之是用ode45无法得到收敛解
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-17 17:40 , Processed in 0.063171 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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