ODE的计算方法(1)

滴水静禅天•2017-04-09•信息计算• 1340

从最简单的微分方程说起

该方程为最简单系统的一个描述,变量X随时间的变化率与X成比例关系.当给定初始条件IC:t0->X₀,该系统可解为:,可以看出X随时间的变化为相关的指数变化关系.

当>0,呈现指数递增(细菌增长);,恒常速率变化,,呈指数衰减(放射性衰变).

细菌案例的扩展

细菌的生长依赖营养物质的供给(基质S),在单一基质条件下有关系,其中叫做伪计量系数(单位细胞产出数所消费基质的量).同时,也不在恒为常数而有可能是自变量S或者X的函数:

方式	增殖模型	定律
1		零阶(常数)定律
2		Monod定律:基质有限约束
3		Haldane定律:基质有限且免疫约束
4		Contois定律:基质有限且生物抑制约束

对于生物体系:

令:

对于基质体系:

细菌和基质体系微分方程的求解

上述各表达式所需的常数的典型值见下面的参数表:

Matlab处理方法

处理方案以ode45为例:

[t,y] = ode45(odefun,tspan,y0,options)

上式中odefun为求解的目标常微分方程(组),tspan为求解时域定义, y0为初始条件向量, options为计算外加选项.

首先建立odefun.m

function dydt=cellSubSys(t,INPT)

global nu mumax Km Ki Kc kinetics;

X=INPT(1);

S=INPT(2);

switch kinetics

case('const')

mu=mumax;

case('monod')

mu=mumax*(S/(Km+S));

case('haldane')

mu=mumax*(S/(Km+S+(S^2/Ki)));

case('contois')

mu=mumax*(S/Kc*X+S);

end

phi=mu*X;

Xt=phi;

St=-nu*phi;

dydt=[Xt;St];

end

在工作窗口进行调用函数进行计算:

clear all

>> global nu mumax Km Ki Kc kinetics

>> %参数赋值

>> nu=0.5e-11;

>> mumax=1.4;

>> Km=12;

>> Ki=3;

>> Kc=3e-11;

Kinetic={'const','monod','haldane','contois'};

%定义研究时间

>> t=0:0.1:20;

>> %定义边界条件

>> X0=1.4e11;

>> S0=9;

>> %装配输入向量INPT

>> INPT=[X0;S0];

>> %调用求解器计算

>> %%定义求解选项

>> options=odeset('RelTol',1e-3,'AbsTol',1e-3,'Stats','on');

%构建解集准备计算(1-3可以使用ode45,动力学4使用ode15s)

OutSet=cell(3,2);

>> for i=1:3

kinetics=Kinetic{i};

[tout,yout]=ode45(@cellSubSys,t,INPT,options);

OutSet{i,1}=tout;

OutSet{i,2}=yout;

End

%结果后处理可视化

figure(1)

>> for i=1:3

subplot(2,3,i);

plot(OutSet{i,1},OutSet{i,2}(:,1),'r');

subplot(2,3,i+3);

plot(OutSet{i,1},OutSet{i,2}(:,2),'b');

end

发表评论