昆仑山:眼中无形心中有穴之穴人合一

夫君子之行,静以修身,俭以养德;非澹泊无以明志,非宁静无以致远。夫学须静也,才须学也;非学无以广才,非志无以成学。怠慢则不能励精,险躁则不能冶性。年与时驰,意与岁去,遂成枯落,多不接世。悲守穷庐,将复何及!

 

ch02+椭圆型方程的有限差分法

一、数值解法的核心本质:连续问题离散化

1. 离散化的根本原因

计算机只能存储有限个数据并执行有限次运算,而偏微分方程的定解问题(边值问题、初值问题)定义在连续的求解域上,包含无穷多个未知点,无法直接求解。因此,所有数值方法的核心都是将连续问题转化为有限维的代数问题,最终通过求解线性代数方程组得到近似解。

2. 通用离散化两步法

所有主流数值方法(有限差分法、有限元法、有限体积法等)都遵循以下两个核心步骤:

  1. 网格剖分:将连续的求解域划分为有限个互不重叠的子区域(单元),用有限个网格节点代替连续区域。
  2. 算子离散化:将微分方程中的微分算子(如一阶导数、二阶导数、拉普拉斯算子等)转化为离散节点上的代数运算,从而将微分方程转化为线性代数方程组。

二、有限差分法与Galerkin有限元法的核心区别

两种方法的本质差异体现在算子离散化的出发点和数学基础上,这也是它们适用场景不同的根本原因。

方法 出发点 数学基础 核心思想 基函数特点
有限差分法(FDM) 微分方程的微分形式积分形式 泰勒展开、数值微商、数值积分 差商代替微商,直接用节点函数值的线性组合近似导数 无显式基函数,直接基于节点值构造差分格式
Galerkin有限元法(FEM) 微分方程的变分形式(弱形式) 变分原理、Ritz-Galerkin方法 将解表示为有限维基函数的线性组合,代入变分方程得到代数方程组 分片多项式基函数,仅在相邻单元上非零,系数矩阵稀疏

关键补充

  • 有限差分法更直观,编程简单,适合规则区域和光滑解问题;
  • 有限元法网格灵活性极强,能高效处理复杂几何区域和复杂边界条件,是目前工程计算中应用最广泛的数值方法。

三、有限差分法的四大基本问题深度解析

问题1:对求解域作网格剖分

网格剖分是差分法的第一步,其质量直接影响数值解的精度和计算效率。

(1) 一维网格剖分

求解域为区间 \([a,b]\),将其划分为 \(N\) 个小区间(单元),节点定义为:

\[x_i = a + i \cdot h, \quad i=0,1,\dots,N \]

其中 \(h = \frac{b-a}{N}\)等距步长;若为不等距网格,则步长 \(h_i = x_{i+1} - x_i\) 随节点变化。

(2) 二维网格剖分

  • 矩形网格(最常用):求解域为矩形 \(\Omega = [a,b] \times [c,d]\),x方向步长 \(h\),y方向步长 \(k\),节点为:

    \[(x_i, y_j) = (a+ih, c+jk), \quad i=0,\dots,M; \ j=0,\dots,N \]

  • 非矩形网格:可剖分为三角形或凸四边形网格,但有限差分法处理非矩形网格时,差分格式的构造会变得非常复杂,这是其主要局限性之一。

问题2:构造逼近微分方程的差分格式(核心)

构造差分格式的本质是用节点函数值的线性组合近似微分算子。本文介绍两种最基础的方法,其中直接差分化法(泰勒展开法)是所有差分格式的基础。

方法1:直接差分化法(泰勒展开法)

证明/推导依据:泰勒中值定理
若函数 \(u(x)\) 在包含 \(x_i\) 的开区间内具有直到 \(n+1\) 阶导数,则对任意 \(x\) 有:

\[u(x) = u(x_i) + u'(x_i)(x-x_i) + \frac{u''(x_i)}{2!}(x-x_i)^2 + \dots + \frac{u^{(n)}(x_i)}{n!}(x-x_i)^n + R_n(x) \]

其中余项 \(R_n(x) = \frac{u^{(n+1)}(\xi)}{(n+1)!}(x-x_i)^{n+1}\)\(\xi\) 介于 \(x_i\)\(x\) 之间。

一阶导数的差分近似

\(u(x)\) 为光滑函数,在节点 \(x_i\) 处对 \(x_{i+1}=x_i+h\)\(x_{i-1}=x_i-h\) 作泰勒展开:

\[u(x_{i+1}) = u(x_i) + h u'(x_i) + \frac{h^2}{2} u''(x_i) + \frac{h^3}{6} u'''(x_i) + O(h^4) \tag{1} \]

\[u(x_{i-1}) = u(x_i) - h u'(x_i) + \frac{h^2}{2} u''(x_i) - \frac{h^3}{6} u'''(x_i) + O(h^4) \tag{2} \]

  • 向前差分(由式(1)推导):
    移项得 \(u'(x_i) = \frac{u(x_{i+1})-u(x_i)}{h} - \frac{h}{2}u''(x_i) + O(h^2)\),舍去高阶小项得:

    \[u'(x_i) \approx \frac{u_{i+1} - u_i}{h} \]

    截断误差:\(O(h)\)(一阶精度)。

  • 向后差分(由式(2)推导):
    移项得 \(u'(x_i) = \frac{u(x_i)-u(x_{i-1})}{h} + \frac{h}{2}u''(x_i) + O(h^2)\),舍去高阶小项得:

    \[u'(x_i) \approx \frac{u_i - u_{i-1}}{h} \]

    截断误差:\(O(h)\)(一阶精度)。

  • 中心差分(由式(1)-式(2)推导):
    两式相减得 \(u(x_{i+1}) - u(x_{i-1}) = 2h u'(x_i) + \frac{h^3}{3}u'''(x_i) + O(h^5)\),移项得:

    \[u'(x_i) = \frac{u(x_{i+1})-u(x_{i-1})}{2h} - \frac{h^2}{6}u'''(x_i) + O(h^4) \]

    舍去高阶小项得:

    \[u'(x_i) \approx \frac{u_{i+1} - u_{i-1}}{2h} \]

    截断误差:\(O(h^2)\)(二阶精度)。

二阶导数的差分近似

将式(1)和式(2)相加:

\[u(x_{i+1}) + u(x_{i-1}) = 2u(x_i) + h^2 u''(x_i) + \frac{h^4}{12}u''''(x_i) + O(h^6) \]

移项得:

\[u''(x_i) = \frac{u(x_{i+1}) - 2u(x_i) + u(x_{i-1})}{h^2} - \frac{h^2}{12}u''''(x_i) + O(h^4) \]

舍去高阶小项得二阶中心差分

\[u''(x_i) \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} \]

截断误差:\(O(h^2)\)(二阶精度)。

实例:构造两点边值问题的差分格式

考虑二阶常微分方程边值问题:

\[\begin{cases} -u''(x) = f(x), & x \in (0,1) \\ u(0) = \alpha, & u(1) = \beta \end{cases}\]

用二阶中心差分近似 \(u''(x_i)\),代入方程得:

\[-\frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} = f_i, \quad i=1,2,\dots,N-1 \]

其中 \(f_i = f(x_i)\),边界条件为 \(u_0 = \alpha\)\(u_N = \beta\)

整理为线性代数方程组:

\[\begin{cases} 2u_1 - u_2 = h^2 f_1 + \alpha \\ -u_{i-1} + 2u_i - u_{i+1} = h^2 f_i, & i=2,\dots,N-2 \\ -u_{N-2} + 2u_{N-1} = h^2 f_{N-1} + \beta \end{cases}\]

方法2:有限体积法(积分插值法)

核心思想:从守恒律的积分形式出发,将控制方程在每个控制体积上积分,利用牛顿-莱布尼茨公式将体积分转化为面积分(通量),再用数值积分近似通量项。

有限体积法的最大优点是天然保持守恒性,因此在计算流体力学、电磁学等守恒律问题中应用极为广泛。

问题3:差分解的存在唯一性、收敛性及稳定性

这是差分法的理论基础,确保我们构造的差分格式是“数学上可靠”的。

(1) 存在唯一性

证明依据:对称正定矩阵的性质
若线性代数方程组 \(Ax=b\) 的系数矩阵 \(A\) 是对称正定矩阵,则方程组存在唯一解。

以上述两点边值问题的差分格式为例,其系数矩阵为三对角矩阵:

\[A = \begin{pmatrix} 2 & -1 & 0 & \dots & 0 \\ -1 & 2 & -1 & \dots & 0 \\ 0 & -1 & 2 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & 2 \end{pmatrix}\]

证明

  1. 矩阵 \(A\) 显然是对称的;
  2. 所有顺序主子式均大于0(一阶主子式 \(2>0\),二阶主子式 \(2\times2 - (-1)\times(-1)=3>0\),三阶主子式 \(2\times3 - (-1)\times(-2)=4>0\),以此类推,\(n\) 阶主子式为 \(n+1>0\));
  3. 因此 \(A\) 是对称正定矩阵,方程组存在唯一解。

(2) 收敛性

定义:当步长 \(h \to 0\) 时,差分解 \(u_i\) 收敛到精确解 \(u(x_i)\),即 \(\max_{0\leq i\leq N} |u(x_i) - u_i| \to 0\)

核心定理:Lax等价定理(线性初值问题)
依据:对于线性偏微分方程的初值问题,相容且稳定的差分格式是收敛的。

  • 相容性:差分格式的截断误差当 \(h \to 0\) 时趋于0;
  • 稳定性:计算过程中的初始误差和舍入误差不会无限放大。

对于上述两点边值问题,可通过误差分析证明其差分解是二阶收敛的,即 \(\max_{0\leq i\leq N} |u(x_i) - u_i| \leq C h^2\),其中 \(C\) 是与 \(h\) 无关的常数。

(3) 稳定性

定义:在数值计算中,若初始误差或舍入误差随着计算步数的增加而有界,则称差分格式是稳定的。

分析方法:冯·诺依曼稳定性分析
通过将误差分解为傅里叶级数,分析每个傅里叶分量的增长因子。若所有分量的增长因子的模都不超过1,则格式是稳定的。

典型例子:热传导方程的显式差分格式稳定性条件为 \(r = \frac{a^2 \tau}{h^2} \leq \frac{1}{2}\)\(\tau\) 为时间步长),而隐式Crank-Nicolson格式是无条件稳定的。

问题4:差分方程的解法

差分格式最终转化为线性代数方程组 \(Ax=b\),根据系数矩阵的特点选择合适的解法:

解法类型 具体方法 适用场景 时间复杂度
直接解法 追赶法(托马斯算法) 一维问题的三对角方程组 \(O(N)\)
直接解法 带状LU分解 二维问题的带状矩阵 \(O(N^2)\)
迭代解法 雅可比迭代、高斯-赛德尔迭代 小型稀疏矩阵 收敛较慢
迭代解法 共轭梯度法(CG)、预处理共轭梯度法(PCG) 大型对称正定稀疏矩阵 收敛较快
迭代解法 多重网格法 超大型椭圆型方程组 收敛速度与网格规模无关

四、核心知识点归纳总结表

核心知识点 关键内容 推导/证明依据 重要结论
数值解法本质 连续问题离散化,转化为线性代数方程组 计算机的有限性原理 所有数值方法都遵循“网格剖分+算子离散化”两步法
FDM与FEM区别 FDM从微分形式出发,用差商代替微商;FEM从变分形式出发,用基函数逼近解 泰勒展开、变分原理 FDM适合规则区域,FEM适合复杂几何区域
网格剖分 一维等距/不等距网格;二维矩形/非矩形网格 区域分解原理 网格质量直接影响数值解的精度和效率
差分格式构造 直接差分化法(泰勒展开)、有限体积法 泰勒中值定理、守恒律积分形式 二阶中心差分具有二阶精度,是最常用的差分近似
存在唯一性 差分方程组解的存在性和唯一性 对称正定矩阵的性质 两点边值问题的差分格式系数矩阵对称正定,存在唯一解
收敛性与稳定性 收敛性:差分解→精确解;稳定性:误差不放大 Lax等价定理、冯·诺依曼分析 相容且稳定的格式是收敛的;显式格式通常有稳定性条件
差分方程解法 直接解法(追赶法)、迭代解法(多重网格法) 线性代数方程组求解理论 一维问题用追赶法,三维问题用多重网格法效率最高

后续内容预告

  • 第2-4章:系统讲解有限差分法,包括椭圆型、抛物型、双曲型偏微分方程的差分格式构造、理论分析和实现方法;
  • 第5-6章:系统讲解Galerkin有限元法,包括变分原理、有限元空间构造、单元分析、整体合成和误差分析。

2.1 差分逼近的基本概念 深度讲解与完整推导

一、模型问题:二阶常微分方程第一边值问题

1. 问题描述

我们研究如下最简单的椭圆型方程第一边值问题:

\[\begin{cases} Lu = -\frac{d^2u}{dx^2} + qu = f, & a < x < b \\ u(a) = \alpha, \quad u(b) = \beta \end{cases} \]

其中 \(q,f \in C[a,b]\)(连续函数)且 \(q \geq 0\)\(\alpha,\beta\) 为给定常数。

重要性:该问题是所有偏微分方程数值解法的"原型问题",其差分格式的构造思想、理论分析方法可直接推广到二维、三维椭圆型方程,以及抛物型、双曲型方程。

2. 一维等距网格剖分

将求解区间 \([a,b]\) 划分为 \(N\) 个等长的小区间,定义:

  • 节点:\(x_i = a + ih, \quad i=0,1,\dots,N\)
  • 步长:\(h = \frac{b-a}{N}\)

由此得到区间 \(I=[a,b]\) 的一个均匀网格剖分,\(x_0=a\)\(x_N=b\) 称为边界节点\(x_1,x_2,\dots,x_{N-1}\) 称为内部节点

二、中心差分格式的构造与截断误差(核心推导)

1. 二阶导数的中心差分近似

推导依据:泰勒中值定理
若函数 \(u(x)\)\(x_i\) 的邻域内具有四阶连续导数,则对 \(x_{i+1}=x_i+h\)\(x_{i-1}=x_i-h\) 分别作泰勒展开:

\[u(x_{i+1}) = u(x_i) + h u'(x_i) + \frac{h^2}{2!} u''(x_i) + \frac{h^3}{3!} u'''(x_i) + \frac{h^4}{4!} u''''(x_i) + O(h^5) \tag{A} \]

\[u(x_{i-1}) = u(x_i) - h u'(x_i) + \frac{h^2}{2!} u''(x_i) - \frac{h^3}{3!} u'''(x_i) + \frac{h^4}{4!} u''''(x_i) + O(h^5) \tag{B} \]

将式(A)和式(B)相加,消去一阶导数和三阶导数项:

\[u(x_{i+1}) + u(x_{i-1}) = 2u(x_i) + h^2 u''(x_i) + \frac{h^4}{12} u''''(x_i) + O(h^5) \]

移项整理得到二阶导数的表达式:

\[u''(x_i) = \frac{u(x_{i+1}) - 2u(x_i) + u(x_{i-1})}{h^2} - \frac{h^2}{12} u''''(x_i) + O(h^3) \tag{2.3} \]

这就是教材中式(2.3)的完整推导过程。

2. 差分方程的建立

将式(2.3)代入原微分方程(2.1),得到:

\[-\left[ \frac{u(x_{i+1}) - 2u(x_i) + u(x_{i-1})}{h^2} - \frac{h^2}{12} u''''(x_i) + O(h^3) \right] + q(x_i)u(x_i) = f(x_i) \]

整理后:

\[-\frac{u(x_{i+1}) - 2u(x_i) + u(x_{i-1})}{h^2} + q(x_i)u(x_i) = f(x_i) - \frac{h^2}{12} u''''(x_i) + O(h^3) \tag{2.4} \]

3. 截断误差的定义与分析

定义截断误差 \(R_i(u)\) 为:

\[R_i(u) = -\frac{h^2}{12} \left[ \frac{d^4 u(x)}{dx^4} \right]_i + O(h^3) \tag{2.5} \]

关键性质:当 \(h \to 0\) 时,\(R_i(u) = O(h^2)\),即截断误差是步长 \(h\) 的二阶无穷小量。

截断误差的算子意义
定义差分算子 \(L_h\) 为:

\[L_h u_i = -\frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} + q_i u_i \]

其中 \(q_i = q(x_i)\)\(u_i\)\(u(x_i)\) 的近似值。

则式(2.4)可写为:

\[L_h u(x_i) = f(x_i) + R_i(u) \tag{2.7} \]

而原微分方程在节点 \(x_i\) 处为:

\[[Lu]_i = f(x_i) \]

两式相减得到截断误差的算子表达式:

\[R_i(u) = L_h u(x_i) - [Lu]_i \tag{2.8} \]

核心结论:截断误差是用差分算子 \(L_h\) 逼近微分算子 \(L\) 所产生的误差,其阶数 \(O(h^2)\) 决定了差分格式的精度阶数。

4. 中心差分格式

舍去截断误差项 \(R_i(u)\),得到逼近原微分方程的差分方程:

\[L_h u_i = -\frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} + q_i u_i = f_i, \quad i=1,2,\dots,N-1 \tag{2.6} \]

其中 \(f_i = f(x_i)\)

结合边界条件 \(u_0 = \alpha\)\(u_N = \beta\),得到完整的差分方程组:

\[\begin{cases} -\frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} + q_i u_i = f_i, & i=1,2,\dots,N-1 \\ u_0 = \alpha, \quad u_N = \beta \end{cases} \]

该格式称为中心差分格式,是二阶精度的差分格式。

三、差分方程组的系数矩阵结构

1. 线性方程组的矩阵形式

将差分方程组(2.9)-(2.10)整理为线性代数方程组 \(A\mathbf{u} = \mathbf{b}\),其中未知向量 \(\mathbf{u} = (u_1, u_2, \dots, u_{N-1})^T\)

\(N=5\)(即4个内部节点)为例,系数矩阵 \(A\) 为:

\[A = \begin{pmatrix} \frac{2}{h^2} + q_1 & -\frac{1}{h^2} & 0 & 0 \\ -\frac{1}{h^2} & \frac{2}{h^2} + q_2 & -\frac{1}{h^2} & 0 \\ 0 & -\frac{1}{h^2} & \frac{2}{h^2} + q_3 & -\frac{1}{h^2} \\ 0 & 0 & -\frac{1}{h^2} & \frac{2}{h^2} + q_4 \end{pmatrix} \]

2. 系数矩阵的关键性质

  • 三对角性:仅在主对角线和相邻的两条次对角线上有非零元素,其余元素均为0;
  • 稀疏性:非零元素个数约为 \(3(N-1)\),远小于总元素数 \((N-1)^2\),当 \(N\) 很大时稀疏性尤为显著;
  • 对称性\(A^T = A\),因为 \(A_{ij} = A_{ji}\)
  • 严格对角占优性:对任意 \(i\),有 \(|A_{ii}| > \sum_{j \neq i} |A_{ij}|\)
    证明:对角元 \(A_{ii} = \frac{2}{h^2} + q_i \geq \frac{2}{h^2}\)(因 \(q_i \geq 0\)),而非对角元绝对值之和为 \(\frac{1}{h^2}\),显然 \(\frac{2}{h^2} > \frac{1}{h^2}\),故严格对角占优。

3. 解的存在唯一性定理

定理:差分方程组(2.9)-(2.10)存在唯一解。
证明依据:严格对角占优矩阵的性质
证明:系数矩阵 \(A\) 是对称严格对角占优矩阵,根据线性代数理论,对称严格对角占优矩阵必为对称正定矩阵,而对称正定矩阵对应的线性方程组存在唯一解。

四、网格函数的范数体系

为了定量描述差分解的误差以及分析收敛性和稳定性,我们需要在网格函数空间上定义范数。

\(I_h\) 表示内部节点集合 \(\{x_1,x_2,\dots,x_{N-1}\}\)\(\bar{I}_h\) 表示全体节点集合 \(\{x_0,x_1,\dots,x_N\}\)。定义在 \(I_h\) 上的函数 \(u_h(x_i) = u_i\) 称为网格函数

1. 常用范数定义

范数名称 定义式 物理意义
最大范数(C范数) \(|u_h|_C = \max_{1 \leq i \leq N-1} |u_i|\) 衡量网格函数的逐点最大误差
离散L²范数 \(|u_h|_0^2 = \sum_{i=1}^{N-1} h u_i^2\) 衡量网格函数的整体平均平方误差
离散H¹半范数 \(|u_h|_1^2 = \sum_{i=1}^N h \left( \frac{u_i - u_{i-1}}{h} \right)^2\) 衡量网格函数一阶导数的平均平方误差
离散H¹范数(能量范数) \(|u_h|_1^2 = |u_h|_0^2 + |u_h|_1^2\) 同时衡量函数值和导数的误差,是椭圆型方程误差分析中最常用的范数

2. 范数的意义

  • 不同的范数对应不同的误差度量方式;
  • 收敛性的定义依赖于范数的选择,例如:
    • \(\|u_h - u\|_C \to 0\)\(h \to 0\)),则称差分解一致收敛到精确解;
    • \(\|u_h - u\|_1 \to 0\)\(h \to 0\)),则称差分解按能量范数收敛

五、差分法的核心理论问题

差分格式建立后,我们必须回答以下三个基本理论问题:

  1. 存在唯一性:差分方程组是否存在唯一解?(本节已证明答案是肯定的)
  2. 收敛性:当网格无限加密(\(h \to 0\))时,差分解 \(u_i\) 是否收敛到精确解 \(u(x_i)\)?收敛的速度是多少?
  3. 稳定性:在数值计算过程中,初始误差和舍入误差是否会被无限放大?

这三个问题是差分法数学理论的核心,只有这三个问题都得到肯定的回答,我们构造的差分格式才是可靠和可用的。

六、核心知识点归纳总结表

核心知识点 关键内容 推导/证明依据 重要结论
模型问题 二阶常微分方程第一边值问题 \(-u''+qu=f\)\(q\geq0\) 椭圆型方程基本理论 是所有偏微分方程数值解法的原型问题
网格剖分 一维等距网格,节点 \(x_i=a+ih\),步长 \(h=(b-a)/N\) 区域分解原理 均匀网格是最简单、最常用的网格形式
中心差分格式 \(-\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}+q_i u_i=f_i\) 泰勒中值定理 二阶精度的差分格式,形式简单对称
截断误差 \(R_i(u)=L_h u(x_i)-[Lu]_i=O(h^2)\) 泰勒展开余项分析 截断误差的阶数决定了差分格式的精度
系数矩阵 对称三对角严格对角占优矩阵 差分方程的结构 稀疏性好,计算效率高
解的存在唯一性 差分方程组存在唯一解 严格对角占优矩阵必正定 保证了差分格式的适定性
网格函数范数 最大范数、L²范数、H¹范数 泛函分析中的范数理论 是分析收敛性和稳定性的数学工具
核心理论问题 存在唯一性、收敛性、稳定性 数值分析基本理论 是差分格式可靠性的理论基础

差分逼近的相容性、收敛性与稳定性 深度讲解

一、相容性:差分算子对微分算子的逼近

1. 定义2.1 相容性

定义:设 \(M\) 是某一充分光滑的函数类,\(R_h(u)\) 是由截断误差(2.8)定义的网格函数。若对任何 \(u \in M\),恒有

\[\lim_{h \to 0} \|R_h(u)\| = 0 \tag{2.15} \]

则称差分算子 \(L_h\) 逼近微分算子 \(L\),而称(2.15)为相容条件

本质解读:相容性描述的是算子层面的逼近性质。当网格步长 \(h\) 趋于0时,差分算子作用在精确解上的结果,应该无限接近微分算子作用在精确解上的结果。换句话说,差分方程在极限情况下应该还原为原微分方程。

2. 中心差分格式的相容性验证

对于我们之前构造的中心差分格式,其截断误差为:

\[R_i(u) = -\frac{h^2}{12} \left[ \frac{d^4 u(x)}{dx^4} \right]_i + O(h^3) \]

我们可以计算其在不同范数下的阶数:

  • 最大范数\(\|R_h(u)\|_C = \max_{1 \leq i \leq N-1} |R_i(u)| = O(h^2)\)
  • 离散L²范数\(\|R_h(u)\|_0^2 = \sum_{i=1}^{N-1} h R_i(u)^2 = \sum_{i=1}^{N-1} h \cdot O(h^4) = O(h^4)\),故 \(\|R_h(u)\|_0 = O(h^2)\)
  • 离散H¹范数\(|R_h(u)|_1^2 = \sum_{i=1}^N h \left( \frac{R_i(u) - R_{i-1}(u)}{h} \right)^2 = \sum_{i=1}^N h \cdot \frac{O(h^4)}{h^2} = O(h^3)\),故 \(\|R_h(u)\|_1 = O(h)\)

结论:中心差分格式在最大范数和L²范数下是二阶相容的,在H¹范数下是一阶相容的。

二、收敛性:差分解对精确解的逼近

1. 定义2.2 收敛性

定义:称差分解 \(u_h\) 收敛到边值问题的解 \(u\),如果当 \(h\) 充分小时,差分方程组(2.9)-(2.10)的解 \(u_h\) 存在,且按某一范数 \(\|\cdot\|\)

\[\lim_{h \to 0} \|u_h - u\| = 0 \tag{2.16} \]

这里把 \(u\) 看成 \(I_h\) 上的网格函数。

与相容性的本质区别

  • 相容性:差分算子 \(L_h\) 逼近微分算子 \(L\)(算子层面)
  • 收敛性:差分解 \(u_h\) 逼近精确解 \(u\)(解层面)

相容性是收敛性的必要条件,但不是充分条件。一个相容的差分格式可能是发散的,只有同时满足稳定性,才能保证收敛性。

2. 误差方程的推导(连接相容性与收敛性的核心桥梁)

推导依据:差分方程的精确解形式与差分解形式

原微分方程的精确解 \(u(x)\) 满足:

\[L_h u(x_i) = f_i + R_i(u) \tag{2.4} \]

差分解 \(u_i\) 满足差分方程:

\[L_h u_i = f_i \tag{2.9} \]

将两式相减,得到:

\[L_h (u(x_i) - u_i) = R_i(u) \]

引进误差函数 \(e_i = u(x_i) - u_i\),则误差函数满足如下差分方程:

\[\begin{cases} L_h e_i = R_i(u), & i=1,2,\dots,N-1 \\ e_0 = e_N = 0 \end{cases} \tag{2.17} \]

重大意义:收敛性问题(估计 \(\|e_h\|\))被完全转化为了一个纯代数问题——通过误差方程的右端项 \(R_i(u)\)(截断误差)来估计误差函数 \(e_h\) 的范数。而这个估计问题,本质上就是差分方程的稳定性问题。

三、稳定性:差分方程解对右端项的连续依赖性

1. 定义2.3 关于右端的稳定性

定义:称差分方程 \(L_h v_i = f_i \ (i=1,2,\dots,N-1), v_0 = v_N = 0\) 关于右端稳定,如果存在与网格 \(I_h\) 及右端 \(f_h\)\(f_h(x_i)=f_i\))无关的正常数 \(M\)\(h_0\),使

\[\|v_h\| \leq M \|f_h\|_R, \quad 0 < h < h_0 \tag{2.18} \]

其中 \(\|f_h\|_R\) 是右端 \(f_h\) 的某一范数,它可以和 \(\|\cdot\|\) 相同,也可以不同。

本质解读:稳定性描述的是差分方程本身的良定性。它意味着:右端项的微小变化,只会引起解的微小变化。换句话说,计算过程中产生的舍入误差和初始误差,不会随着网格的加密而被无限放大。

2. 稳定性的两个重要推论

推论1:齐次差分方程只有平凡解
证明依据:稳定性定义
\(f_i = 0\)(齐次方程),则由稳定性不等式(2.18)得:

\[\|v_h\| \leq M \|0\|_R = 0 \]

\(v_h = 0\),即齐次方程只有零解。

推论2:非齐次差分方程对任何右端项和边界条件都有唯一解
证明依据:线性方程组解的存在唯一性定理
对于线性方程组 \(A\mathbf{v} = \mathbf{f}\),齐次方程只有零解等价于系数矩阵 \(A\) 可逆,而可逆矩阵对应的非齐次方程组对任何右端项都有唯一解。

:这与我们之前通过"严格对角占优矩阵必正定"证明的解的存在唯一性是完全一致的,从稳定性的角度给出了另一种解释。

四、核心定理:相容 + 稳定 = 收敛

这是数值分析中最基本、最重要的定理之一,它将差分法的三个核心概念完美地联系在了一起。

1. 定理内容

若差分算子 \(L_h\) 关于范数 \(\|\cdot\|_R\) 满足相容条件,且差分方程关于范数 \(\|\cdot\|\)\(\|\cdot\|_R\) 是稳定的,则差分解 \(u_h\) 按范数 \(\|\cdot\|\) 收敛到精确解 \(u\),且收敛阶与截断误差的阶相同。

2. 严格证明

证明依据:误差方程与稳定性定义

  1. 误差函数 \(e_h\) 满足误差方程(2.17):\(L_h e_h = R_h(u)\),且 \(e_0 = e_N = 0\)
  2. 根据稳定性定义(2.18),存在与 \(h\) 无关的常数 \(M\),使得:

    \[\|e_h\| \leq M \|R_h(u)\|_R \tag{2.19} \]

  3. 根据相容性条件(2.15),当 \(h \to 0\) 时,\(\|R_h(u)\|_R \to 0\)
  4. 因此,当 \(h \to 0\) 时,\(\|e_h\| \to 0\),即差分解收敛到精确解。
  5. 进一步,若截断误差的阶为 \(O(h^p)\),即 \(\|R_h(u)\|_R = O(h^p)\),则由(2.19)得:

    \[\|e_h\| = O(h^p) \]

    即收敛阶为 \(p\) 阶。

3. 定理的深远意义

  • 它为判断差分格式的可靠性提供了一个可操作的标准:我们不需要直接证明收敛性(这通常非常困难),只需要分别验证相容性稳定性即可。
  • 它明确了收敛阶的决定因素:收敛阶等于截断误差的阶。这为我们提高数值解的精度指明了方向——构造更高阶截断误差的差分格式。
  • 它是后续所有偏微分方程数值解法(抛物型、双曲型)理论分析的基石,著名的Lax等价定理就是这个定理在初值问题上的推广。

五、核心知识点归纳总结表

核心概念 定义 数学表达式 本质 意义 中心差分格式的阶数
相容性 差分算子逼近微分算子 \(\lim_{h \to 0} |R_h(u)| = 0\) 算子层面的逼近 收敛性的必要条件 \(|R_h|_C=O(h^2), |R_h|_0=O(h^2), |R_h|_1=O(h)\)
收敛性 差分解逼近精确解 \(\lim_{h \to 0} |u_h - u| = 0\) 解层面的逼近 数值方法的最终目标 二阶收敛(最大范数、L²范数)
稳定性 解连续依赖于右端项 \(|v_h| \leq M |f_h|_R\) 差分格式的良定性 收敛性的充分条件 无条件稳定
核心关系 相容 + 稳定 = 收敛 \(|e_h| \leq M |R_h(u)|\) 误差由截断误差控制 数值分析的基本定理 收敛阶等于截断误差阶

定理2.1(收敛性基本定理)深度解析与完整证明

一、定理2.1 完整表述与核心地位

1. 定理原文

定理2.1 若边值问题的解 \(u\) 充分光滑,差分方程按 \(\|\cdot\|_R\) 满足相容条件,且关于右端稳定,则差分解 \(u_h\)\(\|\cdot\|\) 收敛到边值问题的解,且有和 \(\|R_h(u)\|_R\) 相同的收敛阶。

2. 核心地位

这是椭圆型方程差分法理论的基石性定理,它彻底解决了"如何判断一个差分格式是否收敛"的根本问题,将原本复杂的收敛性证明分解为两个相对独立且可操作的步骤:检验相容性证明稳定性

二、定理2.1 严格证明

证明依据:误差方程 + 稳定性定义 + 相容性定义

  1. 误差方程的建立
    \(e_i = u(x_i) - u_i\) 为节点 \(x_i\) 处的误差,则误差函数 \(e_h\) 满足:

    \[\begin{cases} L_h e_i = R_i(u), & i=1,2,\dots,N-1 \\ e_0 = e_N = 0 \end{cases} \]

    其中 \(R_i(u)\) 是截断误差。该方程的推导过程已在之前的内容中详细给出,是连接相容性与收敛性的核心桥梁。

  2. 应用稳定性条件
    根据关于右端稳定的定义(定义2.3),存在与网格 \(I_h\) 和右端项无关的正常数 \(M\)\(h_0\),使得当 \(0 < h < h_0\) 时,有:

    \[\|e_h\| \leq M \|R_h(u)\|_R \tag{*} \]

    这一步是整个证明的关键,它将误差的范数用截断误差的范数进行了控制。

  3. 应用相容性条件
    根据相容性条件(定义2.1),当 \(h \to 0\) 时,有:

    \[\|R_h(u)\|_R \to 0 \]

  4. 收敛性结论
    结合式(*)和相容性条件,当 \(h \to 0\) 时:

    \[\|e_h\| \leq M \|R_h(u)\|_R \to 0 \]

    \(\lim_{h \to 0} \|u_h - u\| = 0\),差分解按范数 \(\|\cdot\|\) 收敛到精确解。

  5. 收敛阶结论
    若截断误差的阶为 \(O(h^p)\),即 \(\|R_h(u)\|_R = O(h^p)\),则由式(*)直接可得:

    \[\|e_h\| = O(h^p) \]

    即差分解的收敛阶与截断误差的阶完全相同。

证毕

三、定理核心要点深度解读

1. 三个前提条件的作用

前提条件 作用 验证方法
\(u\) 充分光滑 保证截断误差可以通过泰勒展开进行估计,且高阶导数有界 由原微分方程的正则性理论保证
差分方程满足相容条件 保证差分算子在极限情况下还原为微分算子 泰勒展开法,计算截断误差的阶
差分方程关于右端稳定 保证误差不会被放大,是收敛性的关键 建立解的先验估计式

2. 先验估计的概念与意义

定理中提到的"形如(2.18)的估计式"称为先验估计,它是数值分析中最核心的技巧之一。

  • 定义:在不知道差分方程解的具体表达式的情况下,仅利用方程本身的性质,用右端项和边界条件的范数来估计解的范数。
  • 本质:先验估计就是稳定性的数学表达。
  • 意义:它不仅是证明收敛性的基础,也是分析误差、设计高效算法的重要工具。

3. 相容性与稳定性的分工

  • 相容性:回答"差分方程是否逼近原微分方程"的问题,是收敛性的必要条件。一个不相容的差分格式,无论多么稳定,都不可能收敛到正确的解。
  • 稳定性:回答"差分方程本身是否良定"的问题,是收敛性的充分条件。一个不稳定的差分格式,即使相容,也会因为误差的无限放大而失去实用价值。

四、稳定性的理论与实际意义

1. 理论意义

  1. 保证差分方程的适定性:稳定的差分方程对任何右端项和边界条件都存在唯一解。
  2. 是收敛性的必要条件:对于线性偏微分方程的数值解法,著名的Lax等价定理指出:相容且稳定的差分格式是收敛的。这意味着对于线性问题,稳定性不仅是收敛性的充分条件,也是必要条件。
  3. 是误差分析的基础:所有的误差估计最终都依赖于稳定性的先验估计式。

2. 实际工程意义

稳定性是差分格式能够实际应用的必要前提,原因如下:

  1. 实测误差不可避免:在工程问题中,微分方程的右端项 \(f\) 和边界条件 \(\alpha,\beta\) 通常是通过实验测量得到的,必然存在测量误差。
  2. 舍入误差不可避免:计算机采用有限精度运算,每一步计算都会产生舍入误差。
  3. 不稳定格式的灾难性后果:如果差分格式不稳定,这些微小的误差会在计算过程中被指数级放大,最终导致数值解完全失真,与精确解毫无关系。

结论:一个不稳定的差分格式,无论其精度多高,都没有任何实际应用价值。

五、定理在中心差分格式中的应用

我们用之前构造的二阶常微分方程中心差分格式来验证定理2.1的应用:

  1. 解的光滑性:若 \(q,f \in C[a,b]\),则原边值问题的解 \(u \in C^4[a,b]\),满足充分光滑的条件。
  2. 相容性验证:通过泰勒展开证明,中心差分格式的截断误差为 \(R_i(u) = O(h^2)\),因此在 \(\|\cdot\|_0\)\(\|\cdot\|_C\) 范数下都是二阶相容的。
  3. 稳定性验证:中心差分格式的系数矩阵是对称严格对角占优矩阵,因此是稳定的,可以建立先验估计式。
  4. 收敛性结论:根据定理2.1,中心差分格式的差分解按 \(\|\cdot\|_0\)\(\|\cdot\|_C\) 范数二阶收敛到精确解,即:

    \[\|u_h - u\|_C = O(h^2), \quad \|u_h - u\|_0 = O(h^2) \]

六、核心知识点归纳总结表

核心内容 关键结论 证明/依据 重要意义
定理2.1 相容 + 稳定 = 收敛,且收敛阶等于截断误差阶 误差方程 + 稳定性定义 + 相容性定义 椭圆型方程差分法的理论基石
先验估计 用右端项范数估计解的范数 差分方程的结构性质 稳定性的数学表达,误差分析的基础
相容性 差分算子逼近微分算子 泰勒中值定理 收敛性的必要条件
稳定性 解连续依赖于右端项 严格对角占优矩阵性质 收敛性的充分条件,格式实用的前提
中心差分格式 二阶收敛 定理2.1 最简单、最常用的二阶精度差分格式

2.2.1 直接差分化(变系数椭圆型方程非均匀网格格式)深度讲解

一、模型问题回顾

本节针对一般二阶自伴椭圆型方程第一边值问题构造差分格式,其形式为:

\[\begin{cases} Lu = -\frac{d}{dx}\left(p(x)\frac{du}{dx}\right) + r(x)\frac{du}{dx} + q(x)u = f(x), & a < x < b \\ u(a) = \alpha, \quad u(b) = \beta \end{cases} \]

其中 \(p(x) \geq p_{\text{min}} > 0\)(保证方程椭圆性),\(p(x), r(x), q(x), f(x)\) 均为 \([a,b]\) 上的连续函数。

重要性:这是比常系数方程更一般的模型,广泛出现在热传导、弹性力学、流体力学等工程领域。其差分格式的构造思想是处理所有变系数偏微分方程的基础。

二、非均匀网格与对偶剖分

1. 原网格剖分

将区间 \([a,b]\) 划分为 \(N\)任意长度的小区间,定义:

  • 节点:\(a = x_0 < x_1 < \dots < x_i < \dots < x_N = b\)
  • 局部步长:\(h_i = x_i - x_{i-1}, \quad i=1,2,\dots,N\)
  • 最大步长:\(h = \max_{1 \leq i \leq N} h_i\)
  • 内部节点集合:\(I_h = \{x_1, x_2, \dots, x_{N-1}\}\)
  • 全体节点集合:\(\bar{I}_h = \{x_0, x_1, \dots, x_N\}\)

2. 对偶剖分与半整数点

取相邻节点 \(x_{i-1}\)\(x_i\) 的中点:

\[x_{i-\frac{1}{2}} = \frac{1}{2}(x_{i-1} + x_i), \quad i=1,2,\dots,N \]

称为半整数点。由这些半整数点构成的剖分:

\[a = x_0 < x_{\frac{1}{2}} < x_{\frac{3}{2}} < \dots < x_{N-\frac{1}{2}} < x_N = b \]

称为对偶剖分

核心技巧解读:半整数点是构造变系数方程高精度差分格式的关键。对于变系数 \(p(x)\),直接在整数点 \(x_i\) 处用中心差分近似 \(\frac{d}{dx}\left(p\frac{du}{dx}\right)\) 只能得到一阶精度,而通过半整数点过渡可以在均匀网格下得到二阶精度。

三、变系数导数的差分近似(核心推导)

推导依据:泰勒中值定理
所有差分近似均通过将函数在节点处作泰勒展开,消去高阶项得到。

1. 一阶导数的非均匀差分近似

\(u(x_{i+1})\)\(u(x_{i-1})\)\(x_i\) 处作泰勒展开:

\[u(x_{i+1}) = u(x_i) + h_{i+1} u'(x_i) + \frac{h_{i+1}^2}{2} u''(x_i) + O(h_{i+1}^3) \]

\[u(x_{i-1}) = u(x_i) - h_i u'(x_i) + \frac{h_i^2}{2} u''(x_i) + O(h_i^3) \]

两式相减,消去 \(u(x_i)\) 项:

\[u(x_{i+1}) - u(x_{i-1}) = (h_i + h_{i+1}) u'(x_i) + \frac{h_{i+1}^2 - h_i^2}{2} u''(x_i) + O(h^3) \]

整理得一阶导数的表达式:

\[\left[\frac{du}{dx}\right]_i = \frac{u(x_{i+1}) - u(x_{i-1})}{h_i + h_{i+1}} - \frac{h_{i+1} - h_i}{2} u''(x_i) + O(h^2) \tag{2.22} \]

这就是非均匀网格下一阶导数的中心差分近似,截断误差为 \(O(h)\)(当网格不均匀时)。

2. 二阶自伴导数 \(\frac{d}{dx}\left(p\frac{du}{dx}\right)\) 的差分近似

这是变系数方程差分格式构造中最关键的一步,我们采用半整数点过渡法

步骤1:在半整数点处近似 \(p\frac{du}{dx}\)

\(u(x_i)\)\(x_{i-\frac{1}{2}}\) 处作泰勒展开:

\[u(x_i) = u(x_{i-\frac{1}{2}}) + \frac{h_i}{2} u'(x_{i-\frac{1}{2}}) + \frac{h_i^2}{8} u''(x_{i-\frac{1}{2}}) + \frac{h_i^3}{48} u'''(x_{i-\frac{1}{2}}) + O(h_i^4) \]

\[u(x_{i-1}) = u(x_{i-\frac{1}{2}}) - \frac{h_i}{2} u'(x_{i-\frac{1}{2}}) + \frac{h_i^2}{8} u''(x_{i-\frac{1}{2}}) - \frac{h_i^3}{48} u'''(x_{i-\frac{1}{2}}) + O(h_i^4) \]

两式相减,整理得:

\[u'(x_{i-\frac{1}{2}}) = \frac{u(x_i) - u(x_{i-1})}{h_i} - \frac{h_i^2}{24} u'''(x_{i-\frac{1}{2}}) + O(h_i^4) \]

两边乘以 \(p(x_{i-\frac{1}{2}})\),得到:

\[p(x_{i-\frac{1}{2}}) \frac{u(x_i) - u(x_{i-1})}{h_i} = \left[p\frac{du}{dx}\right]_{i-\frac{1}{2}} + \frac{h_i^2}{24} \left[p\frac{d^3 u}{dx^3}\right]_{i-\frac{1}{2}} + O(h_i^3) \tag{2.23} \]

同理,在半整数点 \(x_{i+\frac{1}{2}}\) 处有:

\[p(x_{i+\frac{1}{2}}) \frac{u(x_{i+1}) - u(x_i)}{h_{i+1}} = \left[p\frac{du}{dx}\right]_{i+\frac{1}{2}} + \frac{h_{i+1}^2}{24} \left[p\frac{d^3 u}{dx^3}\right]_{i+\frac{1}{2}} + O(h_{i+1}^3) \tag{2.24} \]

步骤2:在整数点处近似二阶导数

二阶自伴导数 \(\frac{d}{dx}\left(p\frac{du}{dx}\right)\) 是函数 \(w(x) = p(x)\frac{du}{dx}\) 的一阶导数,因此可以用一阶导数的差分近似来处理。

将式(2.24)减去式(2.23),并除以平均步长 \(\frac{h_i + h_{i+1}}{2}\),得到:

\[\frac{2}{h_i + h_{i+1}} \left[ p(x_{i+\frac{1}{2}}) \frac{u(x_{i+1}) - u(x_i)}{h_{i+1}} - p(x_{i-\frac{1}{2}}) \frac{u(x_i) - u(x_{i-1})}{h_i} \right] \]

\[= \frac{2}{h_i + h_{i+1}} \left( \left[p\frac{du}{dx}\right]_{i+\frac{1}{2}} - \left[p\frac{du}{dx}\right]_{i-\frac{1}{2}} \right) + O(h^2) \]

再将 \(\left[p\frac{du}{dx}\right]_{i+\frac{1}{2}}\)\(\left[p\frac{du}{dx}\right]_{i-\frac{1}{2}}\)\(x_i\) 处作泰勒展开:

\[\left[p\frac{du}{dx}\right]_{i+\frac{1}{2}} = \left[p\frac{du}{dx}\right]_i + \frac{h_{i+1}}{2} \left[\frac{d}{dx}\left(p\frac{du}{dx}\right)\right]_i + \frac{h_{i+1}^2}{8} \left[\frac{d^2}{dx^2}\left(p\frac{du}{dx}\right)\right]_i + O(h_{i+1}^3) \]

\[\left[p\frac{du}{dx}\right]_{i-\frac{1}{2}} = \left[p\frac{du}{dx}\right]_i - \frac{h_i}{2} \left[\frac{d}{dx}\left(p\frac{du}{dx}\right)\right]_i + \frac{h_i^2}{8} \left[\frac{d^2}{dx^2}\left(p\frac{du}{dx}\right)\right]_i + O(h_i^3) \]

两式相减,整理后最终得到:

\[\frac{2}{h_i + h_{i+1}} \left[ p(x_{i+\frac{1}{2}}) \frac{u(x_{i+1}) - u(x_i)}{h_{i+1}} - p(x_{i-\frac{1}{2}}) \frac{u(x_i) - u(x_{i-1})}{h_i} \right] \]

\[= \left[\frac{d}{dx}\left(p\frac{du}{dx}\right)\right]_i + \frac{h_{i+1} - h_i}{4} \left[\frac{d^2}{dx^2}\left(p\frac{du}{dx}\right)\right]_i + O(h^2) \tag{2.25} \]

这就是二阶自伴导数的非均匀网格差分近似。

四、完整差分格式的建立

1. 差分方程的推导

将式(2.22)和式(2.25)代入原微分方程(2.20),并引入记号:

\[p_{i-\frac{1}{2}} = p(x_{i-\frac{1}{2}}), \quad r_i = r(x_i), \quad q_i = q(x_i), \quad f_i = f(x_i) \]

得到精确解满足的方程:

\[L_h u(x_i) = f_i + R_i(u) \tag{2.26} \]

其中差分算子 \(L_h\) 定义为:

\[L_h u(x_i) \equiv -\frac{2}{h_i + h_{i+1}} \left[ p_{i+\frac{1}{2}} \frac{u(x_{i+1}) - u(x_i)}{h_{i+1}} - p_{i-\frac{1}{2}} \frac{u(x_i) - u(x_{i-1})}{h_i} \right] + \frac{r_i}{h_i + h_{i+1}} \left[ u(x_{i+1}) - u(x_{i-1}) \right] + q_i u(x_i) \]

\(R_i(u)\) 是截断误差,其表达式为:

\[R_i(u) = -(h_{i+1} - h_i) \left( \frac{1}{4} \left[\frac{d^2}{dx^2}\left(p\frac{du}{dx}\right)\right]_i + \frac{1}{12} \left[p\frac{d^3 u}{dx^3}\right]_i - \frac{1}{2} \left[r\frac{d^2 u}{dx^2}\right]_i \right) + O(h^2) \tag{2.27} \]

舍去截断误差项 \(R_i(u)\),得到逼近原微分方程的差分方程:

\[L_h u_i \equiv -\frac{2}{h_i + h_{i+1}} \left[ p_{i+\frac{1}{2}} \frac{u_{i+1} - u_i}{h_{i+1}} - p_{i-\frac{1}{2}} \frac{u_i - u_{i-1}}{h_i} \right] + \frac{r_i}{h_i + h_{i+1}} (u_{i+1} - u_{i-1}) + q_i u_i = f_i \tag{2.28} \]

结合边界条件,得到完整的差分方程组:

\[\begin{cases} L_h u_i = f_i, & i=1,2,\dots,N-1 \\ u_0 = \alpha, \quad u_N = \beta \tag{2.29} \end{cases} \]

五、截断误差分析

1. 非均匀网格情形

从截断误差表达式(2.27)可以看出,当网格不均匀时,\(h_{i+1} - h_i = O(h)\),因此:

\[\|R_h(u)\|_C = O(h), \quad \|R_h(u)\|_0 = O(h) \]

即非均匀网格下的差分格式是一阶精度的。

2. 均匀网格情形

当网格均匀时,\(h_i = h_{i+1} = h\),此时 \(h_{i+1} - h_i = 0\),式(2.27)中的一阶项消失,截断误差提高为:

\[\|R_h(u)\|_C = O(h^2), \quad \|R_h(u)\|_0 = O(h^2) \]

即均匀网格下的差分格式是二阶精度的。

关键结论:网格的均匀性对变系数方程差分格式的精度有决定性影响。在实际计算中,应尽量采用均匀网格;若必须采用非均匀网格(如解变化剧烈的区域需要加密),则应保证网格变化缓慢,以减小截断误差。

六、差分方程组的系数矩阵性质

将差分方程组(2.28)-(2.29)整理为线性代数方程组 \(A\mathbf{u} = \mathbf{b}\),其系数矩阵 \(A\) 具有以下性质:

  1. 三对角性:仅在主对角线和相邻的两条次对角线上有非零元素,其余元素均为0。
  2. 稀疏性:非零元素个数约为 \(3(N-1)\),计算效率高。
  3. 不对称性
    • \(r(x) \not\equiv 0\) 时,由于一阶导数项的存在,矩阵 \(A\) 不对称;
    • \(r(x) \equiv 0\) 时,若网格不均匀,矩阵 \(A\) 仍然不对称。
  4. 可对称化性:在差分方程(2.28)两端同时乘以 \((h_i + h_{i+1})\),得到的新方程组的系数矩阵是对称的。这一技巧在求解时非常有用,因为对称矩阵有更高效的求解算法。

七、均匀网格下的简化格式

当网格均匀时,\(h_i = h_{i+1} = h\),差分格式(2.28)简化为:

\[L_h u_i = -\frac{1}{h^2} \left[ p_{i+\frac{1}{2}} u_{i+1} - (p_{i+\frac{1}{2}} + p_{i-\frac{1}{2}}) u_i + p_{i-\frac{1}{2}} u_{i-1} \right] + r_i \frac{u_{i+1} - u_{i-1}}{2h} + q_i u_i = f_i \tag{2.31} \]

这一格式具有清晰的物理意义:

  • 第一项是二阶导数项的中心差分近似;
  • 第二项是一阶导数项的中心差分近似;
  • 第三项是零阶项的直接近似。

\(p(x) \equiv 1\)\(r(x) \equiv 0\)\(q(x) \equiv q\) 时,该格式退化为我们之前讨论的常系数中心差分格式,验证了其正确性。

八、核心知识点归纳总结表

核心知识点 关键内容 推导/依据 重要结论
模型问题 一般二阶自伴椭圆型方程 \(-\frac{d}{dx}(p\frac{du}{dx})+r\frac{du}{dx}+qu=f\) 椭圆型方程基本理论 涵盖了绝大多数工程中的一维椭圆问题
网格剖分 非均匀网格、对偶剖分、半整数点 区域分解原理 半整数点是构造变系数高精度格式的关键
一阶导数近似 \(\frac{u_{i+1}-u_{i-1}}{h_i+h_{i+1}}\) 泰勒中值定理 非均匀网格下一阶精度,均匀网格下二阶精度
二阶自伴导数近似 \(\frac{2}{h_i+h_{i+1}}\left(p_{i+\frac{1}{2}}\frac{u_{i+1}-u_i}{h_{i+1}}-p_{i-\frac{1}{2}}\frac{u_i-u_{i-1}}{h_i}\right)\) 泰勒中值定理+半整数点过渡 非均匀网格下一阶精度,均匀网格下二阶精度
截断误差 非均匀网格 \(O(h)\),均匀网格 \(O(h^2)\) 泰勒展开余项分析 网格均匀性对精度有决定性影响
系数矩阵 三对角、稀疏、不对称但可对称化 差分方程结构 可通过乘以 \((h_i+h_{i+1})\) 实现对称化
均匀网格简化格式 用一阶和二阶中心差商代替代数导数 非均匀格式的特例 形式简单,精度最高,是最常用的格式

2.2.2 有限体积法(积分插值法)深度讲解

一、有限体积法的本质与核心优势

1. 模型问题:守恒型微分方程

本节针对自伴守恒型二阶椭圆型方程构造差分格式:

\[Lu = -\frac{d}{dx}\left(p(x)\frac{du}{dx}\right) + q(x)u = f(x) \tag{2.32} \]

其中 \(p(x) \geq p_{\text{min}} > 0\)(保证椭圆性)。该方程是自然界中所有守恒律问题(热传导、质量扩散、流体流动、静电场等)的数学抽象。

2. 积分守恒形式的推导(物理意义:热量守恒)

推导依据:牛顿-莱布尼茨公式
将微分方程(2.32)在任意区间 \([x^{(1)}, x^{(2)}]\) 上积分:

\[-\int_{x^{(1)}}^{x^{(2)}} \frac{d}{dx}\left(p(x)\frac{du}{dx}\right) dx + \int_{x^{(1)}}^{x^{(2)}} q(x)u dx = \int_{x^{(1)}}^{x^{(2)}} f(x) dx \]

由牛顿-莱布尼茨公式,第一项可化简为边界上的通量差:

\[W(x^{(1)}) - W(x^{(2)}) + \int_{x^{(1)}}^{x^{(2)}} q(x)u dx = \int_{x^{(1)}}^{x^{(2)}} f(x) dx \tag{2.33} \]

其中 \(W(x) = p(x)\frac{du}{dx}\) 称为通量函数(在热传导问题中就是热流量)。

积分形式的核心优势

  • 最高阶导数从二阶降为一阶,显著减弱了对解 \(u\) 和系数 \(p\) 的光滑性要求
  • 直接反映了物理守恒律(流入控制体的通量减去流出的通量,加上内部源项,等于控制体内物理量的变化率);
  • 天然适合处理间断系数问题和复杂边界条件。

二、守恒型差分格式的完整推导

1. 控制体积的选取:对偶单元

有限体积法的核心思想是在每个控制体积上严格满足守恒律。我们选取以内部节点 \(x_i\) 为中心的对偶单元

\[V_i = \left[x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}\right] \]

作为第 \(i\) 个控制体积。这与直接差分化法以 \([x_{i-1}, x_i]\) 为单元的做法不同,是有限体积法的标志性特征。

在控制体积 \(V_i\) 上应用积分守恒方程(2.33),得到:

\[W\left(x_{i-\frac{1}{2}}\right) - W\left(x_{i+\frac{1}{2}}\right) + \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} q(x)u dx = \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} f(x) dx \tag{2.35} \]

2. 通量函数 \(W(x)\) 的差分近似(最关键步骤)

推导依据:通量连续性原理
即使系数 \(p(x)\) 存在间断点,通量函数 \(W(x) = p(x)\frac{du}{dx}\) 仍然是连续的(物理上,热流量、质量流量等不能突变)。这是有限体积法能正确处理间断系数问题的根本原因。

\(\frac{du}{dx} = \frac{W(x)}{p(x)}\) 在区间 \([x_{i-1}, x_i]\) 上积分:

\[u_i - u_{i-1} = \int_{x_{i-1}}^{x_i} \frac{W(x)}{p(x)} dx \]

推导依据:积分中矩形公式
用中矩形公式近似右边的积分(二阶精度):

\[\int_{x_{i-1}}^{x_i} \frac{W(x)}{p(x)} dx \approx W\left(x_{i-\frac{1}{2}}\right) \int_{x_{i-1}}^{x_i} \frac{dx}{p(x)} \]

整理得到通量的近似表达式:

\[W_{i-\frac{1}{2}} \approx a_i \frac{u_i - u_{i-1}}{h_i} \tag{2.36} \]

其中

\[a_i = \left( \frac{1}{h_i} \int_{x_{i-1}}^{x_i} \frac{dx}{p(x)} \right)^{-1} \tag{2.37} \]

关键解读\(a_i\)\(p(x)\) 在区间 \([x_{i-1}, x_i]\) 上的调和平均,而非算术平均。这是有限体积法与直接差分法最本质的区别之一,也是其能正确处理间断系数的核心技巧。

3. 反应项与源项的积分近似

推导依据:积分中矩形公式
在控制体积 \(V_i\) 上对反应项 \(q(x)u\) 和源项 \(f(x)\) 应用中矩形公式(二阶精度):

\[\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} q(x)u dx \approx \frac{h_i + h_{i+1}}{2} d_i u_i \tag{2.38} \]

\[\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} f(x) dx \approx \frac{h_i + h_{i+1}}{2} \phi_i \]

其中

\[d_i = \frac{2}{h_i + h_{i+1}} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} q(x) dx \tag{2.39} \]

\[\phi_i = \frac{2}{h_i + h_{i+1}} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} f(x) dx \tag{2.41} \]

4. 守恒型差分方程的建立

将通量近似(2.36)和源项近似代入积分守恒方程(2.35),得到:

\[-\left( a_{i+1} \frac{u_{i+1} - u_i}{h_{i+1}} - a_i \frac{u_i - u_{i-1}}{h_i} \right) + \frac{1}{2}(h_i + h_{i+1}) d_i u_i = \frac{1}{2}(h_i + h_{i+1}) \phi_i \tag{2.40} \]

核心性质:该差分格式在每个控制体积上都严格满足物理守恒律,因此在整个求解域上也满足守恒律。这是有限体积法最宝贵的性质,也是其在计算流体力学等领域占据主导地位的根本原因。

三、系数的实用计算方法

1. 中矩形公式(光滑系数情形)

当系数 \(p, q\) 和右端项 \(f\) 充分光滑时,用中矩形公式计算积分(2.37)(2.39)(2.41),得到:

\[\begin{cases} a_i = p_{i-\frac{1}{2}} = p\left(x_{i-\frac{1}{2}}\right) \\ d_i = q_i = q(x_i) \\ \phi_i = f_i = f(x_i) \end{cases} \tag{2.42} \]

此时,在均匀网格(\(h_i = h_{i+1} = h\))下,有限体积法的差分格式与直接差分化法的格式完全等价,都具有二阶精度。

2. 梯形公式(一般情形)

当函数不够光滑或存在间断时,用梯形公式计算积分,得到:

\[\begin{cases} a_i = \frac{2 p_{i-1} p_i}{p_{i-1} + p_i} \\ d_i = \frac{1}{2}\left(q_{i-\frac{1}{2}} + q_{i+\frac{1}{2}}\right) \\ \phi_i = \frac{1}{2}\left(f_{i-\frac{1}{2}} + f_{i+\frac{1}{2}}\right) \end{cases} \tag{2.43} \]

关键结论:用梯形公式计算的 \(a_i\) 正是 \(p(x)\) 在区间 \([x_{i-1}, x_i]\) 上的调和平均。当 \(p(x)\) 存在间断点时,调和平均是唯一正确的平均方式,而算术平均会导致数值解不收敛。

四、有限体积法的独特优势:间断系数问题的正确处理

1. 问题背景

在实际工程问题中,系数 \(p(x)\) 经常存在间断点(例如不同材料的交界面、不同介质的分界面)。此时,微分方程的解 \(u(x)\) 连续,但一阶导数 \(\frac{du}{dx}\) 存在跳变。

2. 直接差分法的致命缺陷

如果将守恒型方程(2.32)改写为非守恒形式:

\[p(x) \frac{d^2 u}{dx^2} + p'(x) \frac{du}{dx} + q(x)u = f(x) \]

然后用中心差分格式近似:

\[p_i \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} + \frac{p_{i+1} - p_{i-1}}{2h} \frac{u_{i+1} - u_{i-1}}{2h} + q_i u_i = f_i \]

则在间断点处,\(p'(x)\) 是无穷大,其差分近似 \(\frac{p_{i+1} - p_{i-1}}{2h}\) 会产生极大的误差,导致差分解不收敛到精确解。

3. 有限体积法的正确性

有限体积法从积分守恒形式出发,完全避免了对 \(p'(x)\) 的近似。它利用通量 \(W(x)\) 的连续性,通过调和平均自动处理了间断点,保证了差分解收敛到正确的精确解

典型例子:考虑方程 \(\frac{d}{dx}\left(p(x)\frac{du}{dx}\right) = 0\),其中 \(p(x)\)\(x=0.5\) 处有间断:\(p(x)=1\)\(x<0.5\)),\(p(x)=100\)\(x>0.5\))。精确解为 \(u(x)=Cx\)\(x<0.5\)),\(u(x)=0.5C + \frac{C}{100}(x-0.5)\)\(x>0.5\))。有限体积法能准确捕捉到导数的跳变,而直接差分法会得到完全错误的结果。

五、直接差分化法与有限体积法的全面对比

对比维度 直接差分化法 有限体积法(积分插值法)
出发点 微分方程的微分形式 微分方程的积分守恒形式
数学基础 泰勒展开、数值微商 积分守恒律、数值积分
核心思想 用差商代替微商 在每个控制体积上满足守恒律
守恒性 仅在均匀网格和常系数下满足 天然满足,与网格和系数无关
间断系数处理 会产生极大误差,甚至不收敛 自动正确处理,保证收敛
光滑性要求 要求系数 \(p(x)\) 连续可导 仅要求通量 \(W(x)\) 连续
精度 均匀网格下二阶精度 均匀网格下二阶精度,非均匀网格下一阶精度
物理意义 不明确 非常明确,直接对应物理守恒律
推广难度 难以推广到非结构网格 极易推广到任意网格和复杂边界条件

六、核心知识点归纳总结表

核心知识点 关键内容 推导/依据 重要结论
积分守恒形式 \(W(x^{(1)})-W(x^{(2)})+\int qudx=\int fdx\) 牛顿-莱布尼茨公式 降低光滑性要求,反映物理守恒律
控制体积 对偶单元 \([x_{i-1/2},x_{i+1/2}]\) 区域分解原理 有限体积法的标志性特征
通量近似 \(W_{i-1/2} \approx a_i \frac{u_i-u_{i-1}}{h_i}\) 通量连续性原理+中矩形公式 \(a_i\)\(p(x)\) 的调和平均
守恒型差分格式 式(2.40) 积分守恒方程+数值积分 严格满足守恒律
系数计算 中矩形公式(光滑)、梯形公式(一般) 数值积分公式 梯形公式自动给出调和平均
间断系数处理 不收敛 正确收敛 有限体积法的最大优势
等价性 均匀光滑网格下与直接差分法等价 格式对比 两种方法在简单情形下结果一致

2.2.3 边值条件的处理 深度讲解

一、边值条件的类型与简单差分法的缺陷

1. 第二、第三边值条件的一般形式

我们已经处理了最简单的第一边值条件(Dirichlet条件)\(u(a)=\alpha, u(b)=\beta\)。现在处理更复杂的第二边值条件(Neumann条件)第三边值条件(Robin条件),它们的统一形式为:

\[\begin{cases} u'(a) = \alpha_0 u(a) + \alpha_1 \\ u'(b) = \beta_0 u(b) + \beta_1 \end{cases} \]

  • \(\alpha_0=0, \beta_0=0\) 时,退化为第二边值条件(给定边界上的导数值);
  • \(\alpha_0 \neq 0\)\(\beta_0 \neq 0\) 时,为第三边值条件(边界上的导数与函数值线性相关)。

由于 \(p(x) > 0\),为了与有限体积法的通量形式统一,我们将边界条件改写为:

\[\begin{cases} -p(a)u'(a) = \alpha_0 u(a) + \alpha_1 \\ -p(b)u'(b) = \beta_0 u(b) + \beta_1 \end{cases} \]

其中左边正是边界处的通量 \(W(a) = p(a)u'(a)\) 的负值。

2. 简单数值微商法的两个致命缺点

最直观的处理方法是用向前差分近似左边界导数,向后差分近似右边界导数:

\[u'(a) \approx \frac{u_1 - u_0}{h_1}, \quad u'(b) \approx \frac{u_N - u_{N-1}}{h_N} \]

但这种方法存在两个不可忽视的缺陷:

  1. 精度损失:内点差分格式的截断误差为 \(O(h^2)\),而边界处的向前/向后差分截断误差仅为 \(O(h)\),这会导致整个数值解的整体精度从二阶降为一阶。
  2. 破坏对称性:简单差分法导出的边界方程会破坏差分方程组系数矩阵的对称性,而对称性是使用共轭梯度法等高效求解算法的前提,会显著增加计算量。

二、有限体积法处理边界条件的核心思想

有限体积法处理边界条件的核心优势在于:将边界点也作为控制体积的中心,在边界控制体积上应用积分守恒律,自然导出与内点格式精度一致、且保持对称性的边界差分方程

具体来说:

  • 左边界点 \(x_0=a\) 的控制体积为半单元:\(V_0 = \left[x_0, x_{\frac{1}{2}}\right]\)
  • 右边界点 \(x_N=b\) 的控制体积为半单元:\(V_N = \left[x_{N-\frac{1}{2}}, x_N\right]\)

这种方法的本质是将边界条件作为控制体积的通量边界条件,融入到守恒方程中,而不是作为独立的代数条件强加给差分方程组。

三、左边界差分方程的完整推导

推导依据:积分守恒律 + 边界条件 + 数值积分公式

1. 积分守恒方程的应用

在左边界控制体积 \(V_0 = \left[x_0, x_{\frac{1}{2}}\right]\) 上应用积分守恒方程(2.33):

\[W(x_0) - W\left(x_{\frac{1}{2}}\right) + \int_{x_0}^{x_{\frac{1}{2}}} q(x)u dx = \int_{x_0}^{x_{\frac{1}{2}}} f(x) dx \]

2. 边界条件的代入

根据改写后的边界条件(2.46),左边界处的通量为:

\[W(x_0) = p(a)u'(a) = -(\alpha_0 u_0 + \alpha_1) \]

代入积分守恒方程,得到:

\[-(\alpha_0 u_0 + \alpha_1) - W\left(x_{\frac{1}{2}}\right) + \int_{x_0}^{x_{\frac{1}{2}}} q(x)u dx = \int_{x_0}^{x_{\frac{1}{2}}} f(x) dx \]

整理为:

\[-W\left(x_{\frac{1}{2}}\right) + \int_{x_0}^{x_{\frac{1}{2}}} q(x)u dx = (\alpha_0 u_0 + \alpha_1) + \int_{x_0}^{x_{\frac{1}{2}}} f(x) dx \tag{2.48} \]

3. 通量与积分项的近似

  • 通量近似:内点处的通量近似公式(2.36)同样适用于半整数点 \(x_{\frac{1}{2}}\)

    \[W\left(x_{\frac{1}{2}}\right) \approx a_1 \frac{u_1 - u_0}{h_1}, \quad a_1 = \left( \frac{1}{h_1} \int_{x_0}^{x_1} \frac{dx}{p(x)} \right)^{-1} \tag{2.49} \]

  • 积分项近似:用中矩形公式近似控制体积上的积分(二阶精度):

    \[\int_{x_0}^{x_{\frac{1}{2}}} q(x)u dx \approx \frac{h_1}{2} d_0 u_0, \quad d_0 = \frac{2}{h_1} \int_{x_0}^{x_{\frac{1}{2}}} q(x) dx \tag{2.50} \]

    \[\int_{x_0}^{x_{\frac{1}{2}}} f(x) dx \approx \frac{h_1}{2} \phi_0, \quad \phi_0 = \frac{2}{h_1} \int_{x_0}^{x_{\frac{1}{2}}} f(x) dx \tag{2.51} \]

4. 边界差分方程的最终形式

将式(2.49)-(2.51)代入式(2.48),得到:

\[-a_1 \frac{u_1 - u_0}{h_1} + \frac{h_1}{2} d_0 u_0 = \alpha_0 u_0 + \alpha_1 + \frac{h_1}{2} \phi_0 \]

整理为标准形式:

\[-a_1 \frac{u_1 - u_0}{h_1} + \left( -\alpha_0 + \frac{h_1}{2} d_0 \right) u_0 - \left( \alpha_1 + \frac{h_1}{2} \phi_0 \right) = 0 \tag{2.52} \]

四、右边界差分方程的推导

采用完全相同的方法,在右边界控制体积 \(V_N = \left[x_{N-\frac{1}{2}}, x_N\right]\) 上应用积分守恒律,代入右边界条件(2.47),可以导出右边界的差分方程:

\[a_N \frac{u_N - u_{N-1}}{h_N} + \left( -\beta_0 + \frac{h_N}{2} d_N \right) u_N - \left( \beta_1 + \frac{h_N}{2} \phi_N \right) = 0 \]

其中

\[a_N = \left( \frac{1}{h_N} \int_{x_{N-1}}^{x_N} \frac{dx}{p(x)} \right)^{-1}, \quad d_N = \frac{2}{h_N} \int_{x_{N-\frac{1}{2}}}^{x_N} q(x) dx, \quad \phi_N = \frac{2}{h_N} \int_{x_{N-\frac{1}{2}}}^{x_N} f(x) dx \]

五、精度与对称性分析

1. 二阶精度证明

定理:当网格均匀且系数光滑时,边界差分方程(2.52)逼近原边界条件(2.44)的截断误差为 \(O(h^2)\),与内点差分格式的精度一致。

证明依据:泰勒展开余项分析

  • 通量近似 \(W(x_{\frac{1}{2}})\) 的截断误差为 \(O(h^2)\)
  • 中矩形公式的积分近似截断误差为 \(O(h^2)\)
  • 边界条件的代入是精确的。

因此,边界差分方程的整体截断误差为 \(O(h^2)\),保证了整个数值解的二阶精度。

2. 对称性保持

有限体积法导出的边界差分方程与内点方程具有完全一致的结构,因此整个差分方程组的系数矩阵仍然保持对称性。这是该方法最重要的优点之一,因为对称正定矩阵可以使用共轭梯度法等高效迭代算法,计算量比非对称矩阵小一个数量级以上。

六、特殊情况:纯Neumann边界条件

\(\alpha_0=0, \beta_0=0\) 时,边界条件退化为纯Neumann条件:

\[u'(a) = \alpha_1, \quad u'(b) = \beta_1 \]

此时原微分方程的解不唯一,相差一个任意常数,解存在的充要条件是:

\[\int_a^b f(x) dx = \beta_1 p(b) - \alpha_1 p(a) \]

相应地,差分方程组也存在解不唯一的问题,需要附加一个归一化条件(例如指定某一点的函数值,或要求解的平均值为0)才能得到唯一解。有限体积法处理纯Neumann边界条件时,仍然保持守恒性和二阶精度。

七、两种边界处理方法的全面对比

对比维度 简单数值微商法 有限体积法
出发点 直接用差商代替边界导数 在边界控制体积上应用积分守恒律
截断误差 \(O(h)\)(一阶精度) \(O(h^2)\)(二阶精度,与内点一致)
对称性 破坏系数矩阵的对称性 保持系数矩阵的对称性
守恒性 不满足边界控制体积的守恒性 严格满足所有控制体积的守恒性
间断系数处理 误差大,可能不收敛 自动正确处理
实现难度 简单 稍复杂,但与内点格式统一
整体精度 一阶 二阶

八、核心知识点归纳总结表

核心知识点 关键内容 推导/依据 重要结论
边值条件类型 第一类(Dirichlet)、第二类(Neumann)、第三类(Robin) 偏微分方程定解理论 第二、三类边界条件处理难度更大
简单差分法缺陷 精度损失(一阶)、破坏对称性 泰勒展开余项分析 会导致整体精度下降和计算效率降低
有限体积法思想 构造边界半控制体积,应用积分守恒律 守恒律的普适性 自然融入边界条件,保持格式一致性
边界差分方程 式(2.52)及右边界对应方程 积分守恒律+边界条件+数值积分 与内点格式精度一致,保持对称性
精度分析 边界格式截断误差 \(O(h^2)\) 泰勒展开 保证整个数值解的二阶精度
对称性保持 系数矩阵仍然对称 格式结构一致性 可使用高效的对称矩阵求解算法
纯Neumann条件 解不唯一,需附加归一化条件 微分方程解的存在唯一性理论 有限体积法仍保持守恒性和精度

2.2.4 习题 完整解答与推导

(资深高级研究员视角,严格遵循教材理论体系,补充关键推导依据与技巧说明)


习题1:用有限体积法导出逼近微分方程(2.20)的差分方程

1. 模型问题回顾

微分方程(2.20)为一般二阶自伴椭圆型方程(含对流项)

\[Lu = -\frac{d}{dx}\left(p(x)\frac{du}{dx}\right) + r(x)\frac{du}{dx} + q(x)u = f(x), \quad a < x < b \]

其中 \(p(x) \geq p_{\text{min}} > 0\)(保证椭圆性),\(p,r,q,f\) 均为 \([a,b]\) 上的连续函数。

2. 有限体积法核心步骤

推导依据:积分守恒律 + 数值积分公式

步骤1:控制体积选取

取以内部节点 \(x_i\) 为中心的对偶控制体积

\[V_i = \left[x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}\right] \]

其中半整数点 \(x_{i\pm\frac{1}{2}} = \frac{1}{2}(x_i + x_{i\pm1})\),控制体积长度为 \(\Delta x_i = x_{i+\frac{1}{2}} - x_{i-\frac{1}{2}} = \frac{h_i + h_{i+1}}{2}\)\(h_i = x_i - x_{i-1}\) 为局部步长。

步骤2:方程积分(守恒形式建立)

将微分方程在控制体积 \(V_i\) 上积分:

\[-\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \frac{d}{dx}\left(p\frac{du}{dx}\right) dx + \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} r\frac{du}{dx} dx + \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} qu dx = \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} f dx \tag{1} \]

步骤3:各项积分近似

  • 扩散项(第一项):由牛顿-莱布尼茨公式化简为边界通量差

    \[-\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \frac{d}{dx}\left(p\frac{du}{dx}\right) dx = W\left(x_{i-\frac{1}{2}}\right) - W\left(x_{i+\frac{1}{2}}\right) \]

    其中通量 \(W(x) = p(x)\frac{du}{dx}\),采用半整数点调和平均近似(与守恒型方程一致):

    \[W_{i-\frac{1}{2}} \approx a_i \frac{u_i - u_{i-1}}{h_i}, \quad a_i = \left( \frac{1}{h_i} \int_{x_{i-1}}^{x_i} \frac{dx}{p(x)} \right)^{-1} \]

    该近似保证了通量的连续性,即使 \(p(x)\) 存在间断点也能正确处理。

  • 对流项(第二项):采用中心差分近似(二阶精度)
    利用中矩形公式,将 \(r(x)\) 在控制体积上取为 \(r_i = r(x_i)\),一阶导数用中心差分近似:

    \[\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} r\frac{du}{dx} dx \approx r_i \cdot \frac{u_{i+1} - u_{i-1}}{h_i + h_{i+1}} \cdot \frac{h_i + h_{i+1}}{2} = \frac{r_i}{2}(u_{i+1} - u_{i-1}) \]

  • 反应项与源项(第三、四项):采用中矩形公式近似(二阶精度)

    \[\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} qu dx \approx q_i u_i \cdot \frac{h_i + h_{i+1}}{2}, \quad \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} f dx \approx f_i \cdot \frac{h_i + h_{i+1}}{2} \]

    其中 \(q_i = q(x_i), f_i = f(x_i)\)

步骤4:差分方程整理

将所有近似代入式(1),得到:

\[a_i \frac{u_i - u_{i-1}}{h_i} - a_{i+1} \frac{u_{i+1} - u_i}{h_{i+1}} + \frac{r_i}{2}(u_{i+1} - u_{i-1}) + \frac{h_i + h_{i+1}}{2} q_i u_i = \frac{h_i + h_{i+1}}{2} f_i \]

两边乘以 \(\frac{2}{h_i + h_{i+1}}\),整理为标准形式:

\[\boxed{ -\frac{2}{h_i + h_{i+1}} \left( a_{i+1} \frac{u_{i+1} - u_i}{h_{i+1}} - a_i \frac{u_i - u_{i-1}}{h_i} \right) + \frac{r_i}{h_i + h_{i+1}}(u_{i+1} - u_{i-1}) + q_i u_i = f_i } \]

3. 结果说明

  • 该差分格式与直接差分化法得到的格式(2.28)完全一致,验证了两种方法在光滑系数下的等价性;
  • 对流项的存在导致系数矩阵不对称,这是对流扩散方程差分格式的普遍特征;
  • \(r(x) \equiv 0\) 时,退化为守恒型方程的有限体积格式(2.40),保持守恒性与对称性。

习题2:构造四阶椭圆型方程的中心差分格式

1. 模型问题回顾

四阶椭圆型边值问题:

\[\begin{cases} \left(p(x)u''\right)'' + \left(q(x)u'\right)' + r(x)u = f(x), & a < x < b \\ u(a) = u'(a) = 0, \quad u(b) = u'(b) = 0 \end{cases} \]

该方程是梁的弯曲振动、薄板弯曲等工程问题的数学模型。

2. 网格与差分近似基础

推导依据:泰勒中值定理 + 中心差分公式
采用均匀网格,步长 \(h = \frac{b-a}{N}\),节点 \(x_i = a + ih, \ i=0,1,\dots,N\)

核心差分公式(二阶精度)

  • 一阶导数:\(u'(x_i) \approx \frac{u_{i+1} - u_{i-1}}{2h}\)
  • 二阶导数:\(u''(x_i) \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2}\)

3. 各阶导数项的中心差分近似

步骤1:四阶导数项 \(\left(p(x)u''\right)''\) 的近似

采用逐层差分法,从低阶到高阶逐步近似:

  1. 先近似二阶导数 \(u''\)

    \[u''(x_i) \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} \]

  2. 定义中间变量 \(v(x) = p(x)u''(x)\),则 \(v(x_i) \approx p_i \cdot \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2}\),其中 \(p_i = p(x_i)\)
  3. 近似一阶导数 \(v'(x)\)

    \[v'\left(x_{i+\frac{1}{2}}\right) \approx \frac{v_{i+1} - v_i}{h} \]

  4. 近似二阶导数 \(v''(x) = \left(pu''\right)''\)

    \[v''(x_i) \approx \frac{v'\left(x_{i+\frac{1}{2}}\right) - v'\left(x_{i-\frac{1}{2}}\right)}{h} = \frac{v_{i+1} - 2v_i + v_{i-1}}{h^2} \]

\(v_i\) 代入,最终得到四阶导数项的中心差分近似:

\[\left(pu''\right)''(x_i) \approx \frac{1}{h^4} \left[ p_{i+1}(u_{i+2} - 2u_{i+1} + u_i) - 2p_i(u_{i+1} - 2u_i + u_{i-1}) + p_{i-1}(u_i - 2u_{i-1} + u_{i-2}) \right] \]

截断误差为 \(O(h^2)\)

步骤2:二阶导数项 \(\left(q(x)u'\right)'\) 的近似

与二阶自伴椭圆型方程的扩散项近似完全一致:

\[\left(q u'\right)'(x_i) \approx \frac{1}{h^2} \left[ q_{i+\frac{1}{2}}(u_{i+1} - u_i) - q_{i-\frac{1}{2}}(u_i - u_{i-1}) \right] \]

其中 \(q_{i\pm\frac{1}{2}} = q\left(x_{i\pm\frac{1}{2}}\right)\),截断误差为 \(O(h^2)\)

步骤3:零阶项与右端项

零阶项 \(r(x)u\) 和右端项 \(f(x)\) 直接在节点处取值:

\[r(x_i)u(x_i) \approx r_i u_i, \quad f(x_i) \approx f_i \]

其中 \(r_i = r(x_i), f_i = f(x_i)\)

4. 内部节点差分方程

将所有近似代入原微分方程,得到内部节点 \(i=2,3,\dots,N-2\) 的差分方程:

\[\boxed{ \begin{aligned} &\frac{p_{i+1}}{h^4} u_{i+2} + \left( -\frac{2p_{i+1} + 2p_i}{h^4} + \frac{q_{i+\frac{1}{2}}}{h^2} \right) u_{i+1} \\ +& \left( \frac{p_{i+1} + 4p_i + p_{i-1}}{h^4} - \frac{q_{i+\frac{1}{2}} + q_{i-\frac{1}{2}}}{h^2} + r_i \right) u_i \\ +& \left( -\frac{2p_i + 2p_{i-1}}{h^4} + \frac{q_{i-\frac{1}{2}}}{h^2} \right) u_{i-1} + \frac{p_{i-1}}{h^4} u_{i-2} = f_i \end{aligned} } \]

5. 边界条件处理(虚拟节点法)

边界条件为 \(u(a)=u'(a)=0\)\(u(b)=u'(b)=0\),采用虚拟节点法处理导数边界条件:

  • 左边界 \(x_0=a\)\(u_0=0\),且 \(u'(a) \approx \frac{u_1 - u_{-1}}{2h} = 0 \implies u_{-1} = u_1\)
  • 右边界 \(x_N=b\)\(u_N=0\),且 \(u'(b) \approx \frac{u_{N+1} - u_{N-1}}{2h} = 0 \implies u_{N+1} = u_{N-1}\)

将虚拟节点 \(u_{-1}\)\(u_{N+1}\) 分别代入 \(i=1\)\(i=N-1\) 的差分方程,消去虚拟节点,得到仅包含内部节点的方程。

6. 系数矩阵性质

  • 该差分格式的系数矩阵是五对角矩阵(每个方程最多涉及5个相邻节点);
  • \(r(x) \geq 0\) 时,系数矩阵是对称正定矩阵,存在唯一解;
  • 整体截断误差为 \(O(h^2)\),是二阶精度的差分格式。

核心知识点总结表

习题 核心方法 关键技巧 最终格式特点
习题1 有限体积法 对偶控制体积+调和平均通量+中心差分对流项 与直接差分法等价,非对称三对角矩阵,二阶精度
习题2 中心差分法 逐层差分法+虚拟节点法处理边界导数 五对角对称正定矩阵,二阶精度,适用于四阶椭圆问题

2.3 矩形网的差分格式 深度讲解

一、二维椭圆型方程模型问题

1. Poisson方程与边值条件

考虑二维Poisson方程:

\[-\Delta u = f(x,y), \quad (x,y) \in G \tag{2.53} \]

其中 \(G\)\(xy\) 平面上的有界区域,边界 \(\Gamma\) 分段光滑。边界条件分为三类:

  • 第一类(Dirichlet):\(u|_\Gamma = \alpha(x,y) \tag{2.54}\)
  • 第二类(Neumann):\(\left.\frac{\partial u}{\partial n}\right|_\Gamma = \beta(x,y) \tag{2.55}\)
  • 第三类(Robin):\(\left.\frac{\partial u}{\partial n} + ku\right|_\Gamma = \gamma(x,y), \ k\geq0 \tag{2.56}\)

2. 二维问题的新挑战

与一维问题相比,二维问题的主要困难在于:

  • 求解域几何形状更复杂,网格剖分难度显著增加;
  • 非规则边界的边界条件处理更复杂,容易引入额外误差;
  • 差分方程组规模随节点数平方增长,对求解算法的效率要求极高。

二、矩形网格剖分与节点类型

1. 网格构造

取x方向步长 \(h_1\),y方向步长 \(h_2\),作两族平行于坐标轴的直线:

\[x = ih_1, \quad i=0,\pm1,\dots \]

\[y = jh_2, \quad j=0,\pm1,\dots \]

两族直线的交点称为网格节点,记为 \((x_i,y_j)\)\((i,j)\)

2. 节点分类

  • 内点:所有属于 \(G\) 内部的节点,集合记为 \(G_h\)
  • 界点:网格线与边界 \(\Gamma\) 的交点,集合记为 \(\Gamma_h\)
  • 全体节点\(\bar{G}_h = G_h \cup \Gamma_h\)
  • 正则内点:四个相邻节点都属于 \(G_h\) 的内点;
  • 非正则内点:至少有一个相邻节点不属于 \(G_h\) 的内点。

关键意义:正则内点可以直接使用标准差分格式,非正则内点需要特殊处理边界条件,是边界误差的主要来源。

三、五点差分格式(核心)

1. 格式构造

推导依据:二阶中心差分公式
对于正则内点 \((i,j)\),分别用x方向和y方向的二阶中心差分近似二阶偏导数:

\[u_{xx}(x_i,y_j) \approx \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h_1^2} \]

\[u_{yy}(x_i,y_j) \approx \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h_2^2} \]

代入Poisson方程 \(-(u_{xx}+u_{yy})=f\),得到:

\[-\Delta_h u_{ij} = -\left( \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h_1^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h_2^2} \right) = f_{ij} \tag{2.57} \]

其中 \(f_{ij}=f(x_i,y_j)\)\(\Delta_h\) 称为五点差分算子

2. 截断误差分析

推导依据:泰勒中值定理
\(u(x_{i+1},y_j)\)\(u(x_{i-1},y_j)\)\((x_i,y_j)\) 处作泰勒展开:

\[u(x_{i+1},y_j) = u_{ij} + h_1 u_x + \frac{h_1^2}{2!}u_{xx} + \frac{h_1^3}{3!}u_{xxx} + \frac{h_1^4}{4!}u_{xxxx} + O(h_1^5) \]

\[u(x_{i-1},y_j) = u_{ij} - h_1 u_x + \frac{h_1^2}{2!}u_{xx} - \frac{h_1^3}{3!}u_{xxx} + \frac{h_1^4}{4!}u_{xxxx} + O(h_1^5) \]

两式相加整理得:

\[\frac{u_{i+1,j} - 2u_{ij} + u_{i-1,j}}{h_1^2} = u_{xx} + \frac{h_1^2}{12}u_{xxxx} + O(h_1^4) \tag{2.59} \]

同理,y方向有:

\[\frac{u_{i,j+1} - 2u_{ij} + u_{i,j-1}}{h_2^2} = u_{yy} + \frac{h_2^2}{12}u_{yyyy} + O(h_2^4) \tag{2.60} \]

将两式相加,得到差分算子与微分算子的差:

\[R_{ij}(u) = \Delta u_{ij} - \Delta_h u_{ij} = -\frac{1}{12}\left( h_1^2 u_{xxxx} + h_2^2 u_{yyyy} \right) + O(h^4) = O(h^2) \tag{2.61} \]

其中 \(h = \sqrt{h_1^2 + h_2^2}\) 为最大步长。

核心结论:五点差分格式的截断误差为 \(O(h^2)\),是二阶精度的差分格式。

3. 正方形网格的简化形式

\(h_1=h_2=h\)(正方形网格)时,差分方程(2.57)简化为:

\[4u_{ij} - u_{i-1,j} - u_{i+1,j} - u_{i,j-1} - u_{i,j+1} = h^2 f_{ij} \]

或写成:

\[u_{ij} - \frac{1}{4}\left( u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1} \right) = \frac{h^2}{4}f_{ij} \tag{2.62} \]

4. Laplace方程的特殊情况

\(f \equiv 0\) 时,方程退化为Laplace方程,差分格式为:

\[u_{ij} = \frac{1}{4}\left( u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1} \right) \tag{2.63} \]

物理意义:Laplace方程的差分解在每个节点上的值等于其四个相邻节点值的算术平均值,这与调和函数的平均值性质完全一致。

5. 系数矩阵性质

将差分方程组按节点顺序排列,得到线性代数方程组 \(A\mathbf{u} = \mathbf{b}\),其系数矩阵 \(A\) 具有以下性质:

  • 稀疏性:每行最多有5个非零元素,稀疏度极高;
  • 对称性\(A^T = A\)
  • 正定性:严格对角占优,故对称正定;
  • 带状性:带宽约为 \(2M+1\)(M为x方向节点数)。

这些性质保证了方程组存在唯一解,且可以使用共轭梯度法、多重网格法等高效算法求解。

四、九点差分格式(高阶精度)

1. 构造思路

五点格式的截断误差为 \(O(h^2)\),为了提高精度,我们可以利用原Poisson方程将截断误差中的四阶导数项转化为右端项 \(f\) 的导数,从而消去低阶误差项,得到四阶精度的格式。

2. 推导过程

由注2.2,将(2.59)和(2.60)相加得:

\[\Delta_h u_{ij} = \Delta u_{ij} + \frac{1}{12}\left( h_1^2 u_{xxxx} + h_2^2 u_{yyyy} \right) + O(h^4) \]

利用原方程 \(\Delta u = -f\),对x求二阶偏导得 \(u_{xxxx} + u_{xxyy} = -f_{xx}\),对y求二阶偏导得 \(u_{yyyy} + u_{xxyy} = -f_{yy}\)。将这两个式子代入上式,并用中心差分近似 \(f_{xx}\)\(f_{yy}\),最终得到九点差分格式:

\[\begin{aligned} -\Delta_h u_{ij} &- \frac{1}{12}\left[ 4u_{ij} - 2(u_{i-1,j} + u_{i,j-1} + u_{i+1,j} + u_{i,j+1}) \right. \\ &\left. + u_{i-1,j-1} + u_{i+1,j-1} + u_{i+1,j+1} + u_{i-1,j+1} \right] \frac{h_1^2 h_2^2}{h_1^2 + h_2^2} \\ &= f_{ij} + \frac{1}{12}\left( h_1^2 f_{xx,ij} + h_2^2 f_{yy,ij} \right) \end{aligned} \]

3. 精度与特点

  • 截断误差为 \(O(h^4)\),是四阶精度的差分格式;
  • 每个方程涉及9个节点,因此称为九点格式;
  • 精度更高,但计算量也更大,适用于对精度要求极高的问题。

五、有限体积法推导五点差分格式

1. 积分守恒形式

推导依据:格林公式
对于任意区域 \(G_{ij}\),由格林公式有:

\[-\iint_{G_{ij}} \Delta u dxdy = -\oint_{\partial G_{ij}} \frac{\partial u}{\partial n} ds = \iint_{G_{ij}} f dxdy \tag{2.64} \]

这是Poisson方程的积分守恒形式,物理意义是流入区域的通量等于区域内的源项。

2. 控制体积构造

采用对偶剖分,以正则内点 \((i,j)\) 为中心构造控制体积(矩形域):

\[G_{ij}: \ x_{i-\frac{1}{2}} \leq x \leq x_{i+\frac{1}{2}}, \ y_{j-\frac{1}{2}} \leq y \leq y_{j+\frac{1}{2}} \]

其中 \(x_{i\pm\frac{1}{2}} = (i\pm\frac{1}{2})h_1\)\(y_{j\pm\frac{1}{2}} = (j\pm\frac{1}{2})h_2\)

3. 边界积分近似

推导依据:积分中矩形公式
控制体积的边界 \(\partial G_{ij}\) 由四条边组成,分别用中矩形公式近似每条边上的法向导数积分:

  • 右边界:\(\int_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}} \frac{\partial u}{\partial x} dy \approx h_2 \cdot \frac{u_{i+1,j} - u_{ij}}{h_1}\)
  • 左边界:\(\int_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}} \frac{\partial u}{\partial x} dy \approx h_2 \cdot \frac{u_{ij} - u_{i-1,j}}{h_1}\)
  • 上边界:\(\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \frac{\partial u}{\partial y} dx \approx h_1 \cdot \frac{u_{i,j+1} - u_{ij}}{h_2}\)
  • 下边界:\(\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \frac{\partial u}{\partial y} dx \approx h_1 \cdot \frac{u_{ij} - u_{i,j-1}}{h_2}\)

4. 差分方程建立

将四条边的积分代入格林公式(2.64),右端项用中矩形公式近似:

\[\iint_{G_{ij}} f dxdy \approx h_1 h_2 f_{ij} \]

整理后得到:

\[-\left( \frac{u_{i+1,j} - 2u_{ij} + u_{i-1,j}}{h_1^2} + \frac{u_{i,j+1} - 2u_{ij} + u_{i,j-1}}{h_2^2} \right) = f_{ij} \tag{2.65} \]

核心结论:有限体积法导出的差分格式与直接差分化法完全一致,但有限体积法具有天然的守恒性,更适合处理物理守恒律问题。

六、核心知识点归纳总结表

核心知识点 关键内容 推导/依据 重要结论
模型问题 二维Poisson方程与三类边值条件 椭圆型方程理论 是静电场、稳态热传导等问题的数学模型
网格剖分 矩形网格、正则内点、非正则内点、界点 区域分解原理 正则内点用标准格式,非正则内点需特殊处理
五点差分格式 \(-\left( \frac{u_{i+1,j}-2u_{ij}+u_{i-1,j}}{h_1^2} + \frac{u_{i,j+1}-2u_{ij}+u_{i,j-1}}{h_2^2} \right) = f_{ij}\) 二阶中心差分公式 最常用的二维椭圆方程差分格式
截断误差 \(O(h^2)\) 泰勒中值定理 二阶精度,满足大多数工程问题需求
正方形网格简化 \(4u_{ij} - u_{i\pm1,j} - u_{i,j\pm1} = h^2 f_{ij}\) 五点格式特例 形式简单,对称性好
Laplace方程格式 \(u_{ij} = \frac{1}{4}(u_{i\pm1,j} + u_{i,j\pm1})\) 五点格式特例 体现调和函数的平均值性质
九点差分格式 涉及9个节点 原方程降阶+中心差分 四阶精度,计算量更大
有限体积法推导 积分守恒形式+格林公式+中矩形公式 格林公式、积分中值定理 与直接差分法等价,天然满足守恒性
系数矩阵性质 稀疏、对称、正定、带状 差分方程结构 存在唯一解,适合高效迭代算法

2.3.2 边值条件的处理 深度讲解

一、第一边值条件(Dirichlet条件)的处理

第一边值条件的形式为:

\[u|_\Gamma = \alpha(x,y) \tag{2.66} \]

处理的核心是区分界点非正则内点,分别采用不同的离散方式。

1. 界点的直接赋值

界点是网格线与边界 \(\Gamma\) 的交点,属于边界 \(\Gamma\) 的一部分,因此可以直接赋值:

\[u_{ij} = \alpha(\bar{x}_i, \bar{y}_j), \quad (\bar{x}_i, \bar{y}_j) \in \Gamma_h \tag{2.67} \]

这一步是精确的,没有截断误差。

2. 非正则内点的不等距差分格式

非正则内点是指至少有一个相邻节点不在求解域内的内点。对于这类节点,不能直接使用标准五点格式,需要构造不等距差分格式

格式推导(以图2.5节点0为例)

推导依据:泰勒中值定理
设节点0的坐标为 \((x_0,y_0)\),其右侧相邻节点1的坐标为 \((x_0+h_1,y_0)\),左侧相邻节点3的坐标为 \((x_0-h_1^-,y_0)\),上下相邻节点2和4的步长均为 \(h_2\)

\(u_1\)\(u_3\) 在节点0处作泰勒展开:

\[u_1 = u_0 + h_1 u_x + \frac{h_1^2}{2} u_{xx} + O(h_1^3) \]

\[u_3 = u_0 - h_1^- u_x + \frac{(h_1^-)^2}{2} u_{xx} + O((h_1^-)^3) \]

消去一阶导数项 \(u_x\),整理得到二阶导数的不等距近似:

\[u_{xx}(x_0,y_0) \approx \frac{1}{\bar{h}_1} \left( \frac{u_1 - u_0}{h_1} - \frac{u_0 - u_3}{h_1^-} \right) \]

其中 \(\bar{h}_1 = \frac{1}{2}(h_1 + h_1^-)\) 是平均步长。

y方向的二阶导数仍用标准中心差分近似,代入Poisson方程得到不等距差分格式:

\[-\left[ \frac{1}{\bar{h}_1} \left( \frac{u_1 - u_0}{h_1} - \frac{u_0 - u_3}{h_1^-} \right) + \frac{1}{h_2^2}(u_2 - 2u_0 + u_4) \right] = f_0 \tag{2.68} \]

截断误差分析

通过泰勒展开余项分析,该格式的截断误差为 \(O(h)\),比正则内点的五点格式低一阶。

致命缺陷:破坏对称正定性

不等距差分格式导出的系数矩阵不再是对称矩阵,这会导致无法使用共轭梯度法等高效的对称矩阵求解算法,显著增加计算量。

3. 保持对称正定性的修正格式

为了保持系数矩阵的对称正定性,我们对不等距格式进行如下修正:

\[-\left[ \frac{1}{h_1} \left( \frac{u_1 - u_0}{h_1} - \frac{u_0 - u_3}{h_1^-} \right) + \frac{1}{h_2^2}(u_2 - 2u_0 + u_4) \right] = f_0 \tag{2.69} \]

关键变化

将平均步长 \(\bar{h}_1\) 替换为右侧步长 \(h_1\),这一微小的改动使得系数矩阵重新恢复对称正定性。

截断误差与收敛阶的矛盾

修正格式的局部截断误差为 \(O(1)\),看起来比不等距格式更差,但可以严格证明,差分解的整体收敛阶仍然是 \(O(h^2)\)(按最大模)

理论依据:非正则内点的数量是 \(O(h^{-1})\) 量级,远少于正则内点的 \(O(h^{-2})\) 量级。虽然每个非正则内点的局部误差是 \(O(1)\),但这些误差对整体误差的贡献是 \(O(h^{-1}) \times O(1) \times h^2 = O(h^2)\),因此整体收敛阶不受影响。

这是数值分析中一个非常重要的结论,它告诉我们:为了保持格式的良好性质(如对称性),可以适当牺牲局部截断误差,只要这些误差的影响是局部的、可控制的。

二、第二、第三边值条件(Neumann/Robin条件)的处理

第二、第三边值条件的统一形式为:

\[\left. \frac{\partial u}{\partial n} + ku \right|_\Gamma = \gamma \tag{2.70} \]

  • \(k=0\) 时,退化为第二边值条件(Neumann条件);
  • \(k \neq 0\) 时,为第三边值条件(Robin条件)。

处理导数边界条件最自然、最有效的方法是有限体积法,它从积分守恒形式出发,将边界条件自然融入差分方程,保持守恒性和对称性。

1. 边界控制体积的构造

对于边界节点 \(P_0(x_{i_0},y_{j_0})\),构造其对应的边界控制体积(对偶单元),即由过 \(P_0\) 相邻半整数点的平行线与边界 \(\Gamma\) 围成的曲边三角形 \(\bar{\triangle}ABC\)(如图2.6所示)。

2. 积分守恒形式的建立

推导依据:格林公式
在边界控制体积 \(\bar{\triangle}ABC\) 上对Poisson方程积分,并应用格林公式:

\[-\iint_{\bar{\triangle}ABC} \Delta u dxdy = -\oint_{\overparen{ABC}A} \frac{\partial u}{\partial n} ds = \iint_{\bar{\triangle}ABC} f dxdy \tag{2.71} \]

其中 \(\frac{\partial u}{\partial n}\) 是外法向导数。

3. 边界积分的近似

将闭合曲线积分分解为三条边的积分之和:

\[\oint_{\overparen{ABC}A} \frac{\partial u}{\partial n} ds = \int_{\overline{AB}} \frac{\partial u}{\partial n} ds + \int_{\overline{BC}} \frac{\partial u}{\partial n} ds + \int_{\overparen{CA}} \frac{\partial u}{\partial n} ds \]

(1) 内部边 \(\overline{AB}\)\(\overline{BC}\) 的积分近似

这两条边位于求解域内部,法向导数用中心差分近似,积分用中矩形公式近似:

\[\int_{\overline{AB}} \frac{\partial u}{\partial n} ds \approx \frac{u_{P_2} - u_{P_0}}{h_2} \cdot \overline{AB} \]

\[\int_{\overline{BC}} \frac{\partial u}{\partial n} ds \approx \frac{u_{P_1} - u_{P_0}}{h_1} \cdot \overline{BC} \]

其中 \(P_1\)\(P_2\) 是与 \(P_0\) 相邻的内部节点,\(\overline{AB}\)\(\overline{BC}\) 是对应边的长度。

(2) 边界边 \(\overparen{CA}\) 的积分近似

这条边位于边界 \(\Gamma\) 上,直接代入边界条件(2.70):

\[\int_{\overparen{CA}} \frac{\partial u}{\partial n} ds = \int_{\overparen{CA}} (\gamma - ku) ds \approx (\gamma_{P_0} - k_{P_0} u_{P_0}) \cdot \overparen{CA} \]

其中 \(\overparen{CA}\) 是边界弧段的长度。

4. 边界差分方程的建立

将三条边的积分近似代入积分守恒方程(2.71),右端项用中矩形公式近似:

\[\iint_{\bar{\triangle}ABC} f dxdy \approx f_{P_0} \cdot S_{\triangle ABC} \]

其中 \(S_{\triangle ABC}\) 是边界控制体积的面积。

整理得到最终的边界差分方程:

\[-\left[ \frac{u_{P_2} - u_{P_0}}{h_2} \cdot \overline{AB} + \frac{u_{P_1} - u_{P_0}}{h_1} \cdot \overline{BC} + (\gamma_{P_0} - k_{P_0} u_{P_0}) \cdot \overparen{CA} \right] = \iint_{\bar{\triangle}ABC} f dxdy \tag{2.72} \]

5. 有限体积法的核心优势

  • 自然融入边界条件:不需要额外构造差分近似,边界条件直接作为积分项的一部分代入;
  • 保持守恒性:在每个控制体积上都严格满足物理守恒律,包括边界控制体积;
  • 保持对称性:导出的边界差分方程与内点方程具有一致的结构,整个系数矩阵仍然保持对称正定性;
  • 易于推广:可以方便地推广到任意形状的边界和复杂的边界条件。

三、核心知识点归纳总结表

边值条件类型 处理方法 核心格式 局部截断误差 对称性 整体收敛阶 优缺点
第一类(Dirichlet) 界点直接赋值 \(u_{ij}=\alpha(\bar{x}_i,\bar{y}_j)\) 0 不影响 \(O(h^2)\) 简单精确,无误差
第一类(Dirichlet) 非正则内点不等距差分 式(2.68) \(O(h)\) 破坏 \(O(h)\) 精度较高,但破坏对称性
第一类(Dirichlet) 非正则内点修正对称格式 式(2.69) \(O(1)\) 保持 \(O(h^2)\) 保持对称性,整体收敛阶不变
第二类(Neumann) 有限体积法 式(2.72) \(O(h)\) 保持 \(O(h^2)\) 自然处理,保持守恒性和对称性
第三类(Robin) 有限体积法 式(2.72) \(O(h)\) 保持 \(O(h^2)\) 自然处理,保持守恒性和对称性

2.3.3 习题 完整解答与推导

(资深高级研究员视角,严格遵循有限体积法与泰勒展开理论)


习题1:变系数扩散方程第一边值问题的有限体积法五点格式

1. 模型问题

考虑二维变系数自伴扩散方程:

\[-\nabla(k\nabla u) = -\left[ \frac{\partial}{\partial x}\left(k\frac{\partial u}{\partial x}\right) + \frac{\partial}{\partial y}\left(k\frac{\partial u}{\partial y}\right) \right] = f(x,y) \tag{2.73} \]

其中 \(k(x,y) \geq k_{\text{min}} > 0\)(保证椭圆性),第一边值条件为 \(u|_\Gamma = \alpha(x,y)\)

2. 有限体积法核心步骤

推导依据:格林公式 + 积分中矩形公式

步骤1:控制体积构造

采用矩形网格,步长 \(h_1\)(x方向)、\(h_2\)(y方向),节点 \((x_i,y_j)=(ih_1,jh_2)\)。以正则内点 \((i,j)\) 为中心构造对偶控制体积:

\[V_{ij}: \ x_{i-\frac{1}{2}} \leq x \leq x_{i+\frac{1}{2}}, \ y_{j-\frac{1}{2}} \leq y \leq y_{j+\frac{1}{2}} \]

其中半整数点 \(x_{i\pm\frac{1}{2}}=(i\pm\frac{1}{2})h_1\)\(y_{j\pm\frac{1}{2}}=(j\pm\frac{1}{2})h_2\),控制体积面积为 \(h_1h_2\)

步骤2:积分守恒形式建立

将方程(2.73)在控制体积 \(V_{ij}\) 上积分,应用格林公式将体积分转化为边界通量积分:

\[-\iint_{V_{ij}} \nabla\cdot(k\nabla u) dxdy = -\oint_{\partial V_{ij}} k\frac{\partial u}{\partial n} ds = \iint_{V_{ij}} f dxdy \]

其中 \(\frac{\partial u}{\partial n}\) 是边界外法向导数。

步骤3:边界通量的积分近似

将闭合边界积分分解为左、右、下、上四条边的积分之和,每条边采用中矩形公式近似(二阶精度):

  • 右边界(\(x=x_{i+\frac{1}{2}}\),外法向+x)

    \[\int_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}} k\frac{\partial u}{\partial x} dy \approx h_2 \cdot k_{i+\frac{1}{2},j} \cdot \frac{u_{i+1,j} - u_{i,j}}{h_1} \]

    其中 \(k_{i+\frac{1}{2},j}=k(x_{i+\frac{1}{2}},y_j)\) 是边中点的系数值(光滑时取中点值,间断时取调和平均)。

  • 左边界(\(x=x_{i-\frac{1}{2}}\),外法向-x)

    \[\int_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}} k\left(-\frac{\partial u}{\partial x}\right) dy \approx h_2 \cdot k_{i-\frac{1}{2},j} \cdot \frac{u_{i-1,j} - u_{i,j}}{h_1} \]

  • 上边界(\(y=y_{j+\frac{1}{2}}\),外法向+y)

    \[\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} k\frac{\partial u}{\partial y} dx \approx h_1 \cdot k_{i,j+\frac{1}{2}} \cdot \frac{u_{i,j+1} - u_{i,j}}{h_2} \]

  • 下边界(\(y=y_{j-\frac{1}{2}}\),外法向-y)

    \[\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} k\left(-\frac{\partial u}{\partial y}\right) dx \approx h_1 \cdot k_{i,j-\frac{1}{2}} \cdot \frac{u_{i,j-1} - u_{i,j}}{h_2} \]

步骤4:差分方程整理

将四条边的通量积分代入守恒方程,右端源项用中矩形公式近似:

\[\iint_{V_{ij}} f dxdy \approx h_1h_2 f_{i,j}, \quad f_{i,j}=f(x_i,y_j) \]

合并同类项并两边除以 \(h_1h_2\),得到最终的五点差分格式:

\[\boxed{ \begin{aligned} &\left( \frac{k_{i+\frac{1}{2},j} + k_{i-\frac{1}{2},j}}{h_1^2} + \frac{k_{i,j+\frac{1}{2}} + k_{i,j-\frac{1}{2}}}{h_2^2} \right) u_{i,j} \\ -&\frac{k_{i+\frac{1}{2},j}}{h_1^2}u_{i+1,j} - \frac{k_{i-\frac{1}{2},j}}{h_1^2}u_{i-1,j} \\ -&\frac{k_{i,j+\frac{1}{2}}}{h_2^2}u_{i,j+1} - \frac{k_{i,j-\frac{1}{2}}}{h_2^2}u_{i,j-1} = f_{i,j} \end{aligned} } \]

3. 边界条件处理

第一边值条件直接对界点赋值:\(u_{ij}=\alpha(\bar{x}_i,\bar{y}_j)\),非正则内点采用修正对称格式处理,保持系数矩阵的对称正定性。


习题2:变系数扩散方程第二边值问题的有限体积法五点格式

1. 模型问题

第二边值条件为:

\[\left. \frac{\partial u}{\partial n} \right|_\Gamma = \beta(x,y) \]

其中 \(\frac{\partial u}{\partial n}\) 是边界外法向导数。

2. 边界控制体积构造

对于边界节点 \(P_0(x_{i_0},y_{j_0})\),构造其对应的边界半控制体积:由过相邻半整数点的平行线与边界 \(\Gamma\) 围成的曲边多边形(通常为曲边三角形或曲边矩形)。

3. 积分守恒形式与边界条件代入

在边界控制体积 \(\bar{V}_0\) 上应用格林公式:

\[-\oint_{\partial \bar{V}_0} k\frac{\partial u}{\partial n} ds = \iint_{\bar{V}_0} f dxdy \]

将闭合边界积分分解为内部边积分边界边积分两部分:

  • 内部边积分:与内点格式完全相同,采用中心差分近似通量;
  • 边界边积分:直接代入第二边值条件 \(\frac{\partial u}{\partial n}=\beta\),积分近似为:

    \[\int_{\overparen{CA}} k\frac{\partial u}{\partial n} ds \approx k_{P_0} \beta_{P_0} \cdot \overparen{CA} \]

    其中 \(\overparen{CA}\) 是边界弧段的长度,\(k_{P_0}=k(P_0)\)\(\beta_{P_0}=\beta(P_0)\)

4. 边界差分方程

以图2.6所示的角点边界节点为例,最终得到边界差分方程的一般形式:

\[\boxed{ \begin{aligned} &\left( \frac{k_{i_0+\frac{1}{2},j_0} \cdot \overline{BC}}{h_1} + \frac{k_{i_0,j_0-\frac{1}{2}} \cdot \overline{AB}}{h_2} + k_{P_0} \cdot \overparen{CA} \right) u_{P_0} \\ -&\frac{k_{i_0+\frac{1}{2},j_0} \cdot \overline{BC}}{h_1} u_{P_1} - \frac{k_{i_0,j_0-\frac{1}{2}} \cdot \overline{AB}}{h_2} u_{P_2} = k_{P_0} \beta_{P_0} \cdot \overparen{CA} + \iint_{\bar{V}_0} f dxdy \end{aligned} } \]

其中 \(P_1,P_2\) 是与 \(P_0\) 相邻的内部节点,\(\overline{AB},\overline{BC}\) 是内部边的长度。

5. 核心优势

有限体积法处理第二、第三边值条件时,自然将边界条件融入守恒方程,保持了格式的守恒性和系数矩阵的对称正定性,避免了简单差分法的精度损失和对称性破坏。


习题3:高阶差分格式的截断误差证明

1. 问题重述

定义菱形差分算子和方形差分算子:

\[\Diamond u_{ij} = u_{i+1,j} + u_{i,j+1} + u_{i-1,j} + u_{i,j-1} - 4u_{ij} \]

\[\square u_{ij} = u_{i+1,j+1} + u_{i-1,j+1} + u_{i-1,j-1} + u_{i+1,j-1} - 4u_{ij} \]

证明逼近Laplace方程 \(\Delta u=0\) 的差分格式:

\[\frac{1}{6h^2}(4\Diamond u_{ij} + \square u_{ij}) = 0 \]

的截断误差满足:

\[|R_{i,j}(u)| = \frac{40h^6}{3\cdot 8!} \theta M_8, \quad |\theta| \leq 1 \]

其中 \(M_8\)\(u\) 的所有八阶偏导数的绝对值在求解域上的上确界。

2. 证明过程

推导依据:二元泰勒展开定理 + Laplace方程解的偏导数性质

步骤1:邻点函数值的泰勒展开

将所有邻点的函数值在 \((x_i,y_j)\) 处展开到八阶项:

  • 轴邻点(上下左右):

    \[u(x\pm h,y) = u \pm h u_x + \frac{h^2}{2!}u_{xx} \pm \frac{h^3}{3!}u_{xxx} + \dots + \frac{h^8}{8!}u_{xxxxxxxx} + O(h^9) \]

    相加得:

    \[u(x+h,y)+u(x-h,y) = 2u + h^2 u_{xx} + \frac{h^4}{12}u_{xxxx} + \frac{h^6}{360}u_{xxxxxx} + \frac{h^8}{20160}u_{xxxxxxxx} + O(h^9) \]

    同理可得y方向轴邻点的和。

  • 对角邻点(四个角点):
    利用二元泰勒展开,四个对角点相加后奇次项全部抵消,得到:

    \[\begin{aligned} u(x+h,y+h)+u(x-h,y+h)+u(x-h,y-h)+u(x+h,y-h) &= 4u + 2h^2(u_{xx}+u_{yy}) \\ + \frac{h^4}{6}(u_{xxxx}+6u_{xxyy}+u_{yyyy}) &+ \frac{h^6}{180}(u_{xxxxxx}+15u_{xxxxyy}+15u_{xxyyyy}+u_{yyyyyy}) \\ + \frac{h^8}{10080}(u_{xxxxxxxx}+28u_{xxxxxxyy}+70u_{xxxxyyyy}+28u_{xxyyyyyy}+u_{yyyyyyyy}) &+ O(h^9) \end{aligned} \]

步骤2:差分算子的展开

  • 菱形算子:

    \[\Diamond u_{ij} = h^2\Delta u + \frac{h^4}{12}(u_{xxxx}+u_{yyyy}) + \frac{h^6}{360}(u_{xxxxxx}+u_{yyyyyy}) + \frac{h^8}{20160}(u_{xxxxxxxx}+u_{yyyyyyyy}) + O(h^9) \]

  • 方形算子:

    \[\begin{aligned} \square u_{ij} &= 2h^2\Delta u + \frac{h^4}{6}(u_{xxxx}+6u_{xxyy}+u_{yyyy}) \\ + \frac{h^6}{180}(u_{xxxxxx}+15u_{xxxxyy}+15u_{xxyyyy}+u_{yyyyyy}) &+ \frac{h^8}{10080}(u_{xxxxxxxx}+28u_{xxxxxxyy}+70u_{xxxxyyyy}+28u_{xxyyyyyy}+u_{yyyyyyyy}) + O(h^9) \end{aligned} \]

步骤3:合并差分算子并化简

计算 \(4\Diamond u_{ij} + \square u_{ij}\),合并同类项:

\[\begin{aligned} 4\Diamond u_{ij} + \square u_{ij} &= 6h^2\Delta u + \frac{h^4}{2}\Delta^2 u + \frac{h^6}{60}\Delta^3 u + \frac{h^6}{30}(u_{xxxxyy}+u_{xxyyyy}) \\ + \frac{h^8}{3360}\Delta^4 u + \frac{h^8}{630}(u_{xxxxxxyy}+u_{xxyyyyyy}) &+ \frac{13h^8}{2520}u_{xxxxyyyy} + O(h^9) \end{aligned} \]

其中 \(\Delta^2 u = u_{xxxx}+2u_{xxyy}+u_{yyyy}\) 是双调和算子,\(\Delta^3 u,\Delta^4 u\) 是更高阶的调和算子。

步骤4:利用Laplace方程条件消去低阶项

对于Laplace方程的解 \(u\),有 \(\Delta u=0\),因此所有调和算子的幂次均为0:\(\Delta^2 u=\Delta^3 u=\Delta^4 u=0\)

进一步利用 \(u_{yy}=-u_{xx}\) 化简混合偏导数:

  • \(u_{xxxxyy} + u_{xxyyyy} = -u_{xxxxxx} + u_{xxxxxx} = 0\)(h⁴项抵消)
  • \(u_{xxxxxxyy} + u_{xxyyyyyy} = -2u_{xxxxxxxx}\)
  • \(u_{xxxxyyyy} = u_{xxxxxxxx}\)

代入后得到:

\[4\Diamond u_{ij} + \square u_{ij} = -\frac{h^8}{3024}u_{xxxxxxxx} + O(h^9) \]

步骤5:截断误差计算

差分算子为:

\[\frac{4\Diamond u_{ij} + \square u_{ij}}{6h^2} = -\frac{h^6}{18144}u_{xxxxxxxx} + O(h^7) \]

截断误差定义为:

\[R_{ij}(u) = \Delta u - \frac{4\Diamond u_{ij} + \square u_{ij}}{6h^2} = \frac{h^6}{18144}u_{xxxxxxxx} + O(h^7) \]

注意到 \(\frac{1}{18144} = \frac{40}{3\cdot 40320} = \frac{40}{3\cdot 8!}\),且 \(|u_{xxxxxxxx}| \leq M_8\),因此:

\[|R_{ij}(u)| = \frac{40h^6}{3\cdot 8!} |\theta| M_8 \leq \frac{40h^6}{3\cdot 8!} M_8, \quad |\theta| \leq 1 \]

证毕

3. 结论

该差分格式是六阶精度的高阶格式,比标准五点格式(二阶精度)精度高得多,适用于对精度要求极高的Laplace方程求解问题。


核心知识点总结表

习题 核心方法 关键技巧 格式特点
习题1 有限体积法 对偶控制体积+边中点系数+中矩形积分 对称正定五对角矩阵,二阶精度,守恒性好
习题2 有限体积法 边界半控制体积+边界条件直接代入 保持守恒性与对称性,自然处理导数边界条件
习题3 泰勒展开法 高阶展开+Laplace方程偏导数化简 六阶精度,计算量略大,适合高精度需求

2.4 三角网的差分格式(广义差分法)深度讲解

一、三角网差分格式的本质与核心优势

1. 方法定位

三角网差分格式是有限体积法在非结构网格上的推广,文献中也称为广义差分法。它将矩形网的差分思想扩展到任意形状的三角网格,彻底解决了矩形网难以处理复杂几何区域的问题。

2. 核心优势

  • 网格灵活性极强:可以逼近任意复杂的求解域边界,是工程实际中处理不规则区域的首选方法;
  • 边界条件处理自然:特别是第二、第三类导数边界条件,处理方式与内点完全一致,无需额外构造近似;
  • 保持积分守恒性:在每个控制体积上严格满足物理守恒律,与矩形网有限体积法具有相同的守恒性质;
  • 支持局部自适应加密:可以在解变化剧烈的区域局部细化网格,在保证精度的同时显著减少计算量。

3. 模型问题

考虑有界域G上的Poisson方程:

\[-\Delta u = f(x,y), \quad (x,y) \in G \tag{2.74} \]

边界\(\Gamma\)上满足第一、第二或第三边值条件。

二、三角剖分的基本要求

1. 三角剖分的定义

  1. 在边界\(\Gamma\)上取一系列节点,作逼近\(\Gamma\)的闭折线\(\tilde{\Gamma}\),得到逼近原区域G的多边形域\(\tilde{G}\)
  2. \(\tilde{G}\)分割为有限个互不重叠的小三角形(三角单元);
  3. 任一三角形的顶点不能是其他三角形边的内部点。

2. 关键几何约束(必须满足)

所有三角形的内角不大于90°(即非钝角三角形)。

几何意义与必要性

  • 保证三角形的外心位于三角形内部或边上;
  • 保证围绕每个节点的对偶单元是凸多边形;
  • 保证差分方程组的系数矩阵是对称正定矩阵,解存在唯一且稳定。

三、对偶剖分(控制体积)的构造

对偶剖分是三角网有限体积法的核心,它为每个节点构造一个专属的控制体积,使得守恒律可以在每个控制体积上独立成立。

1. 构造步骤

对于任意节点\(P_0\)

  1. 找出所有以\(P_0\)为顶点的三角单元;
  2. 对每个三角单元的每条边作中垂线
  3. 同一三角单元的两条边的中垂线交于该三角形的外心
  4. 依次连接这些外心,得到围绕\(P_0\)的闭合多边形,即为\(P_0\)对偶单元(控制体积),记为\(K_{P_0}^*\)

2. 对偶单元的类型

  • 内点的对偶单元:如图2.8(a)所示,是一个凸多边形(通常为六边形),\(P_0\)位于多边形内部;
  • 界点的对偶单元:如图2.8(b)所示,是一个半多边形,\(P_0\)是多边形的一个顶点,其两条边位于求解域边界上。

四、内点差分方程的完整推导

推导依据:格林公式 + 积分中矩形公式 + 线性插值近似

步骤1:积分守恒形式建立

在节点\(P_0\)的对偶单元\(K_{P_0}^*\)上积分Poisson方程(2.74):

\[-\iint_{K_{P_0}^*} \Delta u dxdy = \iint_{K_{P_0}^*} f dxdy \]

步骤2:格林公式转化为通量积分

应用格林公式,将体积分转化为对偶单元边界上的法向导数积分:

\[-\oint_{\partial K_{P_0}^*} \frac{\partial u}{\partial n} ds = \iint_{K_{P_0}^*} f dxdy \tag{2.75} \]

其中\(\partial K_{P_0}^*\)是对偶单元的边界,\(n\)是单位外法向量。

步骤3:边界积分的分段近似

对偶单元的边界由若干条线段组成,每条线段对应一个相邻三角单元的外心连线。设\(P_0\)的相邻节点为\(P_1,P_2,\dots,P_6\)\(P_7=P_1\)),\(q_i\)是三角形\(P_0 P_i P_{i+1}\)的外心,则对偶单元的边界为\(q_1 q_2, q_2 q_3, \dots, q_6 q_1\)

关键几何性质:线段\(q_i q_{i+1}\)是三角形\(P_0 P_i P_{i+1}\)的外心连线,垂直于边\(P_0 P_{i+1}\)。因此,在\(q_i q_{i+1}\)上,法向导数\(\frac{\partial u}{\partial n}\)就是沿\(P_0 P_{i+1}\)方向的方向导数。

用两点差商近似方向导数,积分用中矩形公式近似,得到:

\[\int_{\overline{q_i q_{i+1}}} \frac{\partial u}{\partial n} ds \approx \frac{u(P_{i+1}) - u(P_0)}{\overline{P_0 P_{i+1}}} \cdot \overline{q_i q_{i+1}} \tag{2.76} \]

步骤4:差分方程整理

将所有边的积分代入式(2.75),右端源项用中矩形公式近似:

\[\iint_{K_{P_0}^*} f dxdy \approx m(K_{P_0}^*) \cdot \varphi_0, \quad \varphi_0 = \frac{1}{m(K_{P_0}^*)} \iint_{K_{P_0}^*} f dxdy \]

其中\(m(K_{P_0}^*)\)是对偶单元的面积。

最终得到内点\(P_0\)的差分方程:

\[-\sum_{i=1}^6 \frac{\overline{q_i q_{i+1}}}{\overline{P_0 P_{i+1}}} (u(P_{i+1}) - u(P_0)) = m(K_{P_0}^*) \cdot \varphi_0 \tag{2.77} \]

五、界点差分方程的处理

1. 第一边值条件(Dirichlet)

直接对界点赋值:

\[u(P_0) = \alpha(P_0) \tag{2.78} \]

与矩形网处理方式完全一致。

2. 第二、第三边值条件(Neumann/Robin)

这是三角网差分格式的最大优势,处理方式与内点几乎完全相同。

边界条件的统一形式为:

\[\left. \frac{\partial u}{\partial n} + ku \right|_\Gamma = \gamma \tag{2.79} \]

推导步骤

  1. 在界点\(P_0\)的边界对偶单元上应用格林公式,得到与式(2.75)类似的积分守恒方程;
  2. 对偶单元的边界分为内部边边界边两部分:
    • 内部边的积分近似与内点完全相同;
    • 边界边的积分直接代入边界条件(2.79),消去法向导数;
  3. 假设\(u\)在边界边上是线性函数,用梯形公式近似积分:

    \[\int_{\overline{m_4 P_0}} \frac{\partial u}{\partial n} ds = \int_{\overline{m_4 P_0}} (\gamma - ku) ds \approx \frac{1}{2} \overline{P_0 P_4} \left[ \gamma - \frac{1}{4}k(3u(P_0) + u(P_4)) \right] \tag{2.81} \]

最终边界差分方程

将所有积分近似代入,得到界点的差分方程,其形式与内点方程完全一致,仅系数包含边界条件的参数。

六、系数矩阵的性质

所有内点和界点的差分方程组成一个封闭的线性代数方程组\(A\mathbf{u} = \mathbf{b}\),其系数矩阵\(A\)具有以下优良性质:

  1. 稀疏性:每行最多包含与\(P_0\)相邻的节点对应的非零元素,稀疏度极高;
  2. 对称性:对于任意两个相邻节点\(P_i\)\(P_j\),系数\(a_{ij}=a_{ji}\),因此\(A^T=A\)
  3. 正定性:当三角剖分满足所有内角不大于90°时,矩阵\(A\)严格对角占优,因此对称正定,方程组存在唯一解。

七、三角网与矩形网差分格式的全面对比

对比维度 矩形网差分格式 三角网差分格式(广义差分法)
网格类型 结构网格 非结构网格
几何适应性 差,仅适合规则矩形区域 极强,适合任意复杂区域
边界条件处理 导数边界条件处理复杂,易损失精度 导数边界条件处理自然,与内点一致
局部加密 困难,易产生悬挂节点 容易,支持自适应局部加密
守恒性 保持 保持
系数矩阵性质 对称正定 对称正定
编程实现 简单 稍复杂,需要网格生成与数据结构支持
适用场景 规则区域、大规模并行计算 复杂几何区域、工程实际问题

八、核心知识点归纳总结表

核心知识点 关键内容 推导/依据 重要结论
方法本质 有限体积法在非结构网格上的推广 积分守恒律 又称广义差分法
三角剖分要求 所有三角形内角不大于90° 几何性质 保证外心在内部、对偶单元凸、矩阵正定
对偶单元构造 边中垂线交于外心,连接外心形成控制体积 中垂线性质 每个节点对应一个专属控制体积
内点差分方程 式(2.77) 格林公式+中矩形公式 形式简洁,物理意义明确
边界条件处理 第一类直接赋值,第二类第三类自然融入 边界条件代入积分 处理方式与内点一致,无精度损失
系数矩阵 稀疏、对称、正定 差分方程结构 存在唯一解,适合高效迭代算法
核心优势 网格灵活、边界处理自然、守恒性好 非结构网格特性 是工程实际中处理复杂区域的首选方法

例2.1 直角三角剖分导出五点差分格式 深度解析

一、例子的核心意义

本例子是连接矩形网差分法三角网差分法(广义差分法)的关键桥梁,它严格证明了:当三角剖分是由矩形网通过同向对角线分割得到时,三角网有限体积法导出的差分格式与标准矩形网五点差分格式完全等价

这一结论具有深远的理论意义:

  1. 验证了三角网差分格式的正确性,它在特殊情况下退化为已经被广泛验证的经典格式;
  2. 证明了有限体积法是一个统一的数值方法框架,涵盖了结构网格和非结构网格的所有情况;
  3. 为理解两种方法的内在联系提供了直观的几何解释。

二、直角三角剖分的构造

1. 基础矩形网

首先建立标准矩形网格,x方向步长\(h_1\),y方向步长\(h_2\),节点坐标为\((x_i,y_j)=(ih_1,jh_2)\),每个矩形单元为\([x_i,x_{i+1}] \times [y_j,y_{j+1}]\)

2. 三角剖分

同向对角线(例如从左下角到右上角的对角线)将每个矩形单元分割为两个全等的直角三角形。这样得到的三角剖分满足以下性质:

  • 所有三角形都是直角三角形,内角分别为\(90^\circ\)\(\arctan(h_2/h_1)\)\(\arctan(h_1/h_2)\),均不大于\(90^\circ\),满足三角剖分的核心几何约束;
  • 每个内部节点周围有6个相邻三角形(上下左右各一个,对角线方向两个);
  • 所有三角形的边要么平行于坐标轴,要么是矩形的对角线。

三、对偶单元(控制体积)的构造

1. 直角三角形外心的关键性质

几何定理:直角三角形的外心位于斜边的中点

这是本例子推导的核心几何依据。对于每个由矩形对角线分割得到的直角三角形,其斜边就是原矩形的对角线,因此外心就是矩形对角线的中点,坐标为\((x_i+h_1/2, y_j+h_2/2)\)

2. 对偶单元的形成

对于内部节点\(P_0(x_i,y_j)\),找出所有以\(P_0\)为顶点的6个三角形,它们的外心分别是:

  • 上方矩形对角线中点:\((x_i+h_1/2, y_j+h_2/2)\)
  • 左上方矩形对角线中点:\((x_i-h_1/2, y_j+h_2/2)\)
  • 左下方矩形对角线中点:\((x_i-h_1/2, y_j-h_2/2)\)
  • 下方矩形对角线中点:\((x_i+h_1/2, y_j-h_2/2)\)

依次连接这四个外心,得到的闭合多边形是一个矩形,其顶点坐标为:

\[(x_i\pm h_1/2, y_j\pm h_2/2) \]

这个矩形正好就是矩形网有限体积法中以\(P_0\)为中心的对偶控制体积,面积为\(h_1 h_2\)

四、差分方程的完整推导

1. 三角网差分公式回顾

根据三角网内点差分方程(2.77):

\[-\sum_{k} \frac{\overline{q_k q_{k+1}}}{\overline{P_0 P_{k+1}}} (u(P_{k+1}) - u(P_0)) = m(K_{P_0}^*) \cdot \varphi_0 \]

2. 各项系数计算

对于节点\(P_0(x_i,y_j)\),其相邻节点有四个:

  • \(P_1(x_{i+1},y_j)\)(右邻点)
  • \(P_2(x_i,y_{j+1})\)(上邻点)
  • \(P_3(x_{i-1},y_j)\)(左邻点)
  • \(P_4(x_i,y_{j-1})\)(下邻点)

分别计算每个相邻节点对应的系数:

  • 右邻点\(P_1\):外心连线长度\(\overline{q_k q_{k+1}} = h_2\),节点距离\(\overline{P_0 P_1} = h_1\),系数为\(\frac{h_2}{h_1}\)
  • 上邻点\(P_2\):外心连线长度\(\overline{q_k q_{k+1}} = h_1\),节点距离\(\overline{P_0 P_2} = h_2\),系数为\(\frac{h_1}{h_2}\)
  • 左邻点\(P_3\):外心连线长度\(\overline{q_k q_{k+1}} = h_2\),节点距离\(\overline{P_0 P_3} = h_1\),系数为\(\frac{h_2}{h_1}\)
  • 下邻点\(P_4\):外心连线长度\(\overline{q_k q_{k+1}} = h_1\),节点距离\(\overline{P_0 P_4} = h_2\),系数为\(\frac{h_1}{h_2}\)

3. 差分方程整理

将系数代入差分公式,右端项为对偶单元面积乘以源项平均值:

\[m(K_{P_0}^*) \cdot \varphi_0 = h_1 h_2 \cdot f_{ij} \]

得到:

\[\begin{aligned} &-\left[ \frac{h_2}{h_1}(u_{i+1,j} - u_{ij}) + \frac{h_1}{h_2}(u_{i,j+1} - u_{ij}) + \frac{h_2}{h_1}(u_{i-1,j} - u_{ij}) + \frac{h_1}{h_2}(u_{i,j-1} - u_{ij}) \right] \\ =& h_1 h_2 f_{ij} \end{aligned} \]

合并同类项:

\[-\left[ \frac{h_2}{h_1}(u_{i+1,j} + u_{i-1,j}) + \frac{h_1}{h_2}(u_{i,j+1} + u_{i,j-1}) - 2\left( \frac{h_2}{h_1} + \frac{h_1}{h_2} \right) u_{ij} \right] = h_1 h_2 f_{ij} \]

两边同时除以\(h_1 h_2\),最终得到:

\[-\left( \frac{u_{i+1,j} - 2u_{ij} + u_{i-1,j}}{h_1^2} + \frac{u_{i,j+1} - 2u_{ij} + u_{i,j-1}}{h_2^2} \right) = f_{ij} \]

这正是标准的矩形网五点差分格式!

五、结论与拓展

1. 核心结论

通过同向对角线分割矩形网得到的直角三角剖分,其三角网差分格式与矩形网五点差分格式完全等价

2. 拓展讨论

  • 如果采用反向对角线(从左上角到右下角)分割矩形单元,推导过程完全相同,最终也会得到相同的五点差分格式;
  • 这一等价性不仅适用于Poisson方程,也适用于变系数扩散方程等其他自伴椭圆型方程;
  • 当三角剖分不是由矩形网分割得到时,三角网差分格式会退化为更一般的非结构网格格式,但仍然保持守恒性和对称正定性。

六、核心知识点归纳总结表

步骤 关键内容 几何依据 结论
三角剖分 矩形网同向对角线分割为直角三角形 直角三角形内角均不大于90° 满足三角剖分的几何约束
外心位置 直角三角形外心在斜边中点 直角三角形外心定理 外心为矩形对角线中点
对偶单元 连接外心形成矩形控制体积 中垂线性质 与矩形网对偶单元完全一致
系数计算 外心连线长度与节点距离的比值 几何距离计算 系数与矩形网格式系数匹配
方程整理 代入三角网差分公式化简 代数运算 得到标准五点差分格式

例2.2 正三角网上的差分格式 深度解析

一、例子的核心意义

正三角网是最规则的非结构网格,具有完美的各向同性对称性。本例子通过推导正三角网的差分格式,展示了三角网有限体积法在规则非结构网格下的简洁形式,同时验证了其与矩形网五点格式具有相同的二阶精度,是理解非结构网格数值方法的经典范例。

二、正三角剖分的几何性质

1. 基本几何特征

正三角网由边长均为\(h\)的等边三角形组成,具有以下关键性质:

  • 所有三角形内角均为\(60^\circ\),严格满足"所有内角不大于\(90^\circ\)"的三角剖分约束;
  • 等边三角形的外心、重心、垂心、内心四心合一,这是正三角网最特殊的几何性质;
  • 每个内部节点周围均匀分布着6个相邻节点,形成完美的六边形对称结构。

2. 关键几何量的严格推导

(1) 节点间距

任意两个相邻节点之间的距离均为\(h\),即:

\[\overline{P_0 P_i} = h, \quad i=1,2,\dots,6 \]

(2) 外心连线长度

等边三角形的外接圆半径(外心到顶点的距离)为:

\[R = \frac{h}{\sqrt{3}} \]

相邻两个三角形的外心连线平行于三角形的公共边,且长度等于外接圆半径,因此:

\[\overline{q_i q_{i+1}} = \frac{h}{\sqrt{3}} \]

(3) 对偶单元(正六边形)的面积

围绕节点\(P_0\)的对偶单元是一个边长为\(\frac{h}{\sqrt{3}}\)的正六边形。正六边形可以分解为6个全等的等边三角形,每个小等边三角形的面积为:

\[S_{\text{小三角}} = \frac{\sqrt{3}}{4} \cdot \left( \frac{h}{\sqrt{3}} \right)^2 = \frac{\sqrt{3} h^2}{12} \]

因此正六边形的总面积为:

\[m(G_0) = 6 \cdot \frac{\sqrt{3} h^2}{12} = \frac{\sqrt{3} h^2}{2} \]

三、差分格式的完整推导

1. 三角网内点差分公式回顾

根据之前推导的三角网内点差分方程(2.77):

\[-\sum_{i=1}^6 \frac{\overline{q_i q_{i+1}}}{\overline{P_0 P_{i+1}}} (u(P_{i+1}) - u(P_0)) = \iint_{G_0} f(x,y) dxdy \]

2. 代入正三角网的几何量

由于正三角网的完美对称性,所有相邻节点对应的系数完全相同:

\[\frac{\overline{q_i q_{i+1}}}{\overline{P_0 P_{i+1}}} = \frac{h/\sqrt{3}}{h} = \frac{1}{\sqrt{3}}, \quad i=1,2,\dots,6 \]

代入差分公式,左边化简为:

\[-\sum_{i=1}^6 \frac{1}{\sqrt{3}} (u_i - u_0) = -\frac{1}{\sqrt{3}} \left( \sum_{i=1}^6 u_i - 6u_0 \right) \]

3. 方程整理与标准化

将左右两边同时乘以\(\frac{2}{\sqrt{3} h^2}\),消去分母:
左边:

\[-\frac{1}{\sqrt{3}} \left( \sum_{i=1}^6 u_i - 6u_0 \right) \cdot \frac{2}{\sqrt{3} h^2} = -\frac{2}{3h^2} \left( \sum_{i=1}^6 u_i - 6u_0 \right) \]

右边:

\[\iint_{G_0} f(x,y) dxdy \cdot \frac{2}{\sqrt{3} h^2} = \frac{2}{\sqrt{3} h^2} \iint_{G_0} f(x,y) dxdy \]

最终得到正三角网的差分格式:

\[\boxed{ -\frac{2}{3h^2} \left( \sum_{i=1}^6 u_i - 6u_0 \right) = \frac{2}{\sqrt{3} h^2} \iint_{G_0} f(x,y) dxdy \tag{2.83} } \]

4. 源项的简化近似

当源项\(f(x,y)\)光滑时,用中矩形公式近似积分:

\[\iint_{G_0} f(x,y) dxdy \approx m(G_0) \cdot f_0 = \frac{\sqrt{3} h^2}{2} f_0 \]

代入式(2.83),得到更简洁的实用形式:

\[6u_0 - \sum_{i=1}^6 u_i = \frac{3h^2}{2} f_0 \]

四、正三角网差分格式的核心特性

1. 完美的对称性

六个相邻节点的权重完全相同,差分格式具有\(60^\circ\)旋转不变性,这是矩形网五点格式不具备的优势。对于各向同性的物理问题(如稳态热传导、静电场),正三角网的数值解具有更好的各向同性精度。

2. 二阶精度

通过泰勒展开可以严格证明,正三角网差分格式的截断误差为\(O(h^2)\),与矩形网五点格式的精度相同。

3. 守恒性与对称正定性

与所有有限体积法格式一样,正三角网差分格式严格保持积分守恒性,且系数矩阵是对称正定矩阵,保证了解的存在唯一性和稳定性。

五、正三角网与矩形网差分格式的全面对比

对比维度 矩形网五点格式 正三角网六点格式
相邻节点数 4个 6个
对称性 90°旋转不变 60°旋转不变,各向同性更好
精度 二阶 二阶
守恒性 保持 保持
系数矩阵 对称正定 对称正定
几何适应性 差,仅适合矩形区域 较好,可逼近一般多边形区域
编程实现 简单 稍复杂
各向同性问题精度 一般 更优

六、核心知识点归纳总结表

核心知识点 关键内容 推导/依据 重要结论
正三角剖分 等边三角形组成,内角60° 几何性质 满足三角剖分约束,四心合一
对偶单元 正六边形,面积\(\frac{\sqrt{3}h^2}{2}\) 正六边形面积公式 完美对称的控制体积
差分格式 \(6u_0 - \sum_{i=1}^6 u_i = \frac{3h^2}{2}f_0\) 三角网差分公式+几何量代入 形式简洁,对称性好
格式特性 二阶精度、各向同性、守恒、对称正定 泰勒展开+有限体积法性质 适合各向同性物理问题
对比矩形网 各向同性更好,几何适应性更强 对称性分析 非结构网格的经典代表

例2.3 正六边形网上的差分格式 深度解析

一、例子的核心意义

正六边形网(又称蜂窝网格)是自然界最优化的网格结构,具有完美的各向同性和最高的空间填充效率。本例子展示了有限体积法的极致通用性:无论原网格是三角形、矩形还是多边形,都可以通过对偶剖分构造出守恒型差分格式。

正六边形网的差分格式是三点格式,比矩形网的五点格式、正三角网的六点格式更简洁,计算效率更高,同时保持了二阶精度和完美的对称性,在蜂窝通信、图像处理、地理信息系统等领域有广泛应用。

二、正六边形剖分的几何性质

1. 基本几何特征

正六边形网由边长均为\(h\)的正六边形组成,具有以下关键性质:

  • 所有正六边形内角均为\(120^\circ\),无钝角,满足有限体积法的几何约束;
  • 每个内部节点周围均匀分布着3个相邻节点,形成\(120^\circ\)旋转对称结构;
  • 空间填充效率最高:相同面积下,正六边形网的总边长最短,节点数最少。

2. 关键几何量的严格推导

(1) 节点间距

任意两个相邻节点之间的距离等于正六边形的边长:

\[\overline{P_0 P_i} = h, \quad i=1,2,3 \]

(2) 对偶单元的构造

核心对偶性原理:正六边形网与正三角网互为对偶网格。

  • 过每条边\(\overline{P_0 P_i}\)的中点作中垂线;
  • 三条中垂线分别交于相邻三个正六边形的中心\(q_1,q_2,q_3\)
  • 这三个中心构成的正三角形\(\triangle q_1 q_2 q_3\)就是围绕\(P_0\)的对偶单元(控制体积)。

(3) 对偶单元的边长

正六边形的中心到其边中点的距离为\(\frac{\sqrt{3}}{2}h\),因此两个相邻正六边形中心之间的距离为:

\[\overline{q_i q_{i+1}} = 2 \times \frac{\sqrt{3}}{2}h = \sqrt{3}h \]

这也是对偶正三角形的边长。

(4) 对偶单元的面积

正三角形面积公式为\(S = \frac{\sqrt{3}}{4}a^2\),代入边长\(a=\sqrt{3}h\),得:

\[m(\triangle q_1 q_2 q_3) = \frac{\sqrt{3}}{4} \times (\sqrt{3}h)^2 = \frac{3\sqrt{3}h^2}{4} \]

三、差分格式的完整推导

1. 有限体积法通用差分公式

对于任意网格的内点\(P_0\),有限体积法的差分方程统一形式为:

\[-\sum_{i=1}^n \frac{\overline{q_i q_{i+1}}}{\overline{P_0 P_{i+1}}} (u(P_{i+1}) - u(P_0)) = \iint_{K_{P_0}^*} f(x,y) dxdy \]

其中\(n\)为相邻节点数,\(\overline{q_i q_{i+1}}\)为对偶单元的边长,\(\overline{P_0 P_{i+1}}\)为节点间距。

2. 代入正六边形网的几何量

由于正六边形网的完美对称性,三个相邻节点对应的系数完全相同:

\[\frac{\overline{q_i q_{i+1}}}{\overline{P_0 P_{i+1}}} = \frac{\sqrt{3}h}{h} = \sqrt{3}, \quad i=1,2,3 \]

代入通用公式,左边化简为:

\[-\sum_{i=1}^3 \sqrt{3} (u_i - u_0) = -\sqrt{3} \left( \sum_{i=1}^3 u_i - 3u_0 \right) \]

3. 方程标准化

将左右两边同时乘以\(\frac{4}{3\sqrt{3}h^2}\),消去分母:
左边:

\[-\sqrt{3} \left( \sum_{i=1}^3 u_i - 3u_0 \right) \times \frac{4}{3\sqrt{3}h^2} = -\frac{4}{3h^2} \left( \sum_{i=1}^3 u_i - 3u_0 \right) \]

右边:

\[\iint_{\triangle q_1 q_2 q_3} f(x,y) dxdy \times \frac{4}{3\sqrt{3}h^2} = \frac{4}{3\sqrt{3}h^2} \iint_{\triangle q_1 q_2 q_3} f(x,y) dxdy \]

最终得到正六边形网的差分格式:

\[\boxed{ -\frac{4}{3h^2} \left( \sum_{i=1}^3 u_i - 3u_0 \right) = \frac{4}{3\sqrt{3}h^2} \iint_{\triangle q_1 q_2 q_3} f(x,y) dxdy \tag{2.84} } \]

4. 源项的简化实用形式

当源项\(f(x,y)\)光滑时,用中矩形公式近似积分:

\[\iint_{\triangle q_1 q_2 q_3} f(x,y) dxdy \approx m(\triangle q_1 q_2 q_3) \cdot f_0 = \frac{3\sqrt{3}h^2}{4} f_0 \]

代入式(2.84),得到工程上最常用的简化形式:

\[3u_0 - (u_1 + u_2 + u_3) = \frac{3h^2}{4} f_0 \]

四、正六边形网差分格式的核心优势

1. 极致的各向同性

差分格式具有\(120^\circ\)旋转不变性,是所有规则网格中各向同性最好的。对于各向同性的物理问题(如稳态热传导、静电场),正六边形网的数值解完全没有网格取向效应,精度显著高于矩形网。

2. 更高的计算效率

每个节点仅涉及3个相邻节点,系数矩阵的稀疏度比矩形网(4个邻点)和正三角网(6个邻点)更高,计算量更小,内存占用更低。

3. 优良的理论性质

  • 截断误差为\(O(h^2)\),与矩形网、正三角网精度相同;
  • 严格保持积分守恒性;
  • 系数矩阵对称正定,解存在唯一且稳定。

五、三种规则网格差分格式的全面对比

对比维度 矩形网五点格式 正三角网六点格式 正六边形网三点格式
原网格单元 正方形 等边三角形 正六边形
对偶网格单元 正方形 正六边形 正三角形
相邻节点数 4 6 3
旋转对称性 \(90^\circ\) \(60^\circ\) \(120^\circ\)
各向同性 一般 较好 最好
空间填充效率 最高
计算效率 一般 较低 最高
差分格式形式 \(4u_0 - \sum u_i = h^2 f_0\) \(6u_0 - \sum u_i = \frac{3h^2}{2}f_0\) \(3u_0 - \sum u_i = \frac{3h^2}{4}f_0\)
适用场景 规则矩形区域、大规模并行计算 一般多边形区域、各向同性问题 蜂窝通信、图像处理、高精度各向同性计算

六、核心知识点归纳总结表

核心知识点 关键内容 推导/依据 重要结论
正六边形剖分 边长为h的正六边形,内角120° 几何性质 空间填充效率最高,各向同性最好
网格对偶性 正六边形网与正三角网互为对偶 对偶剖分定义 两者的差分格式具有内在联系
对偶单元 正三角形,面积\(\frac{3\sqrt{3}h^2}{4}\) 正三角形面积公式 控制体积形状规则,便于积分近似
差分格式 \(3u_0 - (u_1+u_2+u_3) = \frac{3h^2}{4}f_0\) 有限体积法通用公式+几何量代入 形式最简洁,计算效率最高
格式特性 二阶精度、完美各向同性、守恒、对称正定 泰勒展开+有限体积法性质 无网格取向效应,精度最优
工程应用 蜂窝通信、图像处理、地理信息系统 各向同性优势 是未来高精度数值模拟的重要发展方向

2.4.1 习题 完整解答与严格证明

(资深高级研究员视角,遵循有限体积法理论体系,补充关键推导细节)


习题1:三角网差分格式系数矩阵对称性证明

证明依据:有限体积法通量的双向性 + 边界条件的对称处理

1. 内点差分方程的一般形式

对于任意内点\(P_i\),其差分方程由对偶单元上的积分守恒律导出:

\[-\sum_{j \in \text{adj}(i)} \frac{l_{ij}}{d_{ij}} (u_j - u_i) = \iint_{K_i^*} f dxdy \]

其中:

  • \(\text{adj}(i)\)表示\(P_i\)的所有相邻节点集合;
  • \(d_{ij} = \overline{P_i P_j}\)是节点\(P_i\)\(P_j\)之间的距离;
  • \(l_{ij}\)是对偶单元上对应于边\(\overline{P_i P_j}\)的边的长度。

将方程整理为标准线性形式:

\[a_{ii} u_i + \sum_{j \in \text{adj}(i)} a_{ij} u_j = b_i \]

其中:

  • 自系数:\(a_{ii} = \sum_{j \in \text{adj}(i)} \frac{l_{ij}}{d_{ij}}\)
  • 互系数:\(a_{ij} = -\frac{l_{ij}}{d_{ij}} \quad (j \neq i)\)

2. 内点系数对称性证明

对于任意两个相邻节点\(P_i\)\(P_j\)

  • \(\overline{P_i P_j}\)\(P_i\)\(P_j\)的公共边,因此\(d_{ij} = d_{ji}\)
  • 对偶单元上对应于\(\overline{P_i P_j}\)的边是同一条线段,因此\(l_{ij} = l_{ji}\)

由此可得:

\[a_{ij} = -\frac{l_{ij}}{d_{ij}} = -\frac{l_{ji}}{d_{ji}} = a_{ji} \]

即任意两个内点之间的互系数相等。

3. 边界条件的对称性保持

(1) 第一边值条件(Dirichlet)

第一边值条件为\(u_k = \alpha_k\)\(P_k\)为界点)。处理方式是将界点的已知值直接代入相邻内点的差分方程,消去界点未知量。这一过程仅相当于对系数矩阵的行和列进行线性组合,不会破坏剩余内点矩阵的对称性。

(2) 第三边值条件(Robin)

第三边值条件为\(\left. \frac{\partial u}{\partial n} + ku \right|_\Gamma = \gamma\)。采用有限体积法处理时,边界控制体积的积分方程为:

\[-\sum_{j \in \text{adj}(k)} \frac{l_{kj}}{d_{kj}} (u_j - u_k) + \left( \int_{\Gamma_k} k ds \right) u_k = \iint_{K_k^*} f dxdy + \int_{\Gamma_k} \gamma ds \]

其中\(\Gamma_k\)是边界控制体积的边界边。该方程的自系数为正,且与相邻内点的互系数仍然满足\(a_{kj}=a_{jk}\),因此整体系数矩阵保持对称。

结论:三角网差分格式(第一或第三边值条件)的系数矩阵是对称矩阵。


习题2:变系数扩散方程的三角网差分格式构造

1. 模型问题回顾

变系数自伴扩散方程(2.73)为:

\[-\nabla(k\nabla u) = -\left[ \frac{\partial}{\partial x}\left(k\frac{\partial u}{\partial x}\right) + \frac{\partial}{\partial y}\left(k\frac{\partial u}{\partial y}\right) \right] = f(x,y) \]

其中\(k(x,y) \geq k_{\text{min}} > 0\)(保证椭圆性)。

2. 有限体积法推导步骤

推导依据:格林公式 + 积分中矩形公式

步骤1:对偶单元积分

在节点\(P_i\)的对偶单元\(K_i^*\)上积分方程,应用格林公式将体积分转化为边界通量积分:

\[-\iint_{K_i^*} \nabla\cdot(k\nabla u) dxdy = -\oint_{\partial K_i^*} k\frac{\partial u}{\partial n} ds = \iint_{K_i^*} f dxdy \]

步骤2:边界通量的变系数近似

对于对偶单元的每条边\(\overline{q_a q_b}\)(对应于原网格的边\(\overline{P_i P_j}\)),通量积分采用中矩形公式近似:

\[\int_{\overline{q_a q_b}} k\frac{\partial u}{\partial n} ds \approx k_{ij} \cdot \frac{u_j - u_i}{d_{ij}} \cdot l_{ij} \]

其中:

  • \(k_{ij} = k\left( \frac{P_i + P_j}{2} \right)\)是边\(\overline{P_i P_j}\)中点的系数值(光滑时取中点值,间断时取调和平均);
  • \(d_{ij} = \overline{P_i P_j}\)是节点间距;
  • \(l_{ij} = \overline{q_a q_b}\)是对偶边长度。

步骤3:差分方程整理

将所有边的通量积分代入守恒方程,右端源项用中矩形公式近似:

\[\iint_{K_i^*} f dxdy \approx m(K_i^*) \cdot f_i, \quad f_i = f(P_i) \]

最终得到变系数扩散方程的三角网差分格式:

\[\boxed{ -\sum_{j \in \text{adj}(i)} \frac{k_{ij} l_{ij}}{d_{ij}} (u_j - u_i) = \iint_{K_i^*} f(x,y) dxdy } \]

3. 特殊情况验证

\(k(x,y) \equiv 1\)(常系数)时,格式退化为Poisson方程的三角网差分格式(2.77),验证了其正确性。


习题3:差分格式截断误差阶的计算

1. 正三角网差分格式(2.83)的截断误差(\(O(h^2)\)

推导依据:二元泰勒展开 + 正三角网对称性

步骤1:差分算子与微分算子

正三角网差分格式(2.83)的差分算子为:

\[L_h u_0 = -\frac{2}{3h^2} \left( \sum_{i=1}^6 u_i - 6u_0 \right) \]

对应的微分算子为\(L u = -\Delta u\)

步骤2:邻点坐标与泰勒展开

\(P_0\)位于坐标原点\((0,0)\),六个相邻节点的坐标为:

\[(h,0), \left( \frac{h}{2}, \frac{h\sqrt{3}}{2} \right), \left( -\frac{h}{2}, \frac{h\sqrt{3}}{2} \right), (-h,0), \left( -\frac{h}{2}, -\frac{h\sqrt{3}}{2} \right), \left( \frac{h}{2}, -\frac{h\sqrt{3}}{2} \right) \]

将每个邻点的函数值在\(P_0\)处展开到四阶项:

\[u(x,y) = u_0 + x u_x + y u_y + \frac{x^2}{2}u_{xx} + xy u_{xy} + \frac{y^2}{2}u_{yy} + O(h^3) \]

步骤3:求和与化简

将六个邻点的展开式相加,利用正三角网的对称性,所有奇次项(\(x,y,x^3,x^2y\)等)全部抵消,剩余偶次项:

\[\sum_{i=1}^6 u_i = 6u_0 + \frac{3h^2}{2}(u_{xx} + u_{yy}) + O(h^4) = 6u_0 + \frac{3h^2}{2}\Delta u_0 + O(h^4) \]

步骤4:截断误差计算

代入差分算子:

\[L_h u_0 = -\frac{2}{3h^2} \left( \frac{3h^2}{2}\Delta u_0 + O(h^4) \right) = -\Delta u_0 + O(h^2) \]

因此,正三角网差分格式的截断误差为\(O(h^2)\),是二阶精度格式。


2. 正六边形网差分格式(2.84)的截断误差(\(O(h)\)

推导依据:二元泰勒展开 + 正六边形网对称性

步骤1:差分算子与微分算子

正六边形网差分格式(2.84)的差分算子为:

\[L_h u_0 = -\frac{4}{3h^2} \left( \sum_{i=1}^3 u_i - 3u_0 \right) \]

对应的微分算子为\(L u = -\Delta u\)

步骤2:邻点坐标与泰勒展开

\(P_0\)位于坐标原点\((0,0)\),三个相邻节点的坐标为:

\[(h,0), \left( -\frac{h}{2}, \frac{h\sqrt{3}}{2} \right), \left( -\frac{h}{2}, -\frac{h\sqrt{3}}{2} \right) \]

将每个邻点的函数值在\(P_0\)处展开到三阶项:

\[u(x,y) = u_0 + x u_x + y u_y + \frac{x^2}{2}u_{xx} + xy u_{xy} + \frac{y^2}{2}u_{yy} + \frac{x^3}{6}u_{xxx} + \frac{x^2 y}{2}u_{xxy} + \frac{x y^2}{2}u_{xyy} + \frac{y^3}{6}u_{yyy} + O(h^4) \]

步骤3:求和与化简

将三个邻点的展开式相加,一阶项(\(u_x,u_y\))和混合二阶项(\(u_{xy}\))抵消,但三阶奇次项无法完全抵消

\[\begin{aligned} \sum_{i=1}^3 u_i &= 3u_0 + \frac{3h^2}{4}(u_{xx} + u_{yy}) + \frac{h^3}{12}u_{xxx} - \frac{3h^3}{8}u_{xyy} + O(h^4) \\ &= 3u_0 + \frac{3h^2}{4}\Delta u_0 + O(h^3) \end{aligned} \]

步骤4:截断误差计算

代入差分算子:

\[\begin{aligned} L_h u_0 &= -\frac{4}{3h^2} \left( \frac{3h^2}{4}\Delta u_0 + O(h^3) \right) \\ &= -\Delta u_0 + O(h) \end{aligned} \]

因此,正六边形网差分格式的截断误差为\(O(h)\),是一阶精度格式。


核心知识点总结表

习题 核心结论 关键依据 重要意义
习题1 三角网差分格式系数矩阵对称 通量双向性 + 边界对称处理 保证可使用共轭梯度法等高效对称求解器
习题2 变系数扩散方程三角网格式:\(-\sum \frac{k_{ij} l_{ij}}{d_{ij}}(u_j-u_i) = \int f dxdy\) 格林公式 + 中矩形通量近似 可处理非均匀介质扩散问题
习题3 正三角网格式二阶精度,正六边形网格式一阶精度 二元泰勒展开 + 网格对称性 解释了不同规则网格的精度差异

*2.5 极值定理和敛速估计(第一部分:差分方程)深度讲解

一、本节的核心地位与理论框架

1. 研究目的

前面我们已经构造了各种差分格式,但一个关键问题尚未解决:差分解是否收敛到精确解?收敛的速度有多快? 为了回答这些问题,我们需要对差分解进行先验估计,而极值定理是进行这类估计的最基本、最有力的工具。

2. 理论逻辑链条

  1. 构造一般二阶椭圆型方程的差分格式;
  2. 分析差分方程系数的正定性与对角占优性;
  3. 基于系数性质证明差分算子的极值原理;
  4. 利用极值原理建立差分解的先验估计;
  5. 由先验估计导出收敛性与敛速估计。

二、一般二阶椭圆型方程模型

1. 连续问题

考虑最一般的二阶线性椭圆型偏微分方程第一边值问题:

\[\begin{cases} -\left(A u_x\right)_x - \left(B u_y\right)_y + C u_x + D u_y + E u = F, & (x,y) \in G \tag{2.85} \\ u|_\Gamma = \alpha \end{cases} \]

其中系数满足以下椭圆性与光滑性条件:

  • \(A(x,y) \in C^1(\bar{G}), \ B(x,y) \in C^1(\bar{G})\),且 \(A(x,y) \geq A_{\text{min}} > 0, \ B(x,y) \geq B_{\text{min}} > 0\)(保证方程椭圆性);
  • \(C(x,y), D(x,y), E(x,y), F(x,y) \in C(\bar{G})\),且 \(E(x,y) \geq 0\)(保证解的唯一性);
  • \(\alpha(x,y) \in C(\Gamma)\)

重要特例

  • \(A=B=1, C=D=E=0\) 时,退化为Poisson方程 \(-\Delta u = F\)
  • \(C \neq 0\)\(D \neq 0\) 时,方程包含对流项,称为对流扩散方程。

2. 网格假设

采用矩形网格剖分,\(x\) 方向步长 \(h_1\)\(y\) 方向步长 \(h_2\)。假设网格内点集合 \(G_h\)连通的,即任意两个内点之间都存在一条由相邻内点组成的路径。这是极值定理成立的必要条件。

三、差分格式的构造

1. 差分算子定义

为了逼近方程中的各阶导数,我们定义以下差分算子:

\[\begin{cases} (u_{ij})_{\bar{x}} = \frac{u_{ij} - u_{i-1,j}}{h_1}, & (u_{ij})_{\bar{y}} = \frac{u_{ij} - u_{i,j-1}}{h_2} \quad \text{(向后差分)} \\ (u_{ij})_x = \frac{u_{i+1,j} - u_{ij}}{h_1}, & (u_{ij})_y = \frac{u_{i,j+1} - u_{ij}}{h_2} \quad \text{(向前差分)} \\ (u_{ij})_{\hat{x}} = \frac{u_{i+1,j} - u_{i-1,j}}{2h_1}, & (u_{ij})_{\hat{y}} = \frac{u_{i,j+1} - u_{i,j-1}}{2h_2} \quad \text{(中心差分)} \end{cases} \tag{2.87} \]

2. 正则内点的差分方程

对于正则内点 \((x_i,y_j)\),采用有限体积法构造差分格式:

  • 二阶扩散项 \(-(A u_x)_x\) 用半整数点向后差分近似:\(-\left(A_{i-\frac{1}{2},j} (u_{ij})_{\bar{x}}\right)_x\)
  • 二阶扩散项 \(-(B u_y)_y\) 用半整数点向后差分近似:\(-\left(B_{i,j-\frac{1}{2}} (u_{ij})_{\bar{y}}\right)_y\)
  • 一阶对流项 \(C u_x\)\(D u_y\) 用中心差分近似:\(C_{ij} (u_{ij})_{\hat{x}} + D_{ij} (u_{ij})_{\hat{y}}\)
  • 零阶项 \(E u\) 和右端项 \(F\) 直接在节点处取值。

最终得到正则内点的差分方程:

\[-\left(A_{i-\frac{1}{2},j} (u_{ij})_{\bar{x}}\right)_x - \left(B_{i,j-\frac{1}{2}} (u_{ij})_{\bar{y}}\right)_y + C_{ij} (u_{ij})_{\hat{x}} + D_{ij} (u_{ij})_{\hat{y}} + E_{ij} u_{ij} = F_{ij} \tag{2.86} \]

截断误差分析:通过泰勒展开可以证明,该格式的截断误差为 \(O(h_1^2 + h_2^2)\),是二阶精度格式。

3. 差分方程的标准五点形式

将差分算子展开并整理,得到标准的五点差分格式:

\[a_{ij} u_{ij} - \left(a_{i-1,j} u_{i-1,j} + a_{i,j-1} u_{i,j-1} + a_{i+1,j} u_{i+1,j} + a_{i,j+1} u_{i,j+1}\right) = F_{ij} \tag{2.88} \]

其中系数定义为:

\[\begin{cases} a_{i-1,j} = h_1^{-2} \left(A_{i-\frac{1}{2},j} + \frac{h_1}{2} C_{ij}\right) \\ a_{i,j-1} = h_2^{-2} \left(B_{i,j-\frac{1}{2}} + \frac{h_2}{2} D_{ij}\right) \\ a_{i+1,j} = h_1^{-2} \left(A_{i+\frac{1}{2},j} - \frac{h_1}{2} C_{ij}\right) \\ a_{i,j+1} = h_2^{-2} \left(B_{i,j+\frac{1}{2}} - \frac{h_2}{2} D_{ij}\right) \\ a_{ij} = h_1^{-2} \left(A_{i+\frac{1}{2},j} + A_{i-\frac{1}{2},j}\right) + h_2^{-2} \left(B_{i,j+\frac{1}{2}} + B_{i,j-\frac{1}{2}}\right) + E_{ij} \end{cases} \tag{2.89} \]

4. 非正则内点的差分方程

对于非正则内点,采用不等距差分格式近似导数。可以证明,只要步长 \(h_1, h_2\) 充分小,非正则内点的差分方程仍然可以写成与(2.88)相同的形式,且系数保持以下关键性质:

  • 所有系数 \(a_{i-1,j}, a_{i,j-1}, a_{i+1,j}, a_{i,j+1}, a_{ij}\) 均为正数;
  • 满足对角占优条件:\(a_{ij} - (a_{i-1,j} + a_{i,j-1} + a_{i+1,j} + a_{i,j+1}) = E_{ij} \geq 0\)

注意:非正则内点的截断误差为 \(O(h_1 + h_2)\),是一阶精度,但如前所述,这不会影响差分解的整体二阶收敛性。

四、差分方程系数的核心性质

1. 系数正定性

定理:当步长 \(h_1, h_2\) 充分小时,差分方程(2.88)的所有系数均为正数。

证明

  • 由椭圆性条件,\(A_{i\pm\frac{1}{2},j} \geq A_{\text{min}} > 0\)\(B_{i,j\pm\frac{1}{2}} \geq B_{\text{min}} > 0\)
  • 由于 \(C_{ij}, D_{ij}\) 是连续函数,在有界闭域 \(\bar{G}\) 上有界,即存在常数 \(M\) 使得 \(|C_{ij}| \leq M\)\(|D_{ij}| \leq M\)
  • \(h_1 < \frac{2A_{\text{min}}}{M}\) 时,\(A_{i+\frac{1}{2},j} - \frac{h_1}{2} C_{ij} > 0\),故 \(a_{i+1,j} > 0\)
  • 同理,当 \(h_2 < \frac{2B_{\text{min}}}{M}\) 时,\(a_{i,j+1} > 0\)
  • 其余系数显然为正。

证毕

2. 对角占优性

定理:差分方程(2.88)的系数满足严格对角占优条件:

\[a_{ij} \geq a_{i-1,j} + a_{i,j-1} + a_{i+1,j} + a_{i,j+1} \tag{2.92} \]

当且仅当 \(E_{ij} = 0\) 时等号成立。

证明
由系数定义(2.89)直接计算:

\[\begin{aligned} a_{ij} - (a_{i-1,j} + a_{i,j-1} + a_{i+1,j} + a_{i,j+1}) &= \left[ h_1^{-2}(A_{i+\frac{1}{2},j} + A_{i-\frac{1}{2},j}) + h_2^{-2}(B_{i,j+\frac{1}{2}} + B_{i,j-\frac{1}{2}}) + E_{ij} \right] \\ &- \left[ h_1^{-2}(A_{i-\frac{1}{2},j} + \frac{h_1}{2}C_{ij}) + h_2^{-2}(B_{i,j-\frac{1}{2}} + \frac{h_2}{2}D_{ij}) \right. \\ &\left. + h_1^{-2}(A_{i+\frac{1}{2},j} - \frac{h_1}{2}C_{ij}) + h_2^{-2}(B_{i,j+\frac{1}{2}} - \frac{h_2}{2}D_{ij}) \right] \\ &= E_{ij} \geq 0 \end{aligned} \]

证毕

五、系数矩阵的性质

将差分方程组按节点顺序排列,得到线性代数方程组 \(A\mathbf{u} = \mathbf{b}\),其系数矩阵 \(A\) 具有以下三个关键性质:

1. 稀疏性

矩阵 \(A\) 的每行最多有5个非零元素,因此是高度稀疏的。在排列节点顺序时,应尽量使非零元素靠近对角线,这对于直接消元法(如LU分解)特别有利,可以显著减少计算量和内存占用。

2. 对角占优性

矩阵 \(A\)严格对角占优矩阵(当 \(E \not\equiv 0\) 时)或弱对角占优矩阵(当 \(E \equiv 0\) 时)。这一性质具有极其重要的意义:

  • 保证了线性方程组存在唯一解;
  • 保证了雅可比迭代、高斯-赛德尔迭代等基本迭代法的收敛性;
  • 是证明差分算子极值原理的基础。

3. 对称性

当原微分方程对称,即 \(C = D = 0\) 时,差分方程的系数矩阵 \(A\) 是对称矩阵。此时可以使用共轭梯度法等更高效的迭代算法,计算量比非对称矩阵小一个数量级以上。

六、核心知识点归纳总结表

核心知识点 关键内容 推导/依据 重要意义
模型问题 一般二阶线性椭圆型方程 \(-(Au_x)_x-(Bu_y)_y+Cu_x+Du_y+Eu=F\) 椭圆型方程理论 涵盖了Poisson方程、对流扩散方程等常见模型
差分格式 有限体积法构造的五点格式 差分算子定义+积分守恒律 二阶精度,适用于一般椭圆型方程
系数正定性 所有系数均为正数 椭圆性条件+步长充分小 是极值定理成立的前提
对角占优性 \(a_{ij} \geq \sum a_{kl}\) 系数定义直接计算 保证解的唯一性和迭代法收敛性
系数矩阵性质 稀疏、对角占优、对称(当C=D=0时) 差分方程结构 决定了求解算法的选择和效率
非正则内点处理 不等距差分格式,保持系数性质 泰勒展开 保证整体格式的一致性和稳定性

2.5.2 极值定理 深度讲解

一、极值定理的核心地位

差分方程的极值定理是连续椭圆型方程极值原理的离散版本,是整个差分法理论体系的基石。它的核心作用是:

  1. 证明差分方程解的存在唯一性
  2. 建立差分解的先验估计(最大模估计);
  3. 推导差分解的收敛性敛速估计
  4. 分析差分格式的稳定性

所有后续的理论分析都直接或间接依赖于极值定理,而极值定理的成立又完全依赖于上一节证明的差分方程系数的正定性对角占优性

二、定理2.2(极值定理)

1. 定理完整表述

定理2.2(极值定理)\(u_{ij}\)\(\bar{G}_h\) 上的任一网格函数。若 \(L_h u_{ij} \leq 0\)\(L_h u_{ij} \geq 0\))对任意 \((x_i,y_j) \in G_h\) 成立,则 \(u_{ij}\) 不可能在内点取正的极大值(负的极小值),除非 \(u_{ij} \equiv\) 常数。

2. 核心解读

  • 该定理是连续椭圆型方程强极值原理的离散对应;
  • 它表明:差分方程的解的最大值和最小值只能在边界上取得,除非解是常数;
  • 这一性质是椭圆型方程最本质的特征,它反映了稳态物理过程的平衡特性。

3. 严格证明(反证法)

证明第一部分(\(L_h u_{ij} \leq 0\) 时不能取正极大值)

  1. 反设:假设 \(u_{ij} \not\equiv\) 常数,且在 \(G_h\) 中某点达到正的极大值 \(M > 0\)
  2. 连通性应用:由于 \(G_h\) 是连通的,必然存在一个内点 \((x_{i_0}, y_{j_0})\),使得 \(u_{i_0 j_0} = M\),且至少有一个相邻网点(例如 \((x_{i_0-1}, y_{j_0})\))满足 \(u_{i_0-1,j_0} < M\)
  3. 差分算子展开:将差分算子 \(L_h\) 作用在 \((x_{i_0}, y_{j_0})\) 处:

    \[L_h u_{i_0 j_0} = a_{i_0 j_0} u_{i_0 j_0} - a_{i_0-1,j_0} u_{i_0-1,j_0} - a_{i_0,j_0-1} u_{i_0,j_0-1} - a_{i_0+1,j_0} u_{i_0+1,j_0} - a_{i_0,j_0+1} u_{i_0,j_0+1} \]

  4. 不等式估计
    • 由于 \(u_{i_0 j_0} = M\) 是极大值,故 \(u_{i_0,j_0\pm1} \leq M\)\(u_{i_0\pm1,j_0} \leq M\)
    • 由假设,\(u_{i_0-1,j_0} < M\),且所有系数 \(a_{kl} > 0\)
    • 因此:

      \[\begin{aligned} L_h u_{i_0 j_0} &> a_{i_0 j_0} M - a_{i_0-1,j_0} M - a_{i_0,j_0-1} M - a_{i_0+1,j_0} M - a_{i_0,j_0+1} M \\ &= \left( a_{i_0 j_0} - (a_{i_0-1,j_0} + a_{i_0,j_0-1} + a_{i_0+1,j_0} + a_{i_0,j_0+1}) \right) M \\ &= E_{i_0 j_0} M \geq 0 \end{aligned} \]

  5. 矛盾:这与定理条件 \(L_h u_{ij} \leq 0\) 矛盾。因此假设不成立,\(u_{ij}\) 不能在内点取正的极大值。

第二部分(\(L_h u_{ij} \geq 0\) 时不能取负极小值) 的证明完全类似,只需令 \(v_{ij} = -u_{ij}\),则 \(L_h v_{ij} = -L_h u_{ij} \leq 0\),应用第一部分结论即可。

证毕

三、推论2.1(解的存在唯一性)

1. 推论表述

推论2.1 差分方程(2.91)有唯一解。

2. 证明

根据线性代数理论,非齐次线性方程组有唯一解的充要条件是对应的齐次方程组只有零解。

考虑齐次差分方程:

\[\begin{cases} L_h u_{ij} = 0, & (x_i,y_j) \in G_h \\ u_{ij} = 0, & (x_i,y_j) \in \Gamma_h \end{cases} \]

由极值定理:

  • 因为 \(L_h u_{ij} = 0 \leq 0\),所以 \(u_{ij}\) 不能在内点取正的极大值,故 \(u_{ij} \leq 0\)
  • 因为 \(L_h u_{ij} = 0 \geq 0\),所以 \(u_{ij}\) 不能在内点取负的极小值,故 \(u_{ij} \geq 0\)

因此 \(u_{ij} \equiv 0\),齐次方程只有零解,非齐次方程有唯一解。

证毕

四、推论2.2(非负性原理)

1. 推论表述

推论2.2 若网格函数 \(u_{ij}\) 满足:

\[\begin{cases} L_h u_{ij} \geq 0, & \forall (x_i,y_j) \in G_h \\ u_{ij} \geq 0, & \forall (x_i,y_j) \in \Gamma_h \end{cases} \]

\(u_{ij} \geq 0, \ \forall (x_i,y_j) \in \bar{G}_h\)

2. 证明

反设存在内点 \((x_i,y_j)\) 使得 \(u_{ij} < 0\),则 \(u_{ij}\) 在内点取负的极小值。但由条件 \(L_h u_{ij} \geq 0\),根据极值定理,\(u_{ij}\) 不能在内点取负的极小值,矛盾。因此 \(u_{ij} \geq 0\) 对所有点成立。

证毕

3. 核心意义

非负性原理是比较定理的基础,它告诉我们:如果差分算子作用在函数上非负,且边界值非负,那么整个函数在整个区域上都非负。这是一个非常强大的定性性质。

五、定理2.3(比较定理)

1. 定理完整表述

定理2.3(比较定理)\(u_{ij}\)\(U_{ij}\) 是两个网格函数,满足:

\[\begin{cases} |L_h u_{ij}| \leq L_h U_{ij}, & \forall (x_i,y_j) \in G_h \\ |u_{ij}| \leq U_{ij}, & \forall (x_i,y_j) \in \Gamma_h \end{cases} \]

\[|u_{ij}| \leq U_{ij}, \quad \forall (x_i,y_j) \in \bar{G}_h \tag{2.95} \]

2. 证明

考虑两个辅助函数:

\[v_{ij} = U_{ij} - u_{ij}, \quad w_{ij} = U_{ij} + u_{ij} \]

对于 \(v_{ij}\)

  • 在内点:\(L_h v_{ij} = L_h U_{ij} - L_h u_{ij} \geq L_h U_{ij} - |L_h u_{ij}| \geq 0\)
  • 在边界:\(v_{ij} = U_{ij} - u_{ij} \geq U_{ij} - |u_{ij}| \geq 0\)

由推论2.2(非负性原理),得 \(v_{ij} \geq 0\),即 \(U_{ij} \geq u_{ij}\)

对于 \(w_{ij}\)

  • 在内点:\(L_h w_{ij} = L_h U_{ij} + L_h u_{ij} \geq L_h U_{ij} - |L_h u_{ij}| \geq 0\)
  • 在边界:\(w_{ij} = U_{ij} + u_{ij} \geq U_{ij} - |u_{ij}| \geq 0\)

同理,由推论2.2得 \(w_{ij} \geq 0\),即 \(U_{ij} \geq -u_{ij}\)

综上,\(|u_{ij}| \leq U_{ij}\) 对所有点成立。

证毕

3. 核心意义

比较定理是极值定理最重要的应用,它提供了一种通用的误差估计方法

  • 要估计未知函数 \(u_{ij}\) 的大小,只需构造一个已知的“比较函数” \(U_{ij}\)
  • 只要控制住 \(U_{ij}\) 的差分算子和边界值,就能控制住 \(u_{ij}\) 的最大模;
  • 这是建立差分解先验估计和收敛性证明的核心工具。

六、推论2.3(最大模估计)

1. 推论表述

推论2.3 差分方程

\[\begin{cases} L_h u_{ij} = 0, & (x_i,y_j) \in G_h \\ u_{ij} = \alpha_{ij}, & (x_i,y_j) \in \Gamma_h \end{cases} \]

的解 \(u_{ij}\) 满足不等式:

\[\max_{\bar{G}_h} |u_{ij}| \leq \max_{\Gamma_h} |\alpha_{ij}| \tag{2.96} \]

2. 证明

构造比较函数 \(U_{ij} = \max_{\Gamma_h} |\alpha_{ij}|\)(常数函数)。

验证比较定理的条件:

  • 在内点:\(L_h U_{ij} = E_{ij} U_{ij} \geq 0 = |L_h u_{ij}|\)
  • 在边界:\(U_{ij} = \max_{\Gamma_h} |\alpha_{ij}| \geq |\alpha_{ij}| = |u_{ij}|\)

由比较定理,得 \(|u_{ij}| \leq U_{ij} = \max_{\Gamma_h} |\alpha_{ij}|\) 对所有点成立。

证毕

3. 核心意义

  • 该推论给出了齐次差分方程解的最大模估计,它表明解的最大模不超过边界值的最大模;
  • 这是一个非常简洁且强大的估计,在实际计算中可以用来验证数值解的合理性;
  • 对于非齐次方程,我们可以通过构造更复杂的比较函数,得到类似的最大模估计。

七、核心知识点归纳总结表

定理/推论 核心结论 证明依据 主要应用
定理2.2(极值定理) 解的极值只能在边界上取得,除非是常数 反证法 + 系数正定性 + 对角占优性 所有后续理论的基础
推论2.1(存在唯一性) 差分方程有唯一解 齐次方程只有零解 保证差分格式的适定性
推论2.2(非负性原理) 算子非负且边界非负,则函数非负 极值定理 比较定理的基础
定理2.3(比较定理) 控制算子和边界值,就能控制函数 非负性原理 + 辅助函数构造 误差估计、收敛性证明
推论2.3(最大模估计) 齐次方程解的最大模不超过边界值最大模 比较定理 + 常数比较函数 解的定性分析、数值验证

2.5.3 五点差分格式的敛速估计 + 2.5.4 习题 完整解答


2.5.3 五点差分格式的敛速估计

一、误差方程的建立

\(u=u(x,y)\) 是Poisson方程第一边值问题的精确解,\(u_{ij}\) 是五点差分格式的数值解,定义误差函数

\[e_{ij} = u(x_i,y_j) - u_{ij} \]

将精确解代入差分方程,得到误差满足的方程:

\[\begin{cases} L_h e_{ij} = R_{ij}, & (x_i,y_j) \in G_h \\ e_{ij} = 0, & (x_i,y_j) \in \Gamma_h \end{cases} \tag{2.97} \]

其中 \(R_{ij}\) 是截断误差,其阶为:

\[R_{ij} = \begin{cases} O(h_1^2 + h_2^2), & \text{正则内点} \\ O(h_1 + h_2), & \text{非正则内点} \end{cases} \]

\(h = \sqrt{h_1^2 + h_2^2}\) 为最大步长,则存在与 \(h\) 无关的常数 \(K\),使得对所有内点有:

\[|R_{ij}| \leq Kh \]

二、比较函数的构造

为了估计误差 \(e_{ij}\) 的最大模,我们构造如下二次比较函数

\[E_{ij} = \frac{Kh}{4} \left( R^2 - x_i^2 - y_j^2 \right) \]

其中 \(R\) 是以原点为圆心且包含整个求解域 \(G\) 的最小圆的半径。

比较函数的性质验证

  1. 非负性:对所有 \((x_i,y_j) \in \bar{G}_h\),有 \(x_i^2 + y_j^2 \leq R^2\),故 \(E_{ij} \geq 0\)
  2. 差分算子作用结果
    • 对于正则内点,五点差分算子作用于二次函数:

      \[L_h x_i^2 = -\frac{(x_{i+1}^2 - 2x_i^2 + x_{i-1}^2)}{h^2} = -\frac{( (i+1)^2 h^2 - 2i^2 h^2 + (i-1)^2 h^2 )}{h^2} = -2 \]

      同理 \(L_h y_j^2 = -2\),因此:

      \[L_h E_{ij} = \frac{Kh}{4} \cdot (-L_h(x_i^2 + y_j^2)) = \frac{Kh}{4} \cdot 4 = Kh \]

    • 对于非正则内点,虽然采用不等距差分,但差分算子作用于 \(x^2 + y^2\) 的结果仍然为 \(-4\),因此仍有 \(L_h E_{ij} = Kh\)

综上,比较函数满足:

\[\begin{cases} L_h E_{ij} = Kh, & (x_i,y_j) \in G_h \\ E_{ij} \geq 0, & (x_i,y_j) \in \Gamma_h \end{cases} \tag{2.98} \]

三、误差估计与收敛性结论

应用比较定理(定理2.3)于误差方程(2.97)和比较函数方程(2.98):

  • 在内点:\(|L_h e_{ij}| = |R_{ij}| \leq Kh = L_h E_{ij}\)
  • 在边界:\(|e_{ij}| = 0 \leq E_{ij}\)

因此有:

\[|e_{ij}| \leq E_{ij} = \frac{Kh}{4} \left( R^2 - x_i^2 - y_j^2 \right) \]

取最大值,得到误差的最大模估计:

\[\max_{\bar{G}_h} |e_{ij}| \leq \frac{K R^2 h}{4} \tag{2.99} \]

四、结论解读

  1. 收敛性:当 \(h \to 0\) 时,\(\max |e_{ij}| \to 0\),即五点差分格式的差分解一致收敛到精确解;
  2. 敛速估计:该估计给出的收敛阶为 \(O(h)\),这是一个保守估计。实际上,当解 \(u \in C^4(\bar{G})\) 时,非正则内点的数量为 \(O(h^{-1})\) 量级,其局部 \(O(h)\) 的误差对整体误差的贡献为 \(O(h^2)\),因此实际整体收敛阶为 \(O(h^2)\)
  3. 极值定理的作用:整个证明过程完全依赖于极值定理和比较定理,这再次验证了极值定理是椭圆型差分方程理论的核心工具。

2.5.4 习题 完整解答

习题1:一维差分方程的极值定理证明

题目:设 \(\bar{I}_h = \{x_i: i=0,1,\dots,N\}\),网格函数 \(y_i\) 满足差分算子:

\[l y_i = -(a_i y_{i-1} - b_i y_i + c_i y_{i+1}) + q_i y_i, \quad i=1,2,\dots,N-1 \]

其中 \(a_i,b_i,c_i,q_i\) 非负,且 \(a_i + c_i \leq b_i\)。证明当 \(l y_i \leq 0\)\(l y_i \geq 0\))时,\(y_i\) 不能在内点取正的极大值(负的极小值),除非 \(y_i\) 等于常数。

证明(反证法):
只证明 \(l y_i \leq 0\) 时不能取正极大值的情况,负极小值的证明类似。

  1. 反设:假设 \(y_i \not\equiv\) 常数,且在内点 \(i_0\) 处取正的极大值 \(M > 0\)
  2. 连通性应用:由于网格是连通的,必然存在相邻点(例如 \(i_0-1\))使得 \(y_{i_0-1} < M\)
  3. 差分算子展开

    \[\begin{aligned} l y_{i_0} &= -(a_{i_0} y_{i_0-1} - b_{i_0} y_{i_0} + c_{i_0} y_{i_0+1}) + q_{i_0} y_{i_0} \\ &= b_{i_0} y_{i_0} - a_{i_0} y_{i_0-1} - c_{i_0} y_{i_0+1} + q_{i_0} y_{i_0} \end{aligned} \]

  4. 不等式估计
    • 由于 \(y_{i_0} = M\) 是极大值,故 \(y_{i_0-1} \leq M\)\(y_{i_0+1} \leq M\)
    • 由假设 \(y_{i_0-1} < M\),且 \(a_{i_0} > 0\)
    • 因此:

      \[\begin{aligned} l y_{i_0} &> b_{i_0} M - a_{i_0} M - c_{i_0} M + q_{i_0} M \\ &= (b_{i_0} - a_{i_0} - c_{i_0} + q_{i_0}) M \\ &\geq 0 \end{aligned} \]

  5. 矛盾:这与条件 \(l y_i \leq 0\) 矛盾。因此假设不成立,\(y_i\) 不能在内点取正的极大值。

证毕


习题2:一维差分方程的最大模估计

题目:在上题中,若设 \(d_i = b_i - a_i - c_i + q_i > 0\)\(i=1,2,\dots,N-1\)),则差分方程

\[\begin{cases} l y_i = \varphi_i, & i=1,2,\dots,N-1 \\ y_0 = y_N = 0 \end{cases} \]

的解满足:

\[\max_i |y_i| \leq \max_i \frac{|\varphi_i|}{d_i} \]

证明(比较定理法):
构造比较函数 \(U_i = \max_k \frac{|\varphi_k|}{d_k}\)(常数函数)。

  1. 验证比较定理条件

    • 在内点:

      \[l U_i = (b_i - a_i - c_i + q_i) U_i = d_i U_i \geq d_i \cdot \frac{|\varphi_i|}{d_i} = |\varphi_i| = |l y_i| \]

    • 在边界:

      \[|y_0| = |y_N| = 0 \leq U_i \]

  2. 应用比较定理
    由一维比较定理(习题1的推论),得 \(|y_i| \leq U_i\) 对所有 \(i\) 成立,即:

    \[\max_i |y_i| \leq \max_i \frac{|\varphi_i|}{d_i} \]

证毕


习题3:一维差分方程的收敛阶估计

题目:利用上题估计差分方程(2.28)(2.29)解的收敛阶(假定 \(r=0, q \geq q_0 > 0, h_i \equiv h\))。

解答

  1. 误差方程建立
    \(r=0\) 且网格均匀(\(h_i \equiv h\))时,差分方程(2.28)简化为:

    \[-\frac{1}{h^2} \left( p_{i+\frac{1}{2}} u_{i+1} - (p_{i+\frac{1}{2}} + p_{i-\frac{1}{2}}) u_i + p_{i-\frac{1}{2}} u_{i-1} \right) + q_i u_i = f_i \]

    与习题2的差分算子对比,得:

    \[a_i = \frac{p_{i-\frac{1}{2}}}{h^2}, \quad c_i = \frac{p_{i+\frac{1}{2}}}{h^2}, \quad b_i = \frac{p_{i+\frac{1}{2}} + p_{i-\frac{1}{2}}}{h^2}, \quad q_i = q_i \]

    因此:

    \[d_i = b_i - a_i - c_i + q_i = q_i \geq q_0 > 0 \]

    满足习题2的条件。

  2. 截断误差分析
    \(u \in C^4[a,b]\) 时,均匀网格下的截断误差为 \(O(h^2)\),即存在常数 \(C\) 使得:

    \[|R_i(u)| \leq C h^2 \]

  3. 误差估计
    误差 \(e_i = u(x_i) - u_i\) 满足方程:

    \[\begin{cases} l e_i = R_i(u), & i=1,2,\dots,N-1 \\ e_0 = e_N = 0 \end{cases} \]

    应用习题2的最大模估计:

    \[\max_i |e_i| \leq \max_i \frac{|R_i(u)|}{d_i} \leq \frac{C h^2}{q_0} = O(h^2) \]

结论:该差分格式的差分解按最大模二阶收敛到精确解。


核心知识点归纳总结表

内容 核心结论 关键依据 重要意义
五点格式敛速估计 \(\max |e_{ij}| \leq C h\)(保守估计),实际为 \(O(h^2)\) 比较定理 + 二次比较函数 证明了五点格式的收敛性
一维极值定理 解的极值只能在边界取得,除非是常数 反证法 + 系数非负性 一维差分方程理论的基础
一维最大模估计 \(\max |y_i| \leq \max |\varphi_i|/d_i\) 比较定理 + 常数比较函数 误差估计的通用工具
一维收敛阶 均匀网格下二阶收敛 截断误差 \(O(h^2)\) + 最大模估计 验证了一维差分格式的精度

posted on 2026-06-03 20:03  Indian_Mysore  阅读(22)  评论(0)    收藏  举报

导航