数值抛物型方程的初值问题
我们先从抛物型方程的基本概念开始。抛物型偏微分方程描述的是诸如热传导、扩散等物理过程，其最典型的代表是热传导方程：
\[ \frac{\partial u}{\partial t} = a \frac{\partial^2 u}{\partial x^2} \]
其中 \(u(x, t)\) 是未知函数（例如温度），\(a > 0\) 是扩散系数。一个完整的抛物型方程问题通常包含方程本身、初始条件和边界条件。我们首先聚焦于初值问题，也称为柯西问题。在这个问题中，我们考虑在整个空间（例如一维的整条实数轴）上的演化问题，其定解条件为：
- 方程: \(u_t = a u_{xx}\), 其中 \(x \in (-\infty, \infty), \, t > 0\)
- 初始条件: \(u(x, 0) = u_0(x)\), 其中 \(x \in (-\infty, \infty)\)
这里的核心是，在初始时刻 \(t=0\)，我们已知函数 \(u\) 在整个空间上的分布 \(u_0(x)\)，目标是求解未来任意时刻 \(t>0\) 的函数 \(u(x,t)\)。
基本解与解的积分表示
对于上述在无限区域上的初值问题，我们可以求出其解析解，这为我们理解和检验数值方法提供了基础。这个解可以通过傅里叶变换等方法求得，其形式是一个卷积积分：
\[ u(x, t) = \int_{-\infty}^{\infty} \Phi(x - y, t) \, u_0(y) \, dy \]
其中，核函数 \(\Phi(x, t)\) 称为热传导方程的基本解（或高斯核、热核）：
\[ \Phi(x, t) = \frac{1}{\sqrt{4\pi a t}} \exp\left( -\frac{x^2}{4a t} \right) \quad (t > 0) \]
这个解的物理意义非常直观：在任意点 \(x\) 和时刻 \(t\) 的解 \(u(x,t)\)，是初始时刻所有点 \(y\) 的初始值 \(u_0(y)\) 的加权平均。权重由基本解 \(\Phi(x-y, t)\) 给出，它描述了一个从 \(y\) 点出发的“热源”经过时间 \(t\) 后对 \(x\) 点产生的“影响”。随着时间 \(t\) 增大，基本解 \(\Phi\) 会变得越来越“平坦”和“宽阔”，这反映了扩散过程的抹平效应。
数值离散的核心思想：有限差分法
在实际问题中，初始条件 \(u_0(x)\) 可能很复杂，或者区域不是无限的，导致解析解难以求得。这时我们需要数值方法。最直接的方法是有限差分法。
数值求解初值问题的核心思想是进行离散化：
- 空间离散: 虽然理论上是无限区域 \((-\infty, \infty)\)，但数值计算必须在有限区域内进行。我们选择一个足够大的计算区域 \([X_l, X_r]\)，并在此区域上布置网格点 \(x_j = X_l + j\Delta x\), \(j = 0, 1, \dots, J\)，其中 \(\Delta x\) 是空间步长。
- 时间离散: 将时间也离散化，\(t^n = n\Delta t\), \(n = 0, 1, 2, \dots\)，其中 \(\Delta t\) 是时间步长。
我们的目标是求解网格点 \((x_j, t^n)\) 上的近似解 \(U_j^n \approx u(x_j, t^n)\)。
经典差分格式及其推导
我们用差商来近似偏导数。
- 时间导数： \(\frac{\partial u}{\partial t}(x_j, t^n) \approx \frac{U_j^{n+1} - U_j^n}{\Delta t}\)（向前差商）
- 空间二阶导数： \(\frac{\partial^2 u}{\partial x^2}(x_j, t^n) \approx \frac{U_{j+1}^n - 2U_j^n + U_{j-1}^n}{(\Delta x)^2}\)（中心差商）
将这两个近似代入原方程 \(u_t = a u_{xx}\)，我们得到最经典的显式欧拉格式（或称向前欧拉格式）：
\[ \frac{U_j^{n+1} - U_j^n}{\Delta t} = a \frac{U_{j+1}^n - 2U_j^n + U_{j-1}^n}{(\Delta x)^2} \]
整理后可得：
\[ U_j^{n+1} = U_j^n + \frac{a\Delta t}{(\Delta x)^2} (U_{j+1}^n - 2U_j^n + U_{j-1}^n) \]
这个格式的优点是显式的，意味着 \(U_j^{n+1}\) 可以直接由上一时间层 \(t^n\) 的已知值 \(U_{j-1}^n, U_j^n, U_{j+1}^n\) 简单算术运算得到，无需求解方程组。
稳定性分析：CFL条件的抛物型版本
显式格式虽然简单，但有一个致命的限制：稳定性条件。如果时间步长 \(\Delta t\) 和空间步长 \(\Delta x\) 的选取不满足特定关系，计算误差会指数级增长，导致解完全失真。
对于抛物型方程的显式欧拉格式，其稳定性条件是：
\[ \mu := \frac{a\Delta t}{(\Delta x)^2} \leq \frac{1}{2} \]
其中 \(\mu\) 称为网格傅里叶数或扩散数。
这个条件的物理和数学含义是：信息（扰动）在一个时间步长 \(\Delta t\) 内的扩散范围，必须被单个网格单元所“捕获”。如果 \(\Delta t\) 过大，一个时间步内从 \(x_j\) 点扩散出去的信息会超出其相邻网格点 \(x_{j\pm1}\) 的范围，导致算法无法正确描述物理过程，从而失稳。
隐式格式：无条件稳定的选择
为了克服显式格式的严格稳定性限制，我们可以采用隐式格式。最经典的是隐式欧拉格式（或称向后欧拉格式），它用 \(t^{n+1}\) 时刻的空间导数来近似方程：
\[ \frac{U_j^{n+1} - U_j^n}{\Delta t} = a \frac{U_{j+1}^{n+1} - 2U_j^{n+1} + U_{j-1}^{n+1}}{(\Delta x)^2} \]
整理后：
\[ -\mu U_{j-1}^{n+1} + (1 + 2\mu) U_j^{n+1} - \mu U_{j+1}^{n+1} = U_j^n \]
这个格式是隐式的，因为未知量 \(U^{n+1}\) 同时出现在方程两边。为了求解 \(n+1\) 时间层上所有网格点 \(j\) 的值，我们需要求解一个线性方程组。这个方程组是三对角的，可以用高效算法（如托马斯算法）求解。
隐式欧拉格式最大的优点是无条件稳定，即对于任何 \(\Delta t\) 和 \(\Delta x\) 的取值，计算都是稳定的。这允许我们为了快速达到稳态解而使用较大的时间步长。缺点是每个时间步都需要解方程组，计算量比显式格式大。
Crank-Nicolson格式：平衡精度与稳定性
隐式欧拉格式虽然稳定，但它在时间上只有一阶精度，会引入较大的数值耗散（抹平效应）。为了兼顾稳定性和精度，Crank-Nicolson格式是一个优秀的选择。它将方程在时间层 \(n\) 和 \(n+1\) 之间取平均：
\[ \frac{U_j^{n+1} - U_j^n}{\Delta t} = \frac{a}{2} \left[ \frac{U_{j+1}^n - 2U_j^n + U_{j-1}^n}{(\Delta x)^2} + \frac{U_{j+1}^{n+1} - 2U_j^{n+1} + U_{j-1}^{n+1}}{(\Delta x)^2} \right] \]
整理后：
\[ -\frac{\mu}{2} U_{j-1}^{n+1} + (1 + \mu) U_j^{n+1} - \frac{\mu}{2} U_{j+1}^{n+1} = \frac{\mu}{2} U_{j-1}^{n} + (1 - \mu) U_j^{n} + \frac{\mu}{2} U_{j+1}^{n} \]
Crank-Nicolson格式也是隐式的，需要求解线性方程组。它同样是无条件稳定的，但更重要的是，它在时间和空间上都具有二阶精度，通常能给出比一阶格式精确得多的结果。
从初值问题到初边值问题
最后需要指出，纯粹的初值问题（无限区域）在数值上是通过取一个“足够大”的有限区域来近似实现的。一旦区域是有限的，我们就必须处理边界。此时，初值问题就演变为初边值问题。我们需要在区域边界 \(x = X_l\) 和 \(x = X_r\) 上施加物理上合理的边界条件（如狄利克雷条件、诺伊曼条件等），并将这些边界条件同样进行数值离散，与内部点的差分格式耦合起来进行求解。对边界条件的正确处理对于获得正确的数值解至关重要。