ljf8848 发表于 2007-6-28 10:53

关于 ode45 的编程问题

function dy=MASATU(t,y)
mu=0.5;
e=0.0002;
w=1000;
m=5;
c1=0.1;
c2=50;
k1=15775;
k2=6300*15775;
tao=0.001;
W=sqrt(y(1)^2+y(3)^2);
if W>=tao
    Fcon1=k2*(W-tao)*(y(1)-mu*y(3))/W+c2*y(2);
    Fcon2=k2*(W-tao)*(y(1)+mu*y(3))/W+c2*y(4);
else
    Fcon1=0;
    Fcon2=0;
end
dy=[y(2);
    e*w^2*cos(w*t)-Fcon1/m-c1*y(2)/m-k1*y(1)/m;
    y(4);
    e*w^2*sin(w*t)-Fcon2/m-c1*y(4)/m-k1*y(3)/m];
执行语句
=ode45(@MASATU,,);plot(T,Y);

[ 本帖最后由 eight 于 2007-6-28 11:44 编辑 ]

xjzuo 发表于 2007-6-28 11:07

纵坐标的值太大,显然是不合理的,将原问题用word上传一下看看．

－－－似乎很多人都遇到了这种问题, 看来有时间的话，要好好考虑一下这个问题了.

ljf8848 发表于 2007-6-28 11:16

谢谢大侠，
公式上传。
可能效果不好。先请原谅。

ljf8848 发表于 2007-6-28 13:07

每人给点建议阿！着急啊。还有一个月就毕业了。结果出不赖！

yanyongju 发表于 2007-6-28 13:41

这个问题，你可以用RK来计算，看看单步的值，就知道问题出在哪里了！
ode45的用法没有问题，剩下的就是方程的写法了，自己仔细看看

xjzuo 发表于 2007-6-28 20:44

用另一个自编的刚性求解器算了一下,不过效果仍不太好.
所以建议自己再仔细检查一下参数是否合理,模型是否正确.
另：改变区间或步长求解一下，你就会发现问题所在．

[ 本帖最后由 xjzuo 于 2007-6-28 20:45 编辑 ]

咕噜噜 发表于 2007-6-28 20:48

回复 #2 xjzuo 的帖子

我也遇到这个问题了，一直高不明白什么原因

flybaly 发表于 2007-6-29 13:19

平衡点不稳定吧。。。

appleseed05 发表于 2007-6-29 13:35

RK对初参数很敏感的，建议改成其他动力学数值算法，比如Newmark等试试。
如果非要用ode的话，改成ode23t试试
页: [1]
查看完整版本: 关于 ode45 的编程问题