想请教关于ode函数的一个问题

hfut_dotaman · 发表于 2010-11-23 14:11

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

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

x

本帖最后由 hfut_dotaman 于 2010-11-23 14:28 编辑

我用ode函数解方程组时间设为【0,20】，但方程解到0.0048s就解不下去了，提示如下：
Warning: Failure at t=4.775313e-003. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.696532e-017) at time t.

原方程为：function dx=wei(t,x);
global kv sf v2 u a b m iz k1 k2 pz1 pz2 v1 u2
o1=x(1)+a*x(2)/u-sf;
o2=x(1)-b*x(2)/u;
oi1=k1*tan(o1)/(pz1*u2);
oi2=k2*tan(o2)/(pz2*u2);
fy1=pz1*u2*(oi1-oi1*oi1/3+oi1^3/27);
fy2=pz2*u2*(oi2-oi2*oi2/3+oi2^3/27);
dx(1)=(((fy1+kv*v1)*cos(sf)+(fy2+kv*v2))/m-u*x(2))*u;
dx(2)=(a*(fy1+kv*v1)*cos(sf)-b*(fy2+kv*v2))/iz;
dx=[dx(1);dx(2)];

ode程序为：hold on;global kv sf v2 u a b m iz k1 k2 pz1 pz2 v1 u2
kv=5000;sf=0.0872;v2=0.0174;u=100*1000/3600;a=1.84;b=1.88;m=3018;iz=10437;k1=23000;k2=38000;pz1=1500*9.8;pz2=1518*9.8;v1=0.003;u2=0.8;
N=200;
Tn=[];
options=odeset('RelTol',1e-4);
x0=[0,0];
[t,x]= ode15s('wei',[0,10],x0,options);

ode45 ode15我都试过了都不行时间步长改为0:0.0001:20也不行，RelTol改小也不行，原方程组没错，希望大家多指点，急啊。

hfut_dotaman · 发表于 2010-11-23 14:13

是要做分岔图的，结果oed函数都没搞对，下一步没法做了。
分岔图的纵坐标最后都是10的13次方这个数量级的，请大侠指点。

octopussheng · 发表于 2010-11-23 14:14

本帖最后由 octopussheng 于 2010-11-23 14:14 编辑

这是微分方程的问题，计算出现了刚性，要调整里面的参数了，以前也遇到这种情况的。

hfut_dotaman · 发表于 2010-11-23 14:16

这么快就回复了厉害！

hfut_dotaman · 发表于 2010-11-23 14:17

你是说kv=5000;sf=0.0872;v2=0.0174;u=100*1000/3600;a=1.84;b=1.88;m=3018;iz=10437;k1=23000;k2=38000;pz1=1500*9.8;pz2=1518*9.8;v1=0.003;u2=0.8;
这些参数需要改一些？

octopussheng · 发表于 2010-11-23 20:15

可试着调调吧，这种问题是比较麻烦的。

hfut_dotaman · 发表于 2010-11-23 20:20

多谢我不喜欢字数限制 8个字

gghhjj · 发表于 2010-11-25 09:37

可以尝试ode15s 和 ode23s

账号		自动登录	找回密码
密码			我要加入

[稳定性与分岔] 想请教关于ode函数的一个问题

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