简单的常微分方程(组)初值问题

问题

A，B，C，D，E 五个化合物，A 反应后生成 B，B 反应后生成C，C 反应后生成 D，D 反应后生成 E，即

A—→B—→C—→D—→E

反应动力学表明都是一级反应，已知 k1，k2，k3，k4 是相应的反应速率常数。A，B，C，D，E 的浓度随时间变化的反应速率方程如下：

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪1.2.3.4.5.dCAdt=−k1⋅CAdCBdt=k1⋅CA−k2⋅CBdCCdt=k2⋅CB−k3⋅CCdCDdt=k3⋅CC−k4⋅CDdCEdt=k4⋅CD

当时间 t=0 时，A 的初始浓度为 A0，而 B，C，D，E 的初始浓度都为 0。

求A,B,C,D,E 浓度随时间变化的函数。

解答

简单的常微分方程（组）初值问题。可以顺次求解，也可以统一成方程组求。因为计算比较繁琐，用Maple或mathematica计算比较方便。在微分方程符号解方面，总体上看Maple比Mathematica胜出一筹。但在当前这个问题上，两种软件求解都易如反掌。

Mathematica顺次求解ODE初值问题的代码如下：

ClearAll["Global`*"]a[t_]:=Evaluate[a[t]/.(DSolve[{a'[t]==-k1 a[t],a[0]==A0},a,t]//FullSimplify)[[1,1]]]b[t_]:=Evaluate[b[t]/.(DSolve[{b'[t]==k1*a[t]-k2 b[t],b[0]==0},b,t]//FullSimplify)[[1,1]]]c[t_]:=Evaluate[c[t]/.(DSolve[{c'[t]==k2*b[t]-k3 c[t],c[0]==0},c,t]//FullSimplify)[[1,1]]]d[t_]:=Evaluate[d[t]/.(DSolve[{d'[t]==k3*c[t]-k4 d[t],d[0]==0},d,t]//FullSimplify)[[1,1]]]e[t_]:=Evaluate[e[t]/.(DSolve[{e'[t]==k4*d[t],e[0]==0},e,t]//FullSimplify)[[1,1]]]

得到结果：

CA(t)=CB(t)=CC(t)=CD(t)=CE(t)=A0e−k1tA0k1e−k1t(e(k1−k2)t−1)k1−k2A0k1k2(e−k3t(k1−k2)−e−k2t(k1−k3)+e−k1t(k2−k3))(k1−k2)(k1−k3)(k2−k3)A0k1k2k3(e−k4t(k1−k4)(k4−k2)(k4−k3)+e−k3t(k1−k3)(k2−k3)(k4−k3)+e−k2t(k1−k2)(k3−k2)(k4−k2)+e−k1t(k1−k2)(k1−k3)(k4−k1))A0(k2k3k4e−k1t(k1−k2)(k1−k3)(k1−k4)−k1k3k4e−k2t(k1−k2)(k2−k3)(k2−k4)−k1k2k4e−k3t(k1−k3)(k3−k2)(k3−k4)−k1k2k3e−k4t(k1−k4)(k4−k2)(k4−k3))

最长的 CD(t) 只好贴出来：

CE(t):

这五个浓度之和是常数 A0 ，它们导数之和为 0 可以帮助验算。

思考

这里得到的解析解，乍一看似乎很方便。实际上，如果 ki=kj,(i≠j∈{1,2,3,4}) ，似乎出现了分母为 0 的情形。比如 k2=k3 时…… 难道是该动力学方程要求动力学常数不能相等吗？

如何解决？从极限的角度，或者，直接通过微分方程特定条件下求解的方式