|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
我想求解附件中的一个微分方程,其中L,P,Q,S由参数a ,b ,H, h ,E ,E0 ,mu, et, dens决定,可认为是已知常量,我原本想用龙格库塔法进行求解,编写了如下程序,但是这个程序运行了12个小时还没算完,于是想请教一下高手,求解这个方程还有没有更有效的方法,或者说我编写的程序有没有不当之处。
clear all
global a b H h E E0 mu et dens
%mu表示泊松比,et表示,dens表示密度
tic
a=0.2; b=0.1;H=0.05; h=0.05;E=2*10^10;E0=0.8*10^10;mu=0.4;et=0.35;dens=1500;
x0=[0.1 0];%初值
t_f=[0,10];%积分区间
options=odeset('RelTol',1e-5,'AbsTol',1e-5,'NormControl','on');
[t,X] = ode45(@dfun,t_f,x0,options);
plot(t,X(:,1),'k.','markersize',1)
grid on
xlabel('t'); ylabel('q(t)');
title('图(1) q~t 图')
toc
function dx=dfun(t,x)
global a b H h E E0 mu et dens
L=1/3*et*H^3/4*((pi/a)^4+(pi/b)^4+2*mu*(pi/a)^2*(pi/b)^2)/dens/h;
P=5*((pi/a)^4+(pi/b)^4)*et*H/9/dens/h;
Q=(3*(pi/a)^4+3*(pi/b)^4+2*(mu-1)*(pi/a)^2*(pi/b)^2)*(E*h+E0*H)/9/dens/h;
S=((pi/a)^4+(pi/b)^4+2*mu*(pi/a)^2*(pi/b)^2)*(-E/3*(h^3+3*H^2*h+3*H*h^2)/4-E0/3)/dens/h...
+(-E*(h^3+3*H^2*h+3*H*h^2)/6/(1+mu)-H^3*E0/6/(1+mu))*(pi/a)^2*(pi/b)^2/dens/h;
dx=[x(2); -L*x(2)-P*x(1)^2*x(2)-Q*x(1)^3+S*x(1)]; |
-
|