Question

急!帮我看一下这个matlab程序的问题

是一个微分方程组求解问题
我先建立一个cell.m文件
function dx=cell(t,x)
global d1 b1 k1 k2 k3 d3 b3;
global g4 k4 d4 b5 k5 d5 g6 k6 d6;
global k7 d7 g8 k8 d8 g9 k9 d9 g10 k10 d10 kAb;
dx=zeros(10,1);
dx(1)=x(1)*d1/(x(1)+x(2))^5-d1*x(1)-b1*x(1)*x(3)+k1*x(2)*x(4);
dx(2)=x(2)*d1/(x(1)+x(2))^5-d1*x(2)+b1*x(1)*x(3)-(k1+k2)*x(2)*x(4);
dx(3)=k3*x(2)-d3*x(3)-b3*x(3)*x(7);
dx(4)=(g4+k4*x(2)*x(4)*x(9))*(1-1000*x(4))-d4*x(4);
dx(5)=k5*[d1*x(1)+d1*x(2)+k2*x(2)*x(4)]-d5*x(5);
dx(6)=(g6+k6*x(3)*x(6)*x(10))*(1-10000*x(6))-d6*x(6);
dx(7)=k7*x(6)-d7*x(7)-kAb*b3*x(3)*x(7);
dx(8)=(g8+k8*x(3)*x(8))(1-5000*x(8))-d8*x(8);
dx(9)=[g9+k9*x(3)*(x(8)-g8/d8)]*(1-400*x(9))-d9*x(9);
dx(10)=[g10+k10*x(10)*(x(8)-g8/d8)]*(1-500*x(10))-d10*x(10);
主程序
global d1 b1 k1 k2 k3 d3 b3;
global g4 k4 d4 k5 d5 g6 k6 d6;
global k7 d7 g8 k8 d8 g9 k9 d9 g10 k10 d10 kAb;
d1=0.002;b1=0.01;k1=100;k2=900;b3=0.1;k3=1000;d3=0.4;
g4=1E-10;k4=1E6;d4=0.2;k5=2000;d5=0.25;g6=1E-8;k6=15;
d6=0.4;k7=5.6E4;d7=0.0437;g8=1E-7;k8=10;d8=1.2;
g9=0.5E-10;k9=7;d9=0.2;g10=0.5E-10;k10=3;d10=0.2;kAb=220;V0=1E-7;
x0=[1 0 V0 g4/d4 16 g6/d6 k7*g6/(d6*d7) g8/d8 g9/d9 g10/d10]
t0=[0,300] V0=x0(3)
[t,x]=ode45('cell',t0,x0);
subplot(3,2,1)
plot(t,x(:,2),'-');xlabel('time t');ylabel('感染肝细胞数量Y(t)')
subplot(3,2,2);
plot(t,x(:,3),'-');xlabel('time t');ylabel('病毒数量 V(t)')
subplot(3,2,3)
plot(t,x(:,4),'-');xlabel('time t');ylabel('CTL数量 Z(t)')
subplot(3,2,4)
plot(t,x(:,5),'-');xlabel('time t');ylabel('转氨酶ALT(t)')
subplot(3,2,5)
plot(t,x(:,1),'-');xlabel('time t');ylabel('正常肝细胞数量X(t)')
subplot(3,2,6)
plot(t,x(:,1)+x(:,2),'-');xlabel('time t');ylabel('总肝细胞数量X(t)+Y(t)')
运行结果
??? Error using ==> e:/matlab6p5/toolbox/matlab/funfun/private/odearguments
Solving CELL requires an initial condition vector of length 0.
Error in ==> E:\MATLAB6p5\toolbox\matlab\funfun\ode45.m
On line 155 ==> [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, args, ...
Error in ==> E:\MATLAB6p5\work\cellfangzhen.m
On line 41 ==> [t,x]=ode45('cell',t0,x0);
这是为什么啊

Answer 1

修改了几处：
【1】主程序前加function hahaha，子程序移到后面去了。
【2】用ode23t比ode45快很多
【3】% dx=zeros(10,1)，不用了
【4】dx=[dx1;dx2;....dx10]
最后程序通了。

function hahha
global d1 b1 k1 k2 k3 d3 b3;
global g4 k4 d4 k5 d5 g6 k6 d6;
global k7 d7 g8 k8 d8 g9 k9 d9 g10 k10 d10 kAb;
d1=0.002;b1=0.01;k1=100;k2=900;b3=0.1;k3=1000;d3=0.4;
g4=1E-10;k4=1E6;d4=0.2;k5=2000;d5=0.25;g6=1E-8;k6=15;
d6=0.4;k7=5.6E4;d7=0.0437;g8=1E-7;k8=10;d8=1.2;
g9=0.5E-10;k9=7;d9=0.2;g10=0.5E-10;k10=3;d10=0.2;kAb=220;V0=1E-7;
x0=[1 0 V0 g4/d4 16 g6/d6 k7*g6/(d6*d7) g8/d8 g9/d9 g10/d10]
t0=[0,300],%V0=x0(3)
[t,x]=ode23t(@cell,t0,x0);
subplot(3,2,1)
plot(t,x(:,2),'-');xlabel('time t');ylabel('感染肝细胞数量Y(t)')
subplot(3,2,2);
plot(t,x(:,3),'-');xlabel('time t');ylabel('病毒数量 V(t)')
subplot(3,2,3)
plot(t,x(:,4),'-');xlabel('time t');ylabel('CTL数量 Z(t)')
subplot(3,2,4)
plot(t,x(:,5),'-');xlabel('time t');ylabel('转氨酶ALT(t)')
subplot(3,2,5)
plot(t,x(:,1),'-');xlabel('time t');ylabel('正常肝细胞数量X(t)')
subplot(3,2,6)
plot(t,x(:,1)+x(:,2),'-');xlabel('time t');ylabel('总肝细胞数量X(t)+Y(t)')

function dx=cell(t,x)
global d1 b1 k1 k2 k3 d3 b3;
global g4 k4 d4 b5 k5 d5 g6 k6 d6;
global k7 d7 g8 k8 d8 g9 k9 d9 g10 k10 d10 kAb;
%dx=zeros(10,1);
dx1=x(1)*d1/(x(1)+x(2))^5-d1*x(1)-b1*x(1)*x(3)+k1*x(2)*x(4);
dx2=x(2)*d1/(x(1)+x(2))^5-d1*x(2)+b1*x(1)*x(3)-(k1+k2)*x(2)*x(4);
dx3=k3*x(2)-d3*x(3)-b3*x(3)*x(7);
dx4=(g4+k4*x(2)*x(4)*x(9))*(1-1000*x(4))-d4*x(4);
dx5=k5*[d1*x(1)+d1*x(2)+k2*x(2)*x(4)]-d5*x(5);
dx6=(g6+k6*x(3)*x(6)*x(10))*(1-10000*x(6))-d6*x(6);
dx7=k7*x(6)-d7*x(7)-kAb*b3*x(3)*x(7);
dx8=(g8+k8*x(3)*x(8))*(1-5000*x(8))-d8*x(8);
dx9=[g9+k9*x(3)*(x(8)-g8/d8)]*(1-400*x(9))-d9*x(9);
dx10=[g10+k10*x(10)*(x(8)-g8/d8)]*(1-500*x(10))-d10*x(10);
dx=[dx1;dx2;dx3;dx4;dx5;dx6;dx7;dx8;dx9;dx10]

当t=170时，发生阶跃，不知对不对？