反应扩散方程

反应扩散方程（reaction–diffusion equation）是一类同时描述局部生成、消耗和空间扩散的偏微分方程。浓度较高处的物质向周围输运，改变的是物质的位置；化学反应、种群繁殖或死亡改变的则是当地的数量。这两种过程可以相互促进，也可以产生单独考虑其中一项时没有的传播和失稳现象。

例如，一滴能逐渐分解的染料进入静止液体后，颜色会向外散开，同时总染料量减少。扩散速度与分解速度是两个不同参数：前者控制染色区域变宽的速率，后者控制总量衰减的速率。下面先把这个例子算清楚，再讨论具有增长和多组分耦合的模型。

局部守恒与方程

设 $u(x,t)$ 为位置 $x\in \Omega \subset {ℝ}^n$、时间 $t\ge 0$ 的浓度，$D>0$ 是常扩散系数。常见的标量方程为 $$\frac{\partial u}{\partial t}=D\Delta u+f(u,x,t),\Delta u=\sum _{j=1}^n\frac{{\partial }^{2}u}{\partial x_j^2}.$$ 拉普拉斯算子 $\Delta$ 衡量一个点的值相对邻域的弯曲程度。在一维光滑局部极大点，$u_{xx}\le 0$；严格向下弯曲、二阶导数为负时，扩散项使该点的值降低。局部极小点的符号相反。反应项 $f$ 是单位时间内的局部净生成量，正值表示生成，负值表示消耗。若长度和时间分别以米、秒计，$D$ 的单位为 ${\mathrm{m}}^{\mathrm{2}}/\mathrm{s}$，而 $f$ 的单位为“浓度每秒”。

令 $j$ 为通量，即单位面积、单位时间通过的物质量。采用菲克定律 $j=-D\mathrm{\nabla }u$，其中梯度 $\mathrm{\nabla }u$ 指向浓度增长最快的方向。对任意固定区域 $V$，质量平衡为 $$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{V}u\mathrm{d}x=-{\int }_{\partial V}j\cdot n\mathrm{d}S+{\int }_{V}f\mathrm{d}x.$$ 这里 $n$ 是外法向。第一项减去流出的量，第二项加上内部生成的量。利用散度定理把边界通量化为体积分，得到 $u_t=-\mathrm{\nabla }\cdot j+f=D\Delta u+f$。这说明方程来自一个明确的收支关系。反应扩散模型及其生物学应用可参见 Byrne 的讲义。

在有界区域中，还要给出初值 $u(x,0)=u_0(x)$ 和边界条件。封闭容器常用无通量条件 $\mathrm{\nabla }u\cdot n=0$；与保持固定浓度的外界接触时，可指定边界上的 $u$。无通量只消除了边界收支，此时仍有 $$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }u\mathrm{d}x={\int }_{\Omega }f\mathrm{d}x.$$ 所以发生反应时，封闭容器中的总量未必守恒。

扩散与分解的完整算例

为把空间分散与总量减少分开，考虑整条实线上的模型 $$u_t=D{u}_{xx}-ku,u(x,0)=e^{-x^2},D>0,k\ge 0.$$ 这里已按初始宽度和时间尺度把变量无量纲化；它是一组教学参数，不对应实测染料。$k$ 是线性分解系数。初始曲线是高斯峰，试求保持高斯形状的解 $$u(x,t)=A(t)e^{-B(t)x^2},A(0)=B(0)=1.$$ 逐项求导，除以正的 $u$，可得 $$\frac{u_t}{u}=\frac{A^{\prime }}{A}-B^{\prime }{x}^2,\frac{u_{xx}}{u}=-2B+4B^2{x}^2.$$ 比较常数项和 $x^2$ 项，得到两个常微分方程： $$B^{\prime }=-4D{B}^2,\frac{A^{\prime }}{A}=-2DB-k.$$ 第一个方程给出 $(1/B{)}^{\prime }=4D$，所以 $B=1/(1+4Dt)$。代入第二个并从零积分，得 $$\log A=-kt-\frac{1}{2}\log (1+4Dt).$$ 于是 $$u(x,t)=\frac{e^{-kt}}{\sqrt{1+4Dt}}\exp (-\frac{x^2}{1+4Dt}).$$ 令 $t=0$ 恢复初值；将刚得到的 $A,B$ 代回导数式也直接验证了方程。这一解是热方程的高斯扩散与指数衰减的乘积。

取 D=0.2、k=0.3。曲线变宽由扩散控制，总面积减少由分解控制，峰值降低同时包含两种作用。

在图中取 $D=0.2,k=0.3$。时间 $t=1$ 时，中心值为 $e^{-0.3}/\sqrt{1.8}\approx 0.5522$。但峰值本身不能衡量总量。把解对全空间积分，并令 $z=x/\sqrt{1+4Dt}$，得到 $$M(t)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}u(x,t)\mathrm{d}x=\sqrt{\pi }{e}^{-kt}.$$ 归一化浓度 $u/M$ 的方差是 $(1+4Dt)/2$。因此 $t=1$ 时，总量约为初始的 $0.7408$，方差则从 $0.5$ 增到 $0.9$。若只把 $k$ 改成零，总量保持 $\sqrt{\pi }$，峰依然降低。这是区分扩散和消耗的一种直接检验。

局部增长与传播前沿

当局部反应是受资源限制的增长时，常用 Fisher–KPP 方程 $$u_t=D{u}_{xx}+ru(1-\frac{u}{K}),D,r,K>0.$$ $r$ 是低密度增长率，$K$ 是局部承载水平。如果初值空间均匀，则扩散项为零，方程退化为 Logistic 增长；非均匀初值则允许已占据区域向外推进。Fisher 与 Kolmogorov、Petrovsky、Piskunov 于 1937 年分别研究了这类传播问题。

设向右传播的行波为 $u(x,t)=U(z)$，$z=x-ct$、$c>0$，连接 $U(-\mathrm{\infty })=K$ 与 $U(+\mathrm{\infty })=0$。由链式法则，$u_t=-c{U}^{\prime }$、$u_{xx}=U^{″}$，故 $$D{U}^{″}+c{U}^{\prime }+rU(1-U/K)=0.$$ 在前沿右侧 $U$ 很小时，先略去二次项，再代入衰减形式 $U\approx e^{-\lambda z}$，得到 $$D{\lambda }^2-c\lambda +r=0.$$ 具有正实衰减指数要求判别式非负，即 $c\ge 2\sqrt{Dr}$。这个计算解释了最小波速的来源；它只分析前沿附近，完整的连接轨道存在性还需要行波方程的相平面论证。标准 Fisher–KPP 方程确实在这些波速存在单调行波。

在 $D=r=K=1$ 的无量纲情形，有一个能逐项验证的精确例子： $$U(z)=\frac{1}{(1+e^{z/\sqrt{6}}{)}^2},c=\frac{5}{\sqrt{6}}.$$ 令 $w=e^{z/\sqrt{6}}$，则 $$U^{\prime }=-\frac{2w}{\sqrt{6}(1+w{)}^3},U^{″}=\frac{w(2w-1)}{3(1+w{)}^4}.$$ 把 $U^{″}$、$c{U}^{\prime }$、$U(1-U)$ 通分，分子依次为 $(2w^2-w)/3$、$-5w(1+w)/3$、$2w+w^2$，相加为零。左右无穷远的极限也分别是 $1,0$，所以这确实是一条传播前沿。

精确行波在 t=0、1、2 的位置。这里 c=5/√6≈2.0412，是一个特定波速；最小波速为 2。各曲线在 x=ct 处的高度为 1/4。

多组分系统与扩散驱动失稳

一个组分中扩散倾向于抹平差异，多组分之间的反馈却可能使某些差异增长。设 $𝐮=(u,v{)}^{𝖳}$，系统为 $${𝐮}_t=\left(\begin{array}{cc}{D}_1& 0\\ 0& D_2\end{array}\right)\Delta 𝐮+𝐟(𝐮).$$ 在均匀平衡 ${𝐮}_{\ast }$ 附近令扰动为 $𝜼$，只保留一阶项，得到 $${𝜼}_t=\mathrm{diag}(D_1,D_2)\Delta 𝜼+J𝜼,J={\frac{\partial 𝐟}{\partial 𝐮}|}_{{𝐮}_{\ast }}.$$ 矩阵 $J$ 记录两个组分的小变化如何彼此影响。若一个空间模式满足 $\Delta \varphi =-q\varphi$，取 $𝜼=𝐚(t)\varphi (x)$ 后，其振幅满足 $${𝐚}^{\prime }=(J-q\mathrm{diag}(D_1,D_2))𝐚.$$ 因此应考察这个矩阵的特征值，而不是只检查 $J$。

取一组用于计算的局部线性化参数 $$J=\left(\begin{array}{cc}1& -2\\ 2& -3\end{array}\right),D_1=1,D_2=10.$$ 不考虑空间变化时，$J$ 的特征多项式为 $(\lambda +1{)}^2$，所有线性小扰动最终趋于零，即均匀平衡渐近稳定；这里不要求扰动范数从初始时刻起单调下降。加入波数平方 $q$ 后，模式矩阵的迹和行列式变为 $$\mathrm{tr}=-2-11q,det=(1-q)(-3-10q)+4=1-7q+10q^2.$$ 行列式可分解为 $(1-5q)(1-2q)$。当 $1/5<q<1/2$ 时它为负，两个实特征值异号，于是有一个正的增长率。原来稳定的局部系统因此对某些空间变化失稳，这称为扩散驱动失稳，也是 Turing 斑图理论的线性机制。

最大增长率在 q∈(0.2,0.5) 为正。图只放大零附近的纵轴范围；圆点对应长度 2π、无通量边界下可实现的 q=1/4 模式。

边界会限制可出现的波数。在区间 $[0,2\pi ]$ 上取无通量边界，模式为 $\cos (nx/2)$，$q_n=n^2/4$，其中 $n=0,1,2,\dots$。只有 $n=1$ 的 $q_1=1/4$ 落在失稳区间，其正特征值为 $${\lambda }_{+}=\frac{-19+\sqrt{369}}{8}\approx \mathrm{0.02617.}$$ 所以该模式的小振幅最初近似按 $e^{0.02617t}$ 增长。这里尚未指定完整的非线性反应函数；线性计算不能决定增长后的最终形状，更不能单独保证出现稳定斑图。相比之下，单组分常平衡若满足 $f^{\prime }(u_{\ast })<0$，各模式增长率为 $f^{\prime }(u_{\ast })-Dq<0$，正扩散不会使它以同样机制失稳。

模型边界

常系数扩散假设介质在所研究尺度上可用连续、近似均匀的输运规律表示。若存在流动，应加入对流；若扩散系数随位置或浓度变化，守恒形式通常是 $\mathrm{\nabla }\cdot (D\mathrm{\nabla }u)$，不能直接写成 $D\Delta u$。有些反应项还可能造成有限时间爆破，因此仅写出“扩散加反应”并不自动保证解始终有界。实际建模时，反应规律、边界供给和参数单位需要分别核实。