ode45（）函数请教

lvyongjian

大家好！最近我在做非线性振动，需要用matlab进行数值计算，积分等一些列问题。遇到过很多问题，用matlab求解非线性方程，当参数取正的时候，运行没问题，但改其中一个参数的话，程序就报错。不知为什么，我在网上也查阅了一些帖子，好像是矩阵刚度的原因，矩阵接近奇异，试着改了积分函数ode45改为ode23s，ode15s结果还是出现上面那个问题，一直没弄明白。看了在MATLAB论坛及SimWe仿真论坛发的帖子，特向大家请教。
为了表达清楚，具体问题在附件里。谢谢！

lvyongjian

怎么没人呢？

lin2009

刚性问题（简单的说即方程中各物理量的变化幅值相差很大，可达几个数量级），除了应选用合适的ode函数（如ode15,ode23等），初始值是否选择合理？或需限制步长。可以通过odeset(options, 'MaxStep')来实现。
如果还不行的话，那就应该考虑解是否存在的问题。
或者是参数不合适导致微分方程无解。

也可以通过Simulink利用微分／积分模块进行建模仿真，求解该方程。

lvyongjian

3# lin2009
好的，谢谢lin2009。现在我只是做了论文上的一个例子，具体的参数还不是很确定，我再问问老师参数的范围。

x-lv

3# lin2009
弱弱的问一句，用matlab中的Simulink进行建模仿真，难不难？

lin2009

5# x-lv

入门不难。

lvyongjian

6# lin2009
好的，谢谢！

evantong

楼主 我最近也是在做非线性振动。 用simulink 跟 ode45 算出来的结果竟然 不一样。  也不知道 如何解决。   如果 楼主有兴趣 一起 讨论讨论 可以站内 pm 我啊。

lvyongjian

8# evantong
好的，你有时间加我qq：372100615，共同探讨一下，注明信息！

xyz999

在计算楼主的问题之后才发现自己已经不爱动脑子了，想让软件解决一切问题。
  看给出的杜芬方程：
ddx = -c*dx - k1*x - k2*x^2 - k3*x^3 + F*cos(w*t)，其中ddx表示x的二阶导数，dx表示x的一阶导数。
    因为楼主设置的k2=0，所以方程变为
ddx = -c*dx - k1*x - k3*x^3 + F*cos(w*t)
    k3为负值的情况：当x>1时，x^3 > x。设s = -k3*x^3 - k1*x，那么，不管k3取多大，只要为负值，随着 x 的增大，s 将大于零，并且s的增大速度（即ds/dt）是非线性增大的。造成的结果就是 -c*dx  + s > 0。
    如果设原方程的右端项为r，即 r(x) = -c*dx - k1*x - k3*x^3 + F*cos(w*t) ，则原方程可以简写为
ddx = r(x) ，r(x) 是x的增函数。此时可以明白了：系统是一个不稳定系统，x将随着时间变大而越来越大。这也就是matlab给出错误的原因。不是刚性问题，更不是matlab函数的问题，而是系统本身是不稳定造成的。
clear;
t_final=2.58;
c=0.5;
k1=1;
k2=0;
k3=-0.05;
F=10;
q=0;
w=1;
y0=[0;0];
[t1,x1]=ode45('dafen',[0,t_final],y0,[],c,k1,k2,k3,w,F,q);
figure;
plot(t1,x1,'.-');
legend('x','dx')
grid on;
xlabel('Time');

function x=dafen(t,y,flag,c,k1,k2,k3,w,F,q)
x=[y(2)
F*cos(w.*t+q)-c*y(2)-k1*y(1)-k2*y(1)^2-k3*y(1)^3];

lvyongjian

10# xyz999
非常感谢你的分析，我明白了！:)

s552131269

:'(  也遇到了相识的问题

牧兰

matlab数值计算相关书上讲的很详细啊，我也是最近接触这方面的，看完书再看help很有效果