MATLAB 相关微分方程系统

Question

给定一个微分方程组，例如：

dy/dt = f(t)
dx/dt = g(t)

使用dsolve可以找到解决方案，如：

dsolve(diff(y) == f(t), diff(x) == g(t), y(0) == 1, x(0) == 1);

但是对于所有变量相互依赖的系统呢：

dy/dt = f(y,z)
dx/dt = g(x,y)
dz/dt = h(z,x)

当以相同的方式处理初始条件时，对于确实有解的系统，我找不到解。

我知道我尝试过的系统可以产生解决方案，因为我使用了随机/确定性模拟器 - 认为可能使用了一些奇怪的语法。

如果有帮助，我正在专门寻找导数都为零的解决方案。

编辑：

这是一个例子：

PX/dt = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PZ^n))))/kd_mRNA)-kd_prot*PX;
PY/dt = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PX^n))))/kd_mRNA)-kd_prot*PY;
PZ/dt = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PY^n))))/kd_mRNA)-kd_prot*PZ;

系数：

eff = 20;

KM = 40;

tau_mRNA=2.0;
tau_prot=10;

ps_a=0.5;
ps_0=5.0E-4;

t_ave = tau_mRNA/log(2);

k_tl=eff/t_ave;

a_tr=(ps_a-ps_0)*60;
a0_tr=ps_0*60;

kd_mRNA = log(2)/tau_mRNA;
kd_prot = log(2)/tau_prot;

beta = tau_mRNA/tau_prot; 
alpha = a_tr*eff*tau_prot/(log(2)*KM);
alpha0 = a0_tr*eff*tau_prot/(log(2)*KM);
n=2;

以及初始条件：

PX0 = 20;
PY0 = 0;
PZ0 = 0;

这会产生一个响应：

这显然有一个稳态解（所有导数均为 0）。

在 MATLAB 中我尝试过：

%%
syms PX(t) PY(t) PZ(t);

z = dsolve(diff(PX) == (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PZ^n))))/kd_mRNA)-kd_prot*PX, diff(PY) == (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PX^n))))/kd_mRNA)-kd_prot*PY, diff(PZ)==(k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PY^n))))/kd_mRNA)-kd_prot*PZ,PX(0)==20)

和：

%%
eq1 = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PZ^n))))/kd_mRNA)-kd_prot*PX;
eq2 = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PX^n))))/kd_mRNA)-kd_prot*PY;
eq3 = (k_tl*(a0_tr + ((a_tr*KM^n)/((KM^n) + (PY^n))))/kd_mRNA)-kd_prot*PZ;

dsolve(diff(PX)==eq1,PX(0)==20,diff(PY)==eq2,PY(0)==0,diff(PZ)==eq3,PZ(0)==0)

两者都不产生错误，但返回一个空符号。

Answer 1

您的数值解似乎有一个振荡分量。 “稳态”可以是零幅限循环，这是一个不平凡的解决方案。您绝对不应该期望这样的系统具有易于找到的分析解决方案。您的三个变量之间的循环关系也无济于事。无论如何，Mathematica 10 的 DSolve 也无法找到解决方案。

虽然它不会让您找到解决方案，但您使用符号数学的方式并不是最佳的。当您在符号数学方程式中使用 log(2) 之类的东西时，应首先将 2 转换为符号值。例如，sym(log(2)) 产生近似值6243314768165359/9007199254740992，而log(sym(2)) 返回精确的log(2)。如果存在，后一种形式更有可能导致解决方案。这是您的代码的修改版本，不幸的是仍然返回“警告：找不到显式解决方案”：

eff = 20;

KM = 40;

tau_mRNA=2;
tau_prot=10;

ps_a=1/sym(2);
ps_0=5/sym(10000);

ln2 = log(sym(2));
t_ave = tau_mRNA/ln2;

k_tl=eff/t_ave;

a_tr=(ps_a-ps_0)*60;
a0_tr=ps_0*60;

kd_mRNA = ln2/tau_mRNA;
kd_prot = ln2/tau_prot;

beta = tau_mRNA/tau_prot; 
alpha = a_tr*eff*tau_prot/(ln2*KM);
alpha0 = a0_tr*eff*tau_prot/(ln2*KM);
n=2;

PX0 = 20;
PY0 = 0;
PZ0 = 0;

syms PX(t) PY(t) PZ(t);

eq1 = (k_tl*(a0_tr + a_tr*KM^n/(KM^n + PZ^n))/kd_mRNA)-kd_prot*PX;
eq2 = (k_tl*(a0_tr + a_tr*KM^n/(KM^n + PX^n))/kd_mRNA)-kd_prot*PY;
eq3 = (k_tl*(a0_tr + a_tr*KM^n/(KM^n + PY^n))/kd_mRNA)-kd_prot*PZ;

s = dsolve(diff(PX,t)==eq1,diff(PY,t)==eq2,diff(PZ,t)==eq3,PX(0)==20,PY(0)==0,PZ(0)==0)