1.3.2差分格式的构造+硕博授课

一般椭圆型方程差分格式的构造方法详解

各位同学，上一讲我们建立了一般区域网格剖分、节点分类、控制体的基础框架，这一讲我们聚焦差分方法的核心环节——如何将连续的偏微分方程转化为离散的线性代数方程组。我们会从最基础的差分算子积木讲起，再介绍两种最主流的格式构造方法：直接差商替换法与守恒型有限体积法，最后明确差分格式构造的通用原则与核心要求。

---

一、构造差分格式的基础：差分算子

差分格式的本质是用离散网格节点的函数值的线性组合，逼近连续函数的各阶导数，而差分算子就是实现这一逼近的核心“积木”。我们以x方向的导数为例，所有算子均可直接推广到任意坐标方向。

1. 一阶差分算子

一阶差分算子用于逼近函数的一阶偏导数，根据所用节点的位置，分为向前、向后、中心三类，核心区别在于截断误差阶数与适用场景。

我们统一约定：\(\Delta x\)为x方向的步长，\(v(x,x')\)为多元函数，\(x'\)为其余不变的自变量，所有算子的截断误差均通过泰勒展开推导。

算子类型	算子定义	对应一阶差商（逼近\(\partial v/\partial x\)）	截断误差阶数	核心特点与适用场景
一阶向前差分算子\(\Delta_{+x}\)	\(\Delta_{+x}v(x,x') = v(x+\Delta x,x') - v(x,x')\)	\(\frac{\Delta_{+x}v}{\Delta x}\)	\(O(\Delta x)\)（一阶）	仅用到当前节点与右侧节点，适合左边界处的导数逼近；双曲型方程迎风格式的基础
一阶向后差分算子\(\Delta_{-x}\)	\(\Delta_{-x}v(x,x') = v(x,x') - v(x-\Delta x,x')\)	\(\frac{\Delta_{-x}v}{\Delta x}\)	\(O(\Delta x)\)（一阶）	仅用到当前节点与左侧节点，适合右边界处的导数逼近；与向前差分互补，处理边界条件
半节点一阶中心差分算子\(\delta_x\)	\(\delta_x v(x,x') = v(x+\frac{\Delta x}{2},x') - v(x-\frac{\Delta x}{2},x')\)	\(\frac{\delta_x v}{\Delta x}\)	\(O(\Delta x^2)\)（二阶）	用到半节点的函数值，二阶精度，是构造二阶导数差分算子的核心，适合变系数扩散项的离散
整节点一阶中心差分算子\(\Delta_{0x}\)	\(\Delta_{0x}v = \frac{1}{2}(\Delta_{+x}+\Delta_{-x})v = \frac{v(x+\Delta x)-v(x-\Delta x)}{2}\)	\(\frac{\Delta_{0x}v}{\Delta x} = \frac{v(x+\Delta x)-v(x-\Delta x)}{2\Delta x}\)	\(O(\Delta x^2)\)（二阶）	用到当前节点与左右两个整节点，无需半节点，二阶精度，适合对流项一阶导数的离散

关键推导：中心差分的二阶精度

以半节点中心差分为例，对\(v(x\pm\Delta x/2)\)做四阶泰勒展开：

\[v(x+\frac{\Delta x}{2}) = v(x) + \frac{\Delta x}{2}v'(x) + \frac{(\Delta x)^2}{8}v''(x) + \frac{(\Delta x)^3}{48}v'''(x) + O(\Delta x^4) \]

\[v(x-\frac{\Delta x}{2}) = v(x) - \frac{\Delta x}{2}v'(x) + \frac{(\Delta x)^2}{8}v''(x) - \frac{(\Delta x)^3}{48}v'''(x) + O(\Delta x^4) \]

两式相减后，偶次项完全抵消，得到：

\[\delta_x v = \Delta x \cdot v'(x) + O(\Delta x^3) \implies \frac{\delta_x v}{\Delta x} = v'(x) + O(\Delta x^2) \]

这就是中心差分二阶精度的来源——泰勒展开中的奇次误差项相互抵消，精度比一阶向前/向后差分高一阶。

2. 二阶中心差分算子

二阶差分算子用于逼近函数的二阶偏导数，是椭圆型方程离散的核心（椭圆型方程的主部为二阶导数项）。

我们通过连续两次应用半节点一阶中心差分算子，得到二阶中心差分算子：

\[\delta_x^2 v(x,x') = \delta_x\left(\delta_x v(x,x')\right) = v(x+\Delta x,x') - 2v(x,x') + v(x-\Delta x,x') \]

对应的二阶差商为\(\frac{\delta_x^2 v}{\Delta x^2}\)，用于逼近\(\frac{\partial^2 v}{\partial x^2}\)，截断误差为\(O(\Delta x^2)\)（二阶精度）。

核心验证：上一讲Poisson方程的五点格式，正是用这个二阶中心差分算子分别逼近x、y方向的二阶导数，得到了标准的五点差分模板。

---

二、差分格式的直接构造法：差商代替微商

这是差分格式最基础、最直观的构造方法，核心思想是：在正则内点上，将微分方程中出现的各阶偏导数，直接用对应精度的差分算子（差商）替换，一步得到离散的差分方程。

1. 适用场景

该方法适用于正则内点（差分模板完全在区域内）、系数光滑、解光滑的椭圆型方程，形式简单，易于推导和实现。

2. 实例：二维定常对流扩散方程的离散

我们以二维变系数定常对流扩散方程为例，完整演示直接构造法的过程。方程形式为：

\[-\nabla \cdot (a(x,y)\nabla u(x,y)) + \nabla \cdot (u(x,y)\boldsymbol{v}(x,y)) = f(x,y), \quad \forall(x,y)\in\Omega \tag{1.3.8} \]

其中\(\boldsymbol{v}=(v^1,v^2)\)为流体流速场，方程左侧第一项为变系数扩散项，第二项为对流项，分别对应分子扩散与流体对流的输运效应。

我们取步长\(h_x=\Delta x, h_y=\Delta y\)的矩形网格，节点\((i,j)\)对应坐标\((x_i,y_j)=(i\Delta x,j\Delta y)\)，对两项分别离散：

（1）变系数扩散项的离散

扩散项为散度形式：\(\nabla \cdot (a\nabla u) = \frac{\partial}{\partial x}\left(a\frac{\partial u}{\partial x}\right) + \frac{\partial}{\partial y}\left(a\frac{\partial u}{\partial y}\right)\)。
这里不能直接对\(u\)用二阶差分，因为扩散系数\(a(x,y)\)是空间函数，需要用半节点差分算子保证精度：

\[\left[\frac{\partial}{\partial x}\left(a\frac{\partial u}{\partial x}\right)\right]_{i,j} \approx \frac{\delta_x(a_{i,j}\delta_x u_{i,j})}{(\Delta x)^2} \]

将半节点差分展开，得到：

\[\delta_x(a_{i,j}\delta_x u_{i,j}) = a_{i+\frac{1}{2},j}(u_{i+1,j}-u_{i,j}) - a_{i-\frac{1}{2},j}(u_{i,j}-u_{i-1,j}) \]

其中半节点的扩散系数\(a_{i\pm1/2,j}\)通常取相邻节点的算术平均\(a_{i\pm1/2,j}=\frac{a_{i\pm1,j}+a_{i,j}}{2}\)，保证二阶精度。

y方向的扩散项同理，用\(\delta_y\)算子离散即可。

（2）对流项的离散

对流项为\(\nabla \cdot (u\boldsymbol{v}) = \frac{\partial(uv^1)}{\partial x} + \frac{\partial(uv^2)}{\partial y}\)，是一阶导数，用整节点中心差分算子\(\Delta_{0x}\)离散，保证二阶精度：

\[\left[\frac{\partial(uv^1)}{\partial x}\right]_{i,j} \approx \frac{\Delta_{0x}(uv^1)_{i,j}}{\Delta x} = \frac{(uv^1)_{i+1,j} - (uv^1)_{i-1,j}}{2\Delta x} \]

y方向的对流项同理，用\(\Delta_{0y}\)算子离散。

（3）最终差分格式

将离散后的扩散项、对流项代入原方程，整理得到直接法的差分格式：

\[\begin{aligned} &-\frac{a_{i+\frac{1}{2},j}(U_{i+1,j}-U_{i,j}) - a_{i-\frac{1}{2},j}(U_{i,j}-U_{i-1,j})}{(\Delta x)^2} \\ &-\frac{a_{i,j+\frac{1}{2}}(U_{i,j+1}-U_{i,j}) - a_{i,j-\frac{1}{2}}(U_{i,j}-U_{i,j-1})}{(\Delta y)^2} \\ &+\frac{(Uv^1)_{i+1,j} - (Uv^1)_{i-1,j}}{2\Delta x} + \frac{(Uv^2)_{i,j+1} - (Uv^2)_{i,j-1}}{2\Delta y} = f_{i,j} \end{aligned} \tag{1.3.9} \]

3. 直接构造法的优缺点

优点	缺点
推导简单直观，易于理解和编程实现	无法天然保证物理守恒性，离散格式不满足全局守恒律
对光滑系数、光滑解的问题，能达到设计的精度阶数	系数间断、解有梯度突变时，格式精度会下降，甚至出现非物理振荡
模板固定，适合规则矩形网格的正则内点	非正则内点、曲边边界处无法直接使用，需要额外处理

---

三、守恒型差分格式：有限体积法（积分插值法）

直接构造法的核心缺陷是不保证守恒性，而工程中的对流扩散、弹性力学、流体力学问题，均源于物理守恒律（质量、能量、动量守恒）。守恒型格式的核心思想是从方程的积分守恒形式出发，在控制体上离散，天然保持物理守恒性，也是目前工程数值计算的主流方法。

1. 守恒律的积分形式

原对流扩散方程的物理本质是控制体上的质量/能量守恒，对区域内任意子区域\(\omega\)，积分后利用散度定理，得到守恒型积分形式：

\[\oint_{\partial\omega} a\nabla u\cdot \boldsymbol{\nu} ds - \oint_{\partial\omega} u\boldsymbol{v}\cdot \boldsymbol{\nu} ds + \int_\omega f dxdy = 0 \tag{1.3.10} \]

其中\(\boldsymbol{\nu}\)是子区域边界\(\partial\omega\)的单位外法向量，第一项是扩散流出的总通量，第二项是对流流出的总通量，第三项是源项产生的总量，三者之和为0，完美对应“流入量+产生量=流出量”的守恒律。

2. 守恒型格式的构造步骤

我们以节点\((i,j)\)对应的控制体\(\omega_{i,j}\)为离散单元，完整演示构造过程：

步骤1：选取控制体

对节点\((i,j)\)，取以其为中心的矩形控制体：

\[\omega_{i,j} = \left( (i-\frac{1}{2})h_x, (i+\frac{1}{2})h_x \right) \times \left( (j-\frac{1}{2})h_y, (j+\frac{1}{2})h_y \right) \]

控制体有四条边界：右边界\(x=(i+1/2)h_x\)、左边界\(x=(i-1/2)h_x\)、上边界\(y=(j+1/2)h_y\)、下边界\(y=(j-1/2)h_y\)。

步骤2：对控制体积分守恒方程

将积分守恒形式(1.3.10)应用到控制体\(\omega_{i,j}\)上，把边界积分拆分为四条边的积分之和，体积分保留为源项积分。

步骤3：数值积分与差商逼近

对每条边的线积分，采用中点数值积分公式（二阶精度），边界上的法向导数、通量用相邻节点的差商与平均值逼近：

• 右边界积分：外法向为x正方向，扩散通量积分\(\int a\frac{\partial u}{\partial x} ds \approx h_y \cdot a_{i+1/2,j} \cdot \frac{u_{i+1,j}-u_{i,j}}{h_x}\)；对流通量积分\(\int u v^1 ds \approx h_y \cdot \frac{(u_{i+1,j}+u_{i,j})}{2} \cdot v^1_{i+1/2,j}\)。
• 左、上、下三条边界的积分，按外法向方向同理处理。
• 源项体积分：用中点积分近似为\(\int_{\omega_{i,j}} f dxdy \approx h_x h_y \cdot f_{i,j}\)。

步骤4：整理得到守恒型差分格式

将所有积分的近似结果代入守恒方程，两边除以\(h_x h_y\)，整理得到最终的守恒型差分格式：

\[\begin{aligned} &-\frac{a_{i+\frac{1}{2},j}(U_{i+1,j}-U_{i,j}) - a_{i-\frac{1}{2},j}(U_{i,j}-U_{i-1,j})}{(\Delta x)^2} \\ &-\frac{a_{i,j+\frac{1}{2}}(U_{i,j+1}-U_{i,j}) - a_{i,j-\frac{1}{2}}(U_{i,j}-U_{i,j-1})}{(\Delta y)^2} \\ &+\frac{(U_{i+1,j}+U_{i,j})v^1_{i+\frac{1}{2},j} - (U_{i,j}+U_{i-1,j})v^1_{i-\frac{1}{2},j}}{2\Delta x} \\ &+\frac{(U_{i,j+1}+U_{i,j})v^2_{i,j+\frac{1}{2}} - (U_{i,j}+U_{i,j-1})v^2_{i,j-\frac{1}{2}}}{2\Delta y} = f_{i,j} \end{aligned} \tag{1.3.12} \]

3. 守恒型格式的核心优势

1. 天然满足物理守恒律
   对整个区域的所有控制体求和时，内部相邻控制体的公共边界上，通量会大小相等、符号相反，完全抵消，最终只剩下区域外边界的总通量，与连续的全局守恒律完全一致，不会出现虚假的源汇，数值解能正确反映物理过程。

2. 完美处理间断系数问题
   当扩散系数\(a(x,y)\)在材料界面处间断时，守恒型格式能自动满足界面处的通量连续条件，无需额外的界面处理；而非守恒型格式在间断处会出现精度损失甚至错误，这是工程中多介质问题的核心需求。

3. 适用性极强
   不仅适用于规则矩形网格，还能直接推广到三角形、四面体等非结构化网格，适配复杂几何区域；同时对非线性问题、对流占优问题，守恒型格式的稳定性和收敛性远优于非守恒型格式。

方法定义：这种通过构造数值通量、在控制体上离散守恒律的方法，就是有限体积法；更一般地，通过方程的积分形式结合数值积分、差商逼近构造格式的方法，统称为积分插值法。

---

四、差分格式构造的通用方法与核心原则

1. 通用构造方法：待定系数法

上述两种方法都是针对特定方程的构造方式，而待定系数法是构造任意精度、任意模板差分格式的通用方法。

• 核心思想：假设用节点邻域内\(N\)个节点的函数值的线性组合，逼近某阶导数，通过泰勒展开令低阶误差项的系数为0，建立线性方程组，解出组合系数，得到满足指定精度的差分格式。
• 应用场景：构造高阶紧致格式、非均匀网格格式、非结构化网格格式，以及满足特殊要求（如最大值原理）的格式。

2. 椭圆型方程差分格式的核心构造原则

一个合格的椭圆型方程差分格式，必须满足以下核心要求：

1. 相容性：局部截断误差随步长\(h\to0\)而趋于0，保证离散格式能正确逼近原微分方程，是格式收敛的前提。
2. 守恒性：对守恒型方程，格式必须保持物理守恒律，尤其是间断系数、非线性问题，守恒性是收敛的必要条件。
3. 稳定性：右端项、边界条件的微小扰动，不会导致数值解的无限放大，保证离散方程组是良态的，是格式收敛的核心保障。
4. 离散最大值原理：椭圆型方程的解满足最大值原理，格式也必须满足离散最大值原理，避免出现非物理的数值振荡（如对流占优时，中心差分可能出现振荡，需改用迎风格式）。
5. 计算效率：格式的差分模板尽可能小，系数矩阵尽可能稀疏对称，降低计算与存储开销，适配大规模工程问题。

---

五、核心知识点归纳总结

表1 基础差分算子核心信息汇总

算子类型	算子表达式	逼近导数	截断误差	核心适用场景
一阶向前差分\(\Delta_{+x}\)	\(\Delta_{+x}v = v(x+\Delta x) - v(x)\)	一阶偏导	\(O(\Delta x)\)	左边界处理、双曲型方程迎风格式
一阶向后差分\(\Delta_{-x}\)	\(\Delta_{-x}v = v(x) - v(x-\Delta x)\)	一阶偏导	\(O(\Delta x)\)	右边界处理、与向前差分互补
半节点中心差分\(\delta_x\)	\(\delta_x v = v(x+\Delta x/2) - v(x-\Delta x/2)\)	一阶偏导	\(O(\Delta x^2)\)	变系数扩散项离散、二阶差分算子构造
整节点中心差分\(\Delta_{0x}\)	\(\Delta_{0x}v = \frac{v(x+\Delta x)-v(x-\Delta x)}{2}\)	一阶偏导	\(O(\Delta x^2)\)	对流项一阶导数离散、光滑区域内部导数逼近
二阶中心差分\(\delta_x^2\)	\(\delta_x^2 v = v(x+\Delta x) - 2v(x) + v(x-\Delta x)\)	二阶偏导	\(O(\Delta x^2)\)	椭圆型方程主部离散、Poisson方程五点格式核心

表2 两种差分格式构造方法对比

对比维度	直接构造法（差商代替微商）	守恒型格式（有限体积法）
构造出发点	微分方程的强形式（点态）	微分方程的积分守恒形式
物理守恒性	不天然保证	天然满足全局守恒律
间断系数适配性	差、精度损失	好、自动满足界面通量连续
网格适配性	仅适合规则矩形网格	适配规则/非结构化网格
推导与实现难度	简单、易上手	稍复杂、工程通用性强
核心适用场景	光滑系数、光滑解的规则区域问题	工程实际问题、多介质、非线性、复杂区域问题