ode45求解微分方程组，一直报错，请各位大侠帮我！
（1）求解问题为：解8个方程构成的微分方程组，微分方程的系数是时变的。（2）我的解题思路是：在M文件建立函数“function dx=myfun8(t,x)”，再利用“=ode45('myfun8',,)”求解，初值均为0。
（3）程序问题：（a）用ode45显示结果为“out of memory”,并且用时明显很长；
(b)若荣ode15s等解刚性方程组的命令，显示结果为“Warning: Failure at t=7.619368e-006.Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.706943e-020) at time t.
> In ode15s at 738”
（c）若改变初值，比如把后两个初值改为：0.1、0.1，可以用ode15s求出结果，但初值改变了，也不知道得到的结果是否准确。
纠结了很久了，实在是初学没有经验，麻烦大侠多指导我，不胜感激。
程序代码为：
function dx=myfun8(t,x)
U1=;
%定义初值
Lp1=1e-3; Lp2=1e-3; Lp3=1e-3;
rp1=2.5; rp2=2.5; rp3=2.5;
rc1=2274; rc2=2274; rc3=2274;
Np1=78; Np2=78; Np3=78;
Lc1=3.017; Lc2=3.017; Lc3=3.017; Lc4=3.160; Lc5=3.160; Lc6=3.160; Lc7=3.160;
Ac1=0.5555; Ac2=0.5555; Ac3=0.5555; Ac4=0.5635; Ac5=0.5635; Ac6=0.5635; Ac7=0.5635;
Nv=;
Nf=;
C=;
ud1=1/(eps+0.61425*cosh(4.55*x(7)/Ac1));ud2=1/(eps+0.61425*cosh(4.55*x(8)/Ac2));ud3=1/(eps+0.61425*cosh(4.55*(-x(7)-x(8))/Ac3));
ud4=1/(eps+0.61425*cosh(4.55*x(7)/Ac4));ud5=1/(eps+0.61425*cosh(4.55*x(7)/Ac5));ud6=1/(eps+0.61425*cosh(4.55*(-x(7)-x(8))/Ac6));ud7=1/(eps+0.61425*cosh(4.55*(-x(7)-x(8))/Ac7));
Rmd1=Lc1/(ud1*Ac1+eps);Rmd2=Lc2/(ud2*Ac2+eps);Rmd3=Lc3/(ud3*Ac3+eps);Rmd4=Lc4/(ud4*Ac4+eps);Rmd5=Lc5/(ud5*Ac5+eps);Rmd6=Lc6/(ud6*Ac6+eps);Rmd7=Lc7/(ud7*Ac7+eps);
a1=Rmd1*Rmd2+Rmd1*Rmd3+Rmd1*Rmd6+Rmd1*Rmd7+Rmd2*Rmd3+Rmd2*Rmd6+Rmd2*Rmd7+Rmd4*Rmd2+Rmd4*Rmd3+Rmd4*Rmd6+Rmd4*Rmd7+Rmd5*Rmd2+Rmd5*Rmd3+Rmd5*Rmd6+Rmd5*Rmd7;
Rn=[(Rmd2+Rmd3+Rmd6+Rmd7)/(a1+eps)-Rmd2/(a1+eps); -Rmd2/(a1+eps)(Rmd1+Rmd2+Rmd4+Rmd5)/(a1+eps);];
Ld=Nv*C*Rn*Nf+eps;
r=;
a2=Ld(1,1)*Ld(2,2)*Ld(3,3)-Ld(1,1)*Ld(2,3)*Ld(3,2)-Ld(2,1)*Ld(1,2)*Ld(3,3)+Ld(2,1)*Ld(1,3)*Ld(3,2)+Ld(3,1)*Ld(1,2)*Ld(2,3)-Ld(3,1)*Ld(1,3)*Ld(2,2)+eps;
L=[1/Lp1,0,0,1/Lp1,0,0;0,1/Lp2,0,0,1/Lp2,0;0,0,1/Lp3,0,0,1/Lp3;0,0,0,-(Ld(2,2)*Ld(3,3)-Ld(2,3)*Ld(3,2))/a2,(Ld(1,2)*Ld(3,3)-Ld(1,3)*Ld(3,2))/a2,(-Ld(1,2)*Ld(2,3)+Ld(1,3)*Ld(2,2))/a2;...
0,0,0,(Ld(2,1)*Ld(3,3)-Ld(2,3)*Ld(3,1))/a2,-(Ld(1,1)*Ld(3,3)-Ld(1,3)*Ld(3,1))/a2,(Ld(1,1)*Ld(2,3)-Ld(1,3)*Ld(2,1))/a2;0,0,0,(-Ld(2,1)*Ld(3,2)+Ld(2,2)*Ld(3,1))/a2,(Ld(1,1)*Ld(3,2)-Ld(1,2)*Ld(3,1))/a2,-(Ld(1,1)*Ld(2,2)-Ld(1,2)*Ld(2,1))/a2;];
dx(1:6)=L*(U1-r*);
dx(7:8)=Rn*Nf*;
dx=dx(:);
end
页:
[1]