求解微分方程组

牛仔裤

解 微分方程组
dy1/dx=M-(k2-(1-2*r)*x* Q/v)*y1+A1*y2
dy2/dx=B1*y1-B2*y2
请问：1，怎么用runge_kutta的方法求解出它在(0,1000)内的数值解。原程序怎么编？
2，请问步长能设为1吗？
3，如果不用runge_kutta的方法，用其他方法也可以。
谢谢 帮忙。

[ 本帖最后由 xinyuxf 于 2007-7-22 12:10 编辑 ]

xjzuo

请先给出参数值及初始条件,以方便他人调试.
另:建议自己先动手写一写.

牛仔裤

解 微分方程组
dy1/dx=M-(k2-(1-2*r)*x* Q/v)*y1+A1*y2
dy2/dx=A1*y1-B2*y2
其参数为：M＝0.000002；A1＝0.95；考＝5.0；B2＝1；Q＝0.00005；V＝0.0251；r=1;
初始条件：x=0;y1=20000;y2=300000;
请问：1，怎么用runge_kutta的方法求解出它在(0,1000)内的数值解。原程序怎么编？
2，请问步长能设为1吗？
3，如果不用runge_kutta的方法，用其他方法也可以。
谢谢 帮忙:lol


======================
用ode45可以轻易解决,自己动动手吧.
by xjzuo
======================

[ 本帖最后由 xjzuo 于 2007-3-23 08:51 编辑 ]

牛仔裤

实在不好意思，请问ode45是什么？

eight

help ode45
这是matlab最常用的命令，请紧记

[ 本帖最后由 ChaChing 于 2010-4-2 23:24 编辑 ]

牛仔裤

解 微分方程组
dy1/dx=M-(k2-(1-2*r)*x* Q/v)*y1+A1*y2
dy2/dx=B1*y1-B2*y2
请问：1，怎么用runge_kutta的方法求解出它在(0,1000)内的数值解。原程序怎么编？
2，请问步长能设为1吗？
3，如果不用runge_kutta的方法，用其他方法也可以。
这是我编的程序，但是画出的图很不合理，请指点一下？是不是解的思路不对，我就是用runge_kutta的方法解的。
下面是程序
k1=0.95;k2=5.0;k3=0.5
Q0=0.00005;v1=0.0251;c0=0.001;r=1;n=1
M=Q0*c0/v1;A1=k1;B2=k1+k3;c1=zeros(1,8001);c2=zeros(1,8001);c1(1,1)=20000;c2(1,1)=300000;
r1=zeros(4,201);r2=zeros(4,201);h=0.1;t=1:8001;
for i=1:h:800
    A2=(k2+(1-2*r)*(i+h/2)*Q0/v1)
    A22=(k2+(1-2*r)*(i+h)*Q0/v1)
    r1(1,n)=M-A2*c1(n)+A1*c2(n);
    r2(1,n)=A1*c1(n)-B2*c2(n);
    r1(2,n)=M-A2*(c1(n)+h/2*r1(1,n))+A1*(c2(n)+h/2*r2(1,n));
    r2(2,n)=A1*(c1(n)+h/2*r1(1,n))-B2*(c2(n)+h/2*r2(1,n));
    r1(3,n)=M-A2*(c1(n)+h/2*r1(2,n))+A1*(c2(n)+h/2*r2(2,n));
    r2(3,n)=A1*(c1(n)+h/2*r1(2,n))-B2*(c2(n)+h/2*r2(2,n));
    r1(4,n)=M-A22*(c1(n)+h*r1(3,n))+A1*(c2(n)+h*r2(3,n));
    r2(4,n)=A1*(c1(n)+h*r1(3,n))-B2*(c2(n)+h*r2(3,n));
    c1(n+1)=c1(1,n)+h/6*(r1(1,n)+2*r1(2,n)+2*r1(3,n)+r1(4,n));
    c2(n+1)=c2(1,n)+h/6*(r2(1,n)+2*r2(2,n)+2*r2(3,n)+r2(4,n));
    n=n+1
end
plot(t,c1,'r+');hold on;xlabel('时间/d'),ylabel('渗滤液BOD浓度/(mg/L)')

xjzuo

这样检查比较累,建议直接用ode45求解．
随便取参数画了一下，不知是否是你要的效果．

牛仔裤

谢谢 楼上的答复，这个不是我要画的图，我感觉的我的计算有问题，咳！！