微分方程数值解 第二章（2）

有限体积法求解椭圆型微分方程

守恒型与非守恒型方程形式

对于椭圆型微分方程
$$
Lu = -(pu’)‘+qu = f
$$
守恒形式为 $(p(x)u’)‘$ ，对应散度型算子 。非守恒形式为
$$
p(x)u’’ + p’(x)u’
$$
当 $p(x)$ 有间断点时，若对非守恒形式采用中心差分格式，如
$$
\frac{p_i(u_{i - 1}-2u_i+u_{i + 1})}{h^2}+\frac{p_{i + 1}-p_{i - 1}}{2h}\cdot\frac{u_{i + 1}-u_{i - 1}}{2h}=0
$$
其差分解可能不收敛到真解

以方程
$$
\begin{align}
&\begin{cases}(p(x)u’)’ = 0, & 0 < x < 1\u(0)=1, u(1)=0\end{cases} \
\
&p(x)=\begin{cases}p_1, & 0\leq x\leq\xi\p_2, & \xi < x < 1\end{cases}
\end{align}
$$

为例，精确解为
$$
u(x)=\begin{cases}1 - \alpha x, & 0\leq x\leq\xi\\beta(1 - x), & \xi < x\leq1\end{cases}
$$
其中
$$
\begin{align}
\alpha &=\frac{p_2}{p_1+(p_2 - p_1)\xi}\
\beta &=\alpha\cdot\frac{p_1}{p_2}
\end{align}
$$

有限体积法原理（以恒定温度场为例 ）

在区间 $[a,b]$ 上任取子区间 $[x^{(1)},x^{(2)}]$ ，基于能量守恒有
$$
-\int_{x^{(1)}}^{x^{(2)}}(p(x)u’)‘dx+\int_{x^{(1)}}^{x^{(2)}}qudx=\int_{x^{(1)}}^{x^{(2)}}fdx
$$
即
$$
w(x^{(1)})-w(x^{(2)})+\int_{x^{(1)}}^{x^{(2)}}qudx=\int_{x^{(1)}}^{x^{(2)}}fdx
$$
这里 $[x^{(1)},x^{(2)}]$ 称为控制容积，$w(x)=p(x)u’$ 表示热流量，$p(x)$ 为热传导系数

有限体积法的优势在于，方程中二阶导数不直接出现，从而减弱了对 $p(x)$ 和 $u(x)$ 光滑性的要求，允许 $p(x)$ 存在间断，在处理实际物理问题（如热传导等 ）中更具适用性

有限体积法构建差分格式步骤解析

守恒性体现与网格剖分

守恒型微分方程描述物理量守恒规律，差分格式应反映此特性。对区间进行网格剖分并形成对偶剖分，取对偶单元 $\displaystyle [x_{i - \frac{1}{2}},x_{i + \frac{1}{2}}]$ ，根据能量守恒有
$$
w(x_{i - \frac{1}{2}})-w(x_{i + \frac{1}{2}})+\int_{\displaystyle x_{i - \frac{1}{2}}}^{\displaystyle x_{i + \frac{1}{2}}}qudx=\int_{\displaystyle x_{i - \frac{1}{2}}}^{\displaystyle x_{i + \frac{1}{2}}}fdx
$$

积分近似与导数离散

对 $\displaystyle \int_{\displaystyle x_{i - \frac{1}{2}}}^{\displaystyle x_{i + \frac{1}{2}}}qudx$ 近似为 $\displaystyle u_i\int_{\displaystyle x_{i - \frac{1}{2}}}^{\displaystyle x_{i + \frac{1}{2}}}qdx$ 。由 $w(x)=p(x)u’$ 且 $w(x)$ 连续，改写 $\displaystyle \frac{du}{dx}=\frac{w(x)}{p(x)}$ ，在 $[x_{i - 1},x_i]$ 上积分得
$$
u(x_i)-u(x_{i - 1})=\int_{\displaystyle x_{i - 1}}^{\displaystyle x_i}\frac{w(x)}{p(x)}dx\approx w(x_{i - \frac{1}{2}})\int_{\displaystyle x_{i - 1}}^{\displaystyle x_i}\frac{1}{p(x)}dx
$$
进而
$$
w(x_{i - \frac{1}{2}})\approx\frac{u(x_i)-u(x_{i - 1})}{\displaystyle \int_{x_{i - 1}}^{x_i}\frac{1}{p(x)}dx}=a_i\frac{u(x_i)-u(x_{i - 1})}{h_i}
$$
其中
$$
a_i=\frac{h_i}{\displaystyle \int_{x_{i - 1}}^{x_i}\frac{1}{p(x)}dx}
$$
同理可得 $w(x_{i + \frac{1}{2}})$

差分格式形成与积分计算

代入得到差分格式
$$
-[a_{i + 1}\frac{u_{i + 1}-u_i}{h_{i + 1}}-a_i\frac{u_i - u_{i - 1}}{h_i}]+u_i\int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}}qdx=\int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}}f(x)dx
$$
涉及 $p(x)$、$q(x)$、$f(x)$ 的积分，可用中点公式或梯形公式计算

中点公式下
$$
\begin{align}
a_i&\approx p_{i - \frac{1}{2}}\
\int_{\displaystyle x_{i - \frac{1}{2}}}^{\displaystyle x_{i + \frac{1}{2}}}qdx&\approx\frac{h_i+h_{i + 1}}{2}q_i\
\int_{\displaystyle x_{i - \frac{1}{2}}}^{\displaystyle x_{i + \frac{1}{2}}}f(x)dx&\approx\frac{h_i+h_{i + 1}}{2}f_i
\end{align}
$$
最终差分格式为
$$
-[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}]+\frac{h_i+h_{i + 1}}{2}q_iu_i=\frac{h_i+h_{i + 1}}{2}f_i
$$

注：若 $p$，$q$，$f$ 有第一类间断点，则在间断点的取值应是左右极限的算术平均

有限体积法相关特性及边界条件处理

有限体积法优点

降低光滑性要求：相比直接差分法，有限体积法对函数的光滑性要求更低，在处理 $p(x)$、$q(x)$ 等系数存在间断的情况时更具优势，能有效避免因系数间断导致的差分格式不收敛问题

保持物理量守恒：从物理守恒定律出发构建差分格式，能精确反映物理量在控制容积内的守恒特性，这在模拟实际物理过程（如热传导、流体流动等 ）中至关重要

方便处理网格及边界条件：可灵活适应任意网格形状和结构，在处理复杂区域和不同类型边界条件时更加便捷

边界条件处理

常见边界条件类型

Dirichlet 条件：$u(a)=\alpha$ ，直接给定边界点的函数值

Neumann 条件：$u’(a)=\alpha$ ，给定边界点的导数

Robin 条件：$u’(a)=d_0u(a)+d_1$ ，是函数值和导数的线性组合。守恒型方程常见边界条件为
$$
-p(a)u’(a)=d_0u(a)+d_1
$$

离散化及问题：

以 $\displaystyle u’(a)\approx\frac{u(x_1)-u(x_0)}{h}$ 离散导数，这种离散方式存在不足，一是截断误差比内点低，二是可能破坏差分方程的对称性

特殊控制容积处理：

取特殊控制容积 $[x^{(1)},x^{(2)}]=[x_0,x_{\frac{1}{2}}]$ ，根据守恒方程
$$
w(a)-w(x_{\frac{1}{2}})+\int_{x_0}^{x_{\frac{1}{2}}}qudx=\int_{x_0}^{x_{\frac{1}{2}}}fdx
$$
对 $-p(a)u’(a)=-w(a)$ 进行离散，得到
$$
-a_1\frac{u_1 - u_0}{h}+u_0\cdot\frac{h}{2}d_0-\frac{h}{2}d_0=d_0u_0 + d_1
$$
右边边界同理

当网格均匀且系数光滑时，边界截断误差可达到 $O(h^2)$ ，通过这种方式能更合理地处理边界条件，提高数值解的精度和稳定性

矩形网的差分格式

二维 Poisson 方程及边界条件

方程为
$$
-\Delta u = f(x,y) ，(x,y)\in G ，\Gamma=\partial G
$$
其中
$$
\Delta u = \sum\frac{\partial^2u}{\partial x_i^2}
$$
边界条件：

Dirichlet 条件：$\displaystyle u|_{\Gamma}=\alpha(x,y)$ ，直接给定边界上函数值

Neumann 条件：$\displaystyle \frac{\partial u}{\partial n}|_{\Gamma}=\beta(x,y)$ ，给定边界法向导数

Robin 条件：$\displaystyle (\frac{\partial u}{\partial n}+ku)|_{\Gamma}=\gamma(x,y)$ ，是法向导数与函数值的线性组合

规则区域 $G$ 适合有限差分法，不规则区域适合有限元或有限体积法。面临的困难有网格剖分及其数据结构构建，以及边界条件的妥善处理

**五点差分格式 **

区域离散化：设 $(x_i,y_j)$ 为正则内点，对区域进行离散

导数离散化：对于 $\displaystyle \Delta u = u_{xx}+u_{yy}$ ，记 $u_{ij}$ 为 $u(x_i,y_j)$ 的近似值，$f_{ij}=f(x_i,y_j)$

$u_{xx}$ 离散化为
$$
\displaystyle \frac{u_{i + 1,j}-2u_{ij}+u_{i - 1,j}}{h_1^2}
$$
这里 $h_1$ 是 $x$ 方向步长

$u_{yy}$ 离散化为
$$
\frac{u_{i,j - 1}-2u_{ij}+u_{i,j + 1}}{h_2^2}
$$
$h_2$ 是 $y$ 方向步长

差分方程构建

差分方程
$$
-\Delta_hu_{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}
$$
即五点差分格式，可简写为
$$
-\Delta_hu_h = f_h
$$
当 $h_1 = h_2 = h$ （正方形网格 ）时

格式为
$$
-\frac{1}{h^2}(u_{i - 1,j}+u_{i + 1,j}+u_{i,j - 1}+u_{i,j + 1}-4u_{ij})=f_{ij}
$$

局部截断误差分析

局部截断误差
$$
R_{ij}(u)=(\Delta u)(x_i,y_j)-\Delta_hu(x_i,y_j)
$$
利用真解的 Taylor 展开
$$
\begin{align}
\frac{u(x_{i - 1},y_j)-2u(x_i,y_j)+u(x_{i + 1},y_j)}{h_1^2}=u_{xx}(x_i,y_j)+\frac{u_{x}^{(4)}(x_i,y_j)}{12}h_1^2+O(h_1^4)\
\frac{u(x_i,y_{j - 1})-2u(x_i,y_j)+u(x_i,y_{j + 1})}{h_2^2}=u_{yy}(x_i,y_j)+\frac{u_{y}^{(4)}(x_i,y_j)}{12}h_2^2+O(h_2^4)
\end{align}
$$
综合可得
$$
R_{ij}(u)=-\frac{1}{12}[u_{x}^{(4)}h_1^2+u_{y}^{(4)}h_2^2]|_{(x_i,y_j)}+O(h^4)
$$
表明五点差分格式具有 $O(h^2)$ 阶精度（当 $h_1 = h_2 = h$ 时 ）

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

利用有限体积法推导五点差分格式，涉及矩形网格和对偶剖分，取控制容积 $G_{ij}=ABCD$ (顶点 $A:(x_{i - \frac{1}{2}},y_{i - \frac{1}{2}})$ ，$B:(x_{i + \frac{1}{2}},y_{i - \frac{1}{2}})$ ，$C:(x_{i + \frac{1}{2}},y_{i + \frac{1}{2}})$ ，$D:(x_{i - \frac{1}{2}},y_{i + \frac{1}{2}})$)

在控制容积 $G_{ij}$ 上，根据方程 $-\Delta u = f$ 进行积分：
$$
\int_{G_{ij}}-\Delta u dxdy=\int_{G_{ij}}f dxdy
$$
右侧
$$
\int_{G_{ij}}f dxdy\approx f_{ij}\cdot h_1\cdot h_2
$$
左侧利用 Green 公式
$$
\begin{align}
\int_{G_{ij}}-\Delta u dxdy &=-\int_{\partial G_{ij}}\frac{\partial u}{\partial n}dS
\&=-(\int_{AB}+\int_{BC}+\int_{CD}+\int_{DA}\frac{\partial u}{\partial n}dS)
\end{align}
$$
其中 $\displaystyle \frac{\partial u}{\partial n}$ 是 $u$ 沿 $\partial G_{ij}$ 的外法向导数

以 $AB$ 边为例
$$
\int_{AB}\frac{\partial u}{\partial n}dS\approx\frac{u(x_i,y_{j + 1})-u(x_i,y_j)}{h_2}\cdot h_1
$$
同理计算其他边，整理可得
$$
-\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}
$$
即五点差分格式

矩形区域离散化及数据结构

对于矩形区域 $G=(a,b)\times(c,d)$ ，边界条件 $\displaystyle u|_{\Gamma}=g(x,y)$ ，在 $[a,b]$ 上 $M$ 等分（步长 $h_1$ ），在 $[c,d]$ 上 $N$ 等分（步长 $h_2$ ）进行离散

数据结构：将 $u_{ij}$ 按行从下到上排成一个列向量 $\displaystyle u = [u_{11},u_{21},\cdots,u_{M - 1,1},u_{12},\cdots,u_{M - 1,2},\cdots,u_{M - 1,N}]^T$

按照五点差分格式，对区域内各节点 $(i,j)$ 展开方程。例如：

对于节点 $(1,1)$ ：
$$
-h_1^{-2}(-2u_{11}+u_{21})-h_2^{-2}(-2u_{11}+u_{12})=f_{11}+h_2^{-2}u_{10}+h_1^{-2}u_{01}
$$
这里考虑了边界点对内部点的影响（边界条件代入 ）

对其他节点也按类似方式展开，并也按 $a_{ij}$ 排列，形成一系列方程

矩阵结构分析

$\displaystyle \frac{\partial^2}{\partial x^2}$ 部分

以 $\displaystyle \frac{\partial^2}{\partial x^2}$ 对应的差分部分来看，其矩阵形式具有分块结构。每一块对应 $x$ 方向相邻节点的关系，是一个三对角矩阵形式。整体关于 $\displaystyle \frac{\partial^2}{\partial x^2}$ 部分的矩阵是 $(N - 1)$ 块的分块矩阵 。例如在均匀网格下，对于 $x$ 方向上相邻节点 $i - 1$，$i$，$i + 1$ ，差分格式中 $x$ 方向二阶导数项 $-h_1^{-2}(u_{i - 1,j}-2u_{ij}+u_{i + 1,j})$ 反映在矩阵中对应元素体现相邻节点耦合关系，这些元素按节点顺序排列构成三对角块结构

其结构大致如下
$$
-\frac{1}{h^2}\begin{bmatrix}
A & 0 & \cdots & 0 \
0 & A & \ddots & \vdots \
\vdots & \ddots & \ddots & 0 \
0 & \cdots & 0 & A \
\end{bmatrix}{(N - 1)\times(N - 1)}
$$
其中
$$
A=\begin{bmatrix}
-2 & 1 & 0 & \cdots & 0 \
1 & -2 & 1 & \ddots & \vdots \
0 & 1 & \ddots & \ddots & 0 \
\vdots & \ddots & \ddots & -2 & 1 \
0 & \cdots & 0 & 1 & -2 \
\end{bmatrix}{(M - 1)\times(M - 1)}
$$
$\displaystyle \frac{\partial^2}{\partial y^2}$ 部分

同理分析，其结构大致如下
$$
-\frac{1}{h^2}\begin{bmatrix}
-2I & I & 0 & \cdots & 0 \
I & -2I & I & \ddots & \vdots \
0 & I & \ddots & \ddots & 0 \
\vdots & \ddots & \ddots & -2I & I \
0 & \cdots & 0 & I & -2I \
\end{bmatrix}{(N - 1)\times(N - 1)}
$$
其中
$$
I=\begin{bmatrix}
1 & 0 & \cdots & 0 \
0 & 1 & \ddots & \vdots \
\vdots & \ddots & \ddots & 0 \
0 & \cdots & 0 & 1 \
\end{bmatrix}{(N - 1)\times(N - 1)}
$$

通过克罗内克积（$\otimes$ ）构建相关矩阵：

定义
$$
d_{xx}=\frac{1}{h_1^2}\begin{bmatrix}-2&1\1&-2&\ddots\&\ddots&\ddots&1\&&1&-2\end{bmatrix}{(M - 1)\times(M - 1)}
$$
$I{N - 1}$ 为 $(N - 1)$ 阶单位矩阵
$$
\begin{align}
D_{xx}&=I_{N - 1}\otimes d_{xx}\
&=\frac{1}{h_2^2}\begin{bmatrix}
d_{xx} & 0 & \cdots & 0 \
0 & d_{xx} & \ddots & \vdots \
\vdots & \ddots & \ddots & 0 \
0 & \cdots & 0 & d_{xx} \
\end{bmatrix}
\end{align}
$$
得到 $\displaystyle \frac{\partial^2}{\partial x^2}$ 差分部分的矩阵形式

类似地，对于 $\displaystyle \frac{\partial^2}{\partial y^2}$ ，定义
$$
d_{yy}=\frac{1}{h_2^2}\begin{bmatrix}-2&1\1&-2&\ddots\&\ddots&\ddots&1\&&1&-2\end{bmatrix}{(N - 1)\times(N - 1)}
$$
$I{M - 1}$ 为 $(M - 1)$ 阶单位矩阵
$$
\begin{align}
D_{yy}&=d_{yy}\otimes I_{M - 1}\
&=\frac{1}{h_2^2}\begin{bmatrix}-2I_{M - 1}&I_{M - 1}\I_{M - 1}&-2I_{M - 1}&\ddots\&\ddots&\ddots&I_{M - 1}\&&I_{M - 1}&-2I_{M - 1}\end{bmatrix}
\end{align}
$$

$\displaystyle -(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2})$ 的离散化矩阵

由前面构建的 $D_{xx}$ 和 $D_{yy}$ ，可得
$$
A = -(D_{xx}+D_{yy})
$$
从 $D_{xx}$ 和 $D_{yy}$ 的结构可知，$A$ 是五对角矩阵。这是因为 $D_{xx}$ 体现 $x$ 方向节点耦合，$D_{yy}$ 体现 $y$ 方向节点耦合，二者叠加后，矩阵 $A$ 中每个节点除了与上下左右相邻节点有耦合关系（对应五对角结构 ）外，还可能因边界条件和离散方式存在一些特殊关联

右端项向量 $\overline F$ 的处理

通过特定规则对原右端项向量 $F$ 进行修改得到 $\overline{F}$ ，如
$$
\begin{align}
\overline{F}(1,:)&=F(1,:)+h_1^{-2}[g_{0,1},\cdots,g_{0,N - 1}]\
\overline{F}(:,1)&=F(:,1)+h_2^{-2}[g_{1,0},\cdots,g_{M - 1,0}]^T\
\overline{F}(M-1,:)&=F(M-1,:)+h_1^{-2}[g_{0,1},\cdots,g_{0,N - 1}]\
\overline{F}(:,N-1)&=F(:,N-1)+h_2^{-2}[g_{1,N},\cdots,g_{M - 1,N}]^T
\end{align}
$$
这些操作是为了将边界条件 $g(x,y)$ 合理融入右端项，使得离散后的方程组 $Au = \overline{F}$ 能准确反映原二维 Poisson 方程及边界条件，从而保证数值解的准确性

矩阵形式方程求解

向量形式：将 $u_{ij}$ 按特定方式排列成向量 $u$ ，方程组为 $Au = \overline{F}$

矩阵形式：把 $u_{ij}$ 看作矩阵 $U$ ，以第一列为例，有
$$
-d_{xx}\cdot U(:,1)-U\cdot d_{yy}(:,1)=\overline{F}(:,1)
$$
一般形式为
$$
-d_{xx}\cdot U - U\cdot d_{yy}=\overline{F}
$$
$U$ 是 $(M - 1)\times(N - 1)$ 矩阵 。当边界条件 $g(x,y)\equiv0$ 时，$\overline{F}=F$ 。这种矩阵 - 向量形式的方程便于利用线性代数中的求解方法