求教一道ode解微分方程的问题

amygod

问题描述是这样的:
dc1/dt=exp(-T)*c1*c2;
dc2/dt=exp(-T)*c1*c2;
dc3/dt=-exp(-T)*c1*c2;
dc4/dt=-exp(-T)*c1*c2;
dT/dt=H*exp(-T)*c1*c2/G;
t从0到1;
关于上式的G和H:
A=[1 2 3 4];
B=[2 3 6 5];
C=[1 5 6 9];
D=[4 5 8 9];
x=[c1 c2 c3 c4]/(c1+c2+c3+c4);
对于一个特定的T有：
G=sum((A+B*T+C*T^2+D*T^3).*x)/T;

H在T=200时为20;
H=20+(积分号)(从200到T)GdT;

初始条件t=0时:
c1=1;c2=1.2;c3=0;c4=0;T=800;
由上面的条件可算得G在T=800时的值G0;
求t=1时c1的值.
这程序怎么编啊,高手帮忙啊~~~~~~~~先谢谢了~!

fanghuikeer

G,H是c1,c2,c3,c4,T的简单函数，我不知道matlab中可不可以直接把G,H关于c1,c2,c3,c4,T的表达式直接算出来，但是最起码手算还是可以的，我大概算了一下，形式还比较简单，建议提问题的人也去亲自计算一下．
下面我给你举例说明你这种问题怎么去编程计算：
function xdot=funfun(t,c)
A=[1 2 3 4];
B=[2 3 6 5];
C=[1 5 6 9];
D=[4 5 8 9];
x=[c(1) c(2) c(3) c(4)]/(c(1)+c(2)+c(3)+c(4));
G=c(3);％这里G我是随意设的一个表达式，因为你的题中这个表达式特别的繁琐．
H=20+c(5);％同上,这里用c(5代替T
xdot=[exp(-c(5))*c(1)*c(2);exp(-c(5))*c(1)*c(2);-exp(-c(5))*c(1)*c(2);-exp(-c(5))*c(1)*c(2);H*exp(-c(5))*c(1)*c(2)/G]

[t,c]=ode45(@funfun,[0 1],[1 1.2 0 0 800])

Warning: Divide by zero.
> In funfun at 9
  In funfun\private\odearguments at 110
  In ode45 at 173
Warning: Divide by zero.
> In funfun at 9
  In ode45 at 324

t =

         0
    0.0250
    0.0500
    0.0750
    0.1000
    0.1250
    0.1500
    0.1750
    0.2000
    0.2250
    0.2500
    0.2750
    0.3000
    0.3250
    0.3500
    0.3750
    0.4000
    0.4250
    0.4500
    0.4750
    0.5000
    0.5250
    0.5500
    0.5750
    0.6000
    0.6250
    0.6500
    0.6750
    0.7000
    0.7250
    0.7500
    0.7750
    0.8000
    0.8250
    0.8500
    0.8750
    0.9000
    0.9250
    0.9500
    0.9750
    1.0000


c =

    1.0000    1.2000         0         0  800.0000
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
       NaN       NaN       NaN       NaN       NaN
具体情况自己去试，之所以本题我计算的无解是因为G,H我没有计算具体表达式
假如输入正确的表达式，一定可以得到想要得解．．．

amygod

可问题的关键就是函数H的表达啊,因为它是带有积分号的,积分上限即为变量T,(也就是c(5));
这在函数xdot该怎么定义啊?请LS指教~