关于 ode45 的编程问题

ljf8848

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];
执行语句
[T,Y]=ode45(@MASATU,[0:0.001:1],[1;1;1;1]);plot(T,Y);

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

xjzuo

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

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

ljf8848

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

ljf8848

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

yanyongju

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

xjzuo

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

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

咕噜噜

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

flybaly

平衡点不稳定吧。。。

appleseed05

RK对初参数很敏感的，建议改成其他动力学数值算法，比如Newmark等试试。
如果非要用ode的话，改成ode23t试试