lvyongjian 发表于 2010-10-8 22:13:02

ode45()函数请教

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

lvyongjian 发表于 2010-10-9 09:22:51

怎么没人呢?

lin2009 发表于 2010-10-9 10:04:00

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

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

lvyongjian 发表于 2010-10-9 18:27:16

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

x-lv 发表于 2010-10-9 18:33:16

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

lin2009 发表于 2010-10-9 19:53:37

5# x-lv

入门不难。

lvyongjian 发表于 2010-10-11 20:32:49

6# lin2009
好的,谢谢!

evantong 发表于 2010-10-13 02:32:27

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

lvyongjian 发表于 2010-10-13 16:14:36

8# evantong
好的,你有时间加我qq:372100615,共同探讨一下,注明信息!

xyz999 发表于 2010-10-14 15:03:48

在计算楼主的问题之后才发现自己已经不爱动脑子了,想让软件解决一切问题。
看给出的杜芬方程:
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=;
=ode45('dafen',,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 发表于 2010-10-23 15:27:49

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

s552131269 发表于 2012-4-14 22:00:01

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

牧兰 发表于 2012-4-16 09:56:30

matlab数值计算相关书上讲的很详细啊,我也是最近接触这方面的,看完书再看help很有效果
页: [1]
查看完整版本: ode45()函数请教