在氢气燃烧过程中，有三种成分的浓度值得关注：
- 氧气（反应物的代表）
- 氢原子（中间物的代表，决定了反应速率）
- 水（产物）
根据对各步基元反应的分析，可以得出：反应的整体速率，由最慢的一步反应
\[H+O_2\to OH+O\]
所决定，因此氢原子浓度至关重要。考虑到氢原子的增殖与终止，可以列出一个一阶线性常微分方程表征氢原子浓度的变化：
\[\frac{\mathrm{d}c_H}{\mathrm{d}\tau}=w_H+fc_H-gc_H\]
其中 \(w_H\) 表征了氢原子的初始生成速率，\(f\) 表示了增殖效应，而 \(g\) 表征了销毁效应。
氢原子浓度分析
本问题提法如下：已知 \(\omega_H=0.01\)，\(f=3\)。
- 分别求 \(g=2.9, 3, 3.1\) 的时候，\(c_H\) 随时间变化的曲线。
- \(f=3\)，而 \(t<30\) 时，\(g=2.9\)；\(t\geq30\) 时，\(g=0.1(50-c_H)\)。求 \(c_H\) 随时间变化的曲线。
此处将两部分合而为一，将总计 \(4\) 条曲线绘制在同一幅图中。
解法
课上建议采用最为简便的 Euler 法数值求解这一常微分方程问题。此处打算对这一做法进行最简便的一种优化：将 Euler 法与后退 Euler 法结合，组成一个预估-校正公式。对一个一般的一阶线性方程边值问题
\[\begin{cases}y'(x)=f(x,y)\\y(x_0)=y_0\end{cases}\]
列出如下所示的迭代格式：
\[\begin{cases}p_i=y_i+f(x_i,y_i)\cdot\Delta x\\y_{i+1}=y_i+f(x_{i+1},p_i)\cdot\Delta x\end{cases}\]
其中 \(p_i\) 相当于是用 Euler 法对 \(y_{i+1}\) 做的事先预估，再将其代入到后退 Euler 法中以避免隐式方程的求解。
以下便遵循这个步骤来求解问题。
程序源码与结果
import numpy as np
import matplotlib.pyplot as plt以上引入了必要的 Numpy 库与 Matplotlib 库。
def init():
'初始化各基本参数'
global t, dt, c, w
t = [0]
dt = 0.05
c = [0]
w = .01以上的 init 函数用于在各次模拟之前初始化计算所用的变量与参数。
f = 3 # 增殖因子
g_list = [2.9,3,3.1] # 销毁因子
for num in range(4):
init()
for i in range(750):
t = t + [t[i] + dt]
if num == 3:
if t[i] < 30: g = 2.9
else: g = 0.1 * (50 - c[i])
else: g = g_list[num]
p = c[i] + (0.01 + (f - g) * c[i]) * dt
c = c + [c[i] + (0.01 + (f - g) * p) * dt]
if num == 3: plt.plot(t, c)
else: plt.plot(t, c, linewidth=1, linestyle='--')
plt.title('$c_H$ versus $t$')
plt.xlabel('$t$')
plt.ylabel('$c_H$')
plt.legend(['$g=2.9$', '$g=3.0$', '$g=3.1$', '$g=0.1(50-c_H)$'])
plt.show()从结果可以看出，当 \(g=2.9\) 时生成速率较高，氢原子浓度呈指数上升；当 \(g=3.0\) 时增殖与销毁持平，氢原子仅按初始生成速率线性生成；当 \(g=3.1\) 时销毁速率高于增殖速率，氢原子浓度几乎没有变化。
特别地，若服从问题 (2) 中提出的销毁因子变化规律，氢原子的浓度将在 \(g\) 的突变点处发生转折，迅速降至接近于 \(0\) 的水平。推敲此种规律的提出背景，可能来自于反应物的输入突然停止、反应温度骤降等诸多可能。作为体现总包反应速率的中间产物，氢原子浓度的骤降直接意味着反应的终结。
对规律的修正
仔细推敲氢原子从生成、增殖到全部销毁（尽管从数学上与实在上来说，其都不可能全部销毁）的过程，可以看出：在反应物没有持续供应的情况下，导致氢原子最终散尽的不是 \(g\) 值的上升，而应该是 \(f\) 值的下降——更进一步的说，是反应物的耗尽，导致 \(f\) 从 \(3\) 很快的掉到 \(1\) 左右——与消耗持平，变为不分支的链反应。
为了模拟 \(f\) 值的下降，考虑以下两种规律：
- \(f\) 值在 \(10\) 秒内从 \(3\) 线性的降低到 \(0\)；
- \(f\) 值以 \(\tau=\frac{20}{\mathrm{e}}\) 的时间常数指数下降，趋向于 \(0\) 。（该常数由尝试得出，可参考以下的结果）
以下对这两种过程分别进行模拟。首先是线性下降：
g = 1 # 销毁因子固定为 1
# 线性下降规律
init()
for i in range(1000):
t = t + [t[i] + dt]
if t[i] <= 10: f = 3 - .3 * t[i]
else: f = 0
p = c[i] + (w + (f - g) * c[i]) * dt
c = c + [c[i] + (w + (f - g) * p) * dt]
plt.plot(t, c)
# 指数下降规律
init()
tau = 20 / np.e
for i in range(1000):
t = t + [t[i] + dt]
f = 3 * np.exp(- t[i] / tau)
p = c[i] + (w + (f-g) * c[i]) * dt
c = c + [c[i] + (w + (f - g) * p) * dt]
plt.plot(t, c)
plt.title('$c_H$ versus $t$ (with $f$ decreasing)')
plt.xlabel('$t$')
plt.ylabel('$c_H$')
plt.legend(['$f=0.3(10-t)$', '$f=3\exp(-t\cdot\mathrm{e}/20)$'])
plt.show()从绘图结果可以看出，经过精心取定的两个规律大体上产生了相近的结果：在 \(t=10\) 以内的时间之中，氢原子浓度很快地攀升，又很快的由于反应物的耗尽而回落到接近于 \(0\) 的水平。
产物浓度分析
以上仅分析了中间物——氢原子的浓度变化规律。除此以外，产物——水的浓度变化，也是非常重要的一个指标，值得研究。
根据对各基元反应情况的分析，水的生成速率可按
\[\frac{\mathrm{d}c_{H_2O}}{\mathrm{d}t}=2kc_Hc_{O_2}\]
计算，其中反应常数 \(k\) 在等温条件下维持不变；而 \(c_{O_2}\) 为反应物之一的消耗情况，若我们认为反应时氧气过量，且反应物均按指数规律衰减——与之前氢原子浓度分析时的最后一种情景吻合，则其将很快以指数规律衰减到一个正数（将氢气耗尽后的剩余部分）\(c_{O_2,\text{res}}\)。以下不妨对反应做这样的假设：
- 氢原子浓度按上一节中最后 \(f\) 按 \(\tau=20/\mathrm{e}\) 指数衰减的结果来考虑，不计入中间物与反应物之间的相互影响（用相同的衰变规律表示值这一影响足矣）；
- 氧气浓度按同样的时间常数衰减，设初值为 \(10\)，末值为 \(5\)；
- 常数 \(2k=K\) 设为 \(0.002\)。
接下来据以上假设进行模拟——这种模拟只能推算出产物浓度的总体变化规律，具体的数值并无参考价值。
# 部分参数直接用上面已有的结果
c_w = [0]
for i in range(1000):
c_o = 5 + 5 * np.exp(-t[i] / tau)
c_w = c_w + [c_w[i] + .002 * c[i] * c_o]
plt.title('Concentration of water over time')
plt.xlabel('$t$')
plt.ylabel('$c_{H_2O}$')
plt.plot(t, c_w)
plt.show()可以看到，在这种假设下，水的浓度呈现「由慢转快——稳定生成——由快转慢——停止生成」的增长速率变化规律，这与我们的预期相符合。以上也可以说明，之前对氢原子浓度变化的假设基本合理，应该与某种情形之下的反应条件相吻合——只是具体的细节需要修正与改进。