Question

用matlab求解二阶微分方程数值解，程序出现错误，求大神指点

先附上自己要求解的实际问题

这是一个超生波凝聚液体中颗粒的问题其中：Pp=940;%为了好看一点，用P代替/rouP1=1000;u1=0.001;dp=40*10^-6;Vp=pi*dp^3/6;k=2;%波数为 2E=1*10^6;%能量密度为 1KW/Lg=9.8;c1=1500;cp=2005.5;下面是我的function filefunction xxxx=ttt(t,x)Pp=940;P1=1000;u1=0.001;dp=40*10^-6;Vp=3.1416*dp^3/6;k=2;%波数为 2E=1*10^6;%能量密度为 1KW/Lg=9.8;c1=1500;cp=2005.5;G=1/3*((5*Pp-2*P1)/(P1-2*Pp)-P1*c1^2/Pp*cp^2);M=P1*Vp/2;fb=g*(P1-Pp)*Vp;xxxx=[x(2);(-3*Vp*k*E*G.*sin(2*k.*x(1))+fb-3*pi*u1*dp.*x(2))/(Pp*Vp+M)];%ode方程这是主程序，（名字是tring2.m）clear all;c1=1500;r=c1/(24.1*10^3);sol=ode15s(@ttt,[0 400],[3/5*r/4;0])；%x(0)=3/5*r/4,x'(0)=0运行之后提示busy;中断程序，matlab提示Operation terminated by user during ode45 (line 264)In trying2 (line 5)sol=ode45(@ttt,[0 400],[3/5*r/4;0]);起初我以为是求解区间太大的问题，但是改小之后还是不行在网上查找了好久，一直不知道怎么解决。不知道是哪里出现问题了

Answer 1

1、把

G=1/3*((5*Pp-2*P1)/(P1-2*Pp)-P1*c1^2/Pp*cp^2);

改成

G=1/3*((5*Pp-2*P1)/(P1+2*Pp)-P1*c1^2/(Pp*cp^2));

表达式的错误导致出现高频振荡，积分步长非常小（大约在10^-8量级），所以很长时间算不出。

 

2、从结果看，仿真时间取50秒就足够了（已进入稳态）。

 

3、在出现busy的情况下，可以考虑让 ode* 函数不返回参数，此种情况下，会自动绘制响应曲线，而且可以随时通过GUI停止运行，便于观察求解过程的具体情况。

Answer 2

G 代码里面第一个分式的分母错了。第二个分式也错了，少括号。