求常微分方程初值问题

yxlnbu

大家好，我有一个问题，就是求解一个常微分方程组含2个变量r,vr分别是时间t的函数，只是vr的表达式有点复杂。t从0到10^-4s,
r和vr的初始值分别为1.60，和0
我编制解法时是用Simpson规则来的，可是算出来就是不正确，请各位指点。谢谢了。下面是我的程序:
function dr=myfun_1(t,r)
br=-( -4.6861.*r(1).^5 + 72.6026.*r(1).^4 - 447.6495.*r(1).^3 + 1376.6343.*r(1).^2-2121.6799.*r(1)+1323.3641);
I2=7.05e19.*(exp(-1.92e5.*t)).*sin(t.*1.57e5);
I1=-3.60e19.*(exp(-1.92e5.*t)).*sin(t.*1.57e5);
A=0.01;
c=3e10;
ro=8.924;
Y0=69;
K=36;
m=0.45;
mr=0.1842;
dr=zeros(2,1);
dr(1)=r(2);
dr(2)=(br.*I1.*I2./(ro.*A.*c^2))+(((4*pi/c^2)*(log(8*r(1)/ 0.0564)-1.8064)).*I2^2./(2.*mr))-(Y0.*(1+K.*dr(2)./r(1)).^m)./ro.*r(1);
以上是函数方程。另外我也尝试过把红色的字代入到方程中。可是一样还是解不出来
[t,x]=ode45(@myfun_1,[0,1e-4],[1.6;0]) 这个是解法
想请大侠指教，以上表示是否正确，如果不正确应该如何呢。对于这么复杂的函数是否是编程错误？

meiyongyuandeze

没有报错信息吗?

yxlnbu

回复 2 # meiyongyuandeze 的帖子

没有报错。我是两个分开来运行的。就是在看x值的时候，第一行有初始值，可是下面的都是NaN

meiyongyuandeze

回复 3 # yxlnbu 的帖子

我看了你的方程，其中某些项的差别比较大，应该属于刚性方程的范畴，不宜采用ode45，而应选取ode15s来求解，我求解了下，没问题，有收敛解，只不过计算时间比较长。

1. [t,x]=ode15s(@myfun_1,[0,1.0e-4],[1.6;0]);

复制代码

yxlnbu

回复 4 # meiyongyuandeze 的帖子

谢谢您，我后来在纸上也测算了下，最后一项的数量级相差很大有e10左右，才开始用matlab，不懂的地方有很多，还请多指点了。（运行速度确实慢……数据非常多，电脑比较差，差不多matlab不能用了，内存都到2G了……）

yxlnbu

回复 4 # meiyongyuandeze 的帖子

你好，我把方程改了下，变成有个四个方程的常微分方程组，定义如下，该方程可以求解，可是解基本上都是NAN,我想请您帮我看看，是什么问题。我想得到的是时间从e-6—e-4这一段的数值
谢谢您了
function xdot=dianci(t,x)
A=0.01;
c=3e10;
ro=8.924;
Y0=69e10;
K=36;
rw=0.05642;
Rs=0.004;
Rc=4.73e14;
Lc=1.542e-18;
U=13.34256384;
br=-( -4.6861.*x(1).^5 + 72.6026.*x(1).^4 - 447.6495.*x(1).^3 + 1376.6343.*x(1).^2-2121.6799.*x(1)+1323.3641);
mr=(396.51-413.32.*x(1)+192.46.*x(1).^2-46.45.*x(1).^3+5.69.*x(1).^4-0.28.*x(1).^5);
L=(4*pi/c^2)*(log(8*x(1)/rw)-1.75);
dl=(4*pi/c^2)*(log(8*x(1)/rw)-0.75);%dl/dr
M=(mr.*2.*pi./c.^2);
dm=((2.*pi./c.^2).*br.*x(1));%dm/dr
xdot=[ x(2);%r
      br.*x(3).*x(4)./(ro.*A.*c^2)+(x(4)^2./(4.*pi.*A.*x(1)))*dl-(Y0.*(1+K.*x(2)./x(1)).^0.45)./(ro.*x(1));%vr
      (1/(M^2-L*Lc))*(L*(Rc.*x(3)+x(4).*x(2).*dm+13.33)-M*(Rs.*x(4)+x(3).*x(2).*dm+x(4)*x(2).*dl));%i1
      (1/(M^2-L*Lc))*(Lc*(Rs.*x(4)+x(3).*x(2).*dm+x(4).*x(2).*dl)-M.*(Rc.*x(3)+x(4).*x(2).*dm+U));%i2
      ];%函数方程

[t y]=ode45(@dianci,[0,1],[1.6;0;0;0;])%解法及初值

meiyongyuandeze

变步长积分有些时间点不能达到收敛的条件！不知道积分步长怎么会这么小啊！最好还是用ode15s吧！

yxlnbu

回复 7 # meiyongyuandeze 的帖子

您好，我把初值改了下，原先的初值有问题，我把t=0代入方程求出的结果 当作初值了。这么大意就把错误的发上来，麻烦你了。我想这个问题是不是可以从函数整体上考虑，把时间扩大e4呢？这样来说可能步长会大一点