椭圆型方程有限差分法求解

椭圆型方程及有限差分法基本步骤

对于椭圆型方程
$$
\begin{cases}
-\Delta u = f & \Omega\subset\mathbb{R}^d\u = g & \partial\Omega
\end{cases}
$$
其中
$$
\Delta u=\sum\frac{\partial^2u}{\partial x_i^2}
$$
有限差分法的基本步骤：

1、区域离散化：划分网格，将求解区域离散成节点

2、微分算子离散化：用差分近似导数

3、离散方程组求解：求解得到的离散方程组

一维情形直接离散化方法

区域离散化

对于方程
$$
\begin{cases}
Lu=-u’’ + q(x)u = f(x),\quad a < x < b,q(x)\geq0\
u(a)=\alpha, u(b)=\beta \quad (Dirichlet 边值条件)

\end{cases}
$$

将 $[a,b]$ 分成 $N$ 等分，$\displaystyle h=\frac{b - a}{N}$ ，节点 $x_i = a+ih$ ，$i = 0,1,\cdots,N$

微分算子离散化

利用 Taylor 展开
$$
\begin{align}
u(x_{i + h}) &=u(x_i)+u’(x_i)h+\frac{u’‘(x_i)}{2!}h^2+\frac{u’‘’(x_i)}{3!}h^3+\frac{u^{(4)}(x_i)}{4!}h^4+O(h^5)\
u(x_{i - h})&=u(x_i)-u’(x_i)h+\frac{u’‘(x_i)}{2!}h^2-\frac{u’‘’(x_i)}{3!}h^3+\frac{u^{(4)}(x_i)}{4!}h^4+O(h^5)\

\end{align}
$$
两式相加并整理得
$$
\frac{u(x_{i + h})+u(x_{i - h})-2u(x_i)}{h^2}=u’‘(x_i)+\frac{u^{(4)}(x_i)}{12}h^2+O(h^3)
$$
记 $R_i(u)$ 为截断误差。原方程在 $x_i$ 处写成
$$
-\frac{u(x_{i - h})-2u(x_i)+u(x_{i + h})}{h^2}+q(x_i)u(x_i)=f(x_i)+R_i(u)
$$
舍去 $R_i(u)$ 得差分方程（中心差分格式）
$$
-\frac{u_{i - 1}-2u_i+u_{i + 1}}{h^2}+q_iu_i = f_i\
q_i = q(x_i)\
f_i = f(x_i)
$$
微分算子 $L$ 记为：
$$
Lu=-u’‘+q(x)u
$$
微分算子 $L_h$ 记为：
$$
L_h[u(x_i)]=-\frac{u(x_{i - h})-2u(x_i)+u(x_{i + h})}{h^2}+q(x_i)u(x_i)
$$
截断误差 $R_i(u)$ 记为：
$$
\begin{align}
R_i(u) &=L_h[u(x_i)]-Lu\
&=-\frac{u(x_{i - h})-2u(x_i)+u(x_{i + h})}{h^2}+u’’
\end{align}
$$

离散方程组求解

离散化方程组为
$$
\begin{cases}\displaystyle -\frac{u_{i - 1}-2u_i+u_{i + 1}}{h^2}+q_iu_i = f_i, & i = 1,\cdots,N - 1\u_0=\alpha, u_N=\beta\end{cases}
$$
写成矩阵向量形式
$$
AU = F\
U=(u_1,u_2,\cdots,u_{N - 1})^T
$$
其中
$$
A=\begin{bmatrix}\displaystyle
\frac{2}{h^2}+q_1 &\displaystyle -\frac{1}{h^2} & 0 & \cdots & 0 \
\displaystyle -\frac{1}{h^2} & \displaystyle \frac{2}{h^2}+q_2 &\displaystyle -\frac{1}{h^2} & \ddots & \vdots \
0 & \displaystyle -\frac{1}{h^2} & \ddots & \ddots & 0 \
\vdots& \ddots& \ddots &\displaystyle \frac{2}{h^2}+q_{N-2}& \displaystyle -\frac{1}{h^2} \
0 & \cdots& 0 & \displaystyle -\frac{1}{h^2} & \displaystyle \frac{2}{h^2}+q_{N-1}
\end{bmatrix}
$$
即
$$
-\frac{1}{h^2}\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}
\begin{bmatrix}
u_1 \
u_2 \
\vdots \
\vdots \
\vdots \
u_{N-1}
\end{bmatrix}+\begin{bmatrix}
q_1 & 0 & \cdots & 0 \
0 & q_2 & \ddots & \vdots \
\vdots & \ddots & \ddots & 0 \
0 & \cdots & 0 & q_{N-1}
\end{bmatrix}\begin{bmatrix}
u_1 \
u_2 \
\vdots \
\vdots \
u_{N-1}
\end{bmatrix}=\begin{bmatrix}
\displaystyle f_1+\frac{\alpha}{h^2} \
\displaystyle f_2 \
\vdots \
\vdots \
f_{N-2}\
\displaystyle f_{N-1}+\displaystyle \frac{\beta}{h^2}
\end{bmatrix}
$$

方法一、追赶法（Thomas 法）

对于三对角方程组
$$
-a_iU_{i - 1}+b_iU_i - c_iU_{i + 1}=d_i \quad \quad i = 1,\cdots,N - 1\
u_0 = 0 \quad u_N = 0
$$
且
$$
a_i>0 \quad b_i>0 \quad c_i>0 \quad b_i>a_i + c_i \quad (对角占优 )
$$
通过逐步消元得到
$$
U_i - e_iU_{i + 1}=f_i
$$
其中
$$
\begin{align}
e_i=\frac{c_i}{b_i - a_ie_{i - 1}}\
f_i=\frac{d_i + a_if_{i - 1}}{b_i - a_ie_{i - 1}}\
\end{align}
$$
$i = 1,\cdots,N - 1$ ，从 $i = 1$ 开始计算，令 $e_0 = 0$ ，$f_0 = 0$ ，最后回代求解 $U_i$ 。这种方法利用三对角矩阵的结构特点，高效求解离散方程组

方法二、迭代法

除追赶法外，还可使用迭代法（如 Jacobi、Gauss - Seidel 等）求解离散方程组

以方程
$$
\begin{cases}-u’’ + xu=(x - 1)e^x\u(0)=1, u(1)=e\end{cases}
$$

（真解 $u(x)=e^x$ ）为例

通过中心差分格式得到
$$
-\frac{u_{i - 1}-2u_i+u_{i + 1}}{h^2}+x_iu_i=(x_i - 1)e^{x_i}
$$
整理为
$$
-a_iu_{i - 1}+b_iu_i - c_iu_{i + 1}=d_i
$$
其中
$$
\begin{align}
a_i = \frac{1}{h^2} ,\quad b_i=\frac{2}{h^2}+x_i ,\quad c_i=\frac{1}{h^2} \
d_i=(x_i - 1)e^{x_i}+\frac{u_0}{h^2} \quad (i = 1 时)\
d_{N - 1}=(x_{N - 1} - 1)e^{x_{N - 1}}+\frac{u_N}{h^2}
\end{align}
$$
迭代法通过不断迭代更新未知量的值，逐步逼近方程组的解

理论问题

解的存在唯一性：分析差分方程是否存在唯一解，常通过研究离散方程组系数矩阵的性质（如对角占优、正定等）来判断。例如，若系数矩阵严格对角占优，可利用相关定理证明方程组存在唯一解

收敛性及收敛速度：研究当 $h \to 0$ 时，数值解 $u_i$ 是否收敛到精确解 $u(x_i)$ 以及收敛速度。涉及在不同范数意义下的收敛情况。通过分析截断误差、稳定性等因素，推导收敛阶数，比如中心差分格式在合适条件下对二阶导数的逼近具有 $O(h^2)$ 阶的收敛速度

记号与概念

网格点：$I_h$ 表示网格内点，$\overline{I}_h$ 表示网格内点 + 边界点

网格函数及其范数：定义网格函数 $u_h$

$C$ 范数
$$
\displaystyle |u_h|C=\max{1\leq i\leq N - 1}|u_i|
$$
用于衡量网格函数在离散点上的最大绝对值

$L^2$ 范数
$$
\displaystyle |u_h|0=(\sum{i = 1}^{N - 1}h u_i^2)^{\frac{1}{2}}
$$
从离散点加权平方和角度衡量函数大小

$H^1$ 范数
$$
\displaystyle |u(x)|_{H^1}=(\int_a^b u^2(x)dx+\int_a^b (u’(x))^2dx)^{\frac{1}{2}}
$$
用于衡量函数及其一阶导数在区间 $[a,b]$ 上的综合 “大小”

$H^1$ 半范数
$$
\displaystyle |u_h|1=(\sum{i = 1}^{N - 1}h\cdot(\frac{u_i - u_{i + 1}}{h})^2)^{\frac{1}{2}}
$$
反映网格函数离散导数的某种度量

对于网格函数 $u_h$ ，其 $H^1$ 范数
$$
\displaystyle |u_h|_1=(|u_h|_0^2+|u_h|_1^2)^{\frac{1}{2}}
$$
是对离散函数在离散意义下类似的度量

相容性、收敛性与稳定性定义

相容性：对于某一类光滑函数 $u \in \mathcal{M}$ ，若 $\displaystyle \lim_{h \to 0}|R_h(u)| = 0$ ，则称有限差分格式具有相容性。这里 $R_h(u)$ 是截断误差，反映了差分格式对原微分方程的逼近程度，当步长 $h$ 趋于 0 时，截断误差趋于 0 ，意味着差分格式在极限意义下趋近原方程

收敛性：当 $h$ 充分小后，若差分方程的解 $u_h$ 存在，且在某范数下 $\displaystyle \lim_{h \to 0}|u_h - u| = 0$ ，则称差分格式收敛。即随着步长减小，数值解在特定范数度量下趋近于精确解

其中真解满足的离散方程为
$$
L_h[u(x_i)]=f_i+R_i(u)
$$

数值解满足的离散方程为
$$
L_hu_i=f_i
$$

误差满足的离散方程为
$$
L_he_i=R_i(u)
$$
稳定性：对于差分方程 $$\displaystyle \begin{cases}L_hu_i = f_i\u_0 = u_N = 0\end{cases}$$ ，若存在与网格 $I_h$ 及 $f_h$ 无关的常数 $M$ 和 $h_0$ ，当 $0 < h < h_0$ 时，有 $|u_h| \leq M|f_h|$ ，则称差分方程关于右端稳定。稳定性保证了在一定条件下，右端项的扰动不会引起解的过大变化，且可推出差分方程解唯一，以及误差 $e_h$ 满足 $\displaystyle |e_h| \leq M|R_h(u)|$

重要定理

定理：若解充分光滑，相容性 + 稳定性 $\Rightarrow$ 收敛性（收敛阶与截断误差阶相同）。这是有限差分法的关键结论，表明在函数光滑性满足要求时，只要差分格式同时具备相容性和稳定性，就必然收敛，且收敛阶数与截断误差阶数一致

一维差分格式

直接差分法

一维方程及问题设定

考虑方程
$$
\begin{cases}
\displaystyle Lu = -\frac{d}{dx}(p(x)\frac{du}{dx})+r(x)\frac{du}{dx}+q(x)u = f(x),&a < x < b\u(a)=\alpha, u(b)=\beta
\end{cases}
$$
其中 $p(x)\geq p_{min}>0$

直接差分法步骤

区域离散化：可采用不均匀网格，将区间 $[a,b]$ 离散为 $a = x_0 < x_1 < \cdots < x_N = b$ ，记 $h_i = x_i - x_{i - 1}$ ，$h=\max h_i$ ，引入对偶剖分 $\displaystyle a = x_0 < x_{\frac{1}{2}} < x_1 < \cdots < x_{N - \frac{1}{2}} < x_N = b$ ，$x_{i - \frac{1}{2}}=\frac{1}{2}(x_{i - 1}+x_i)$

导数项离散化尝试：

尝试 1：令 $v = p(x)u’(x)$ ，用
$$
\displaystyle \frac{v(x_{i + 1})-v(x_{i - 1})}{h_{i + 1}+h_i}=\frac{p(x_{i + 1})u’(x_{i + 1})-p(x_{i - 1})u’(x_{i - 1})}{h_{i + 1}+h_i}
$$
近似 $v’(x_i)$ ，但 $u’(x_{i + 1})$ ，$u’(x_{i - 1})$ 还需进一步离散，会涉及 $x_{i + 2}$ ，$x_{i - 2}$ 等节点

尝试 2：引入中点 $x_{i\pm\frac{1}{2}}$ ，用
$$
\frac{v(x_{i + \frac{1}{2}})-v(x_{i - \frac{1}{2}})}{\frac{1}{2}h_{i + 1}+\frac{1}{2}h_i}=\frac{p(x_{i +\frac{1}{2} })u’(x_{i + \frac{1}{2}})-p(x_{i - \frac{1}{2}})u’(x_{i - \frac{1}{2}})}{\frac{1}{2}(h_{i + 1}+h_i)}
$$
近似 $v’(x_i)$ ，其中
$$
u’(x_{i + \frac{1}{2}})\approx\frac{u(x_{i + 1})-u(x_i)}{h_{i + 1}}
$$
$\displaystyle u’(x_{i - \frac{1}{2}})$ 同理。这种方式在对 $(p(x)u’)'(x_i)$ 进行差分离散时，主要用到 $x_{i - 1}$ ，$x_i$ ，$x_{i + 1}$ 三个节点，相比尝试 1 更具优势，后续可在此基础上进一步构建差分格式，将原微分方程转化为离散的差分方程，进而求解

差分格式构建

基于前面的离散化尝试，定义
$$
[u]i\approx\frac{u{i + 1}-u_{i - 1}}{h_{i + 1}+h_i}\
[(p(x)u’)‘]i\approx\frac{p{i+\frac{1}{2}}[u]{i+\frac{1}{2}}-p{i - \frac{1}{2}}[u]{i - \frac{1}{2}}}{\frac{1}{2}(h_i+h{i + 1})}
$$
进一步得到原微分方程的差分格式：
$$
\begin{align}
L_hu_i = -\frac{2}{h_i+h_{i + 1}}[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{r_i}{h_i+h_{i + 1}}(u_{i + 1}-u_{i - 1})+q_iu_i = f_i\
i &= 1,\cdots,N - 1
\end{align}
$$
边界条件为 $u_0=\alpha$ ，$u_N=\beta$ 。该差分格式可写成矩阵形式 $AU = F$ ，其中 $A = A_1+A_2+A_3$ 是三对角矩阵。一般情况下
$$
-(p(x)u’)‘+ru’+qu
$$
不是自伴算子，$A$ 为非对称矩阵；若 $r(x)\equiv0$ ，则是自伴算子，此时若网格均匀，$A$ 对称，若网格不均匀，$A$ 不对称但可对称化

局部截断误差分析

对 $u’(x_i)$ 进行 Taylor 展开分析：

由
$$
\begin{align}
u(x_{i + h_{i + 1}})&=u(x_i)+u’(x_i)h_{i + 1}+\frac{u’‘(x_i)}{2!}h_{i + 1}^2+\frac{u’‘’(x_i)}{3!}h_{i + 1}^3+\cdots\
u(x_{i - h_i})&=u(x_i)-u’(x_i)h_i+\frac{u’‘(x_i)}{2!}h_i^2-\frac{u’‘’(x_i)}{3!}h_i^3+\cdots
\end{align}
$$
可得
$$
u’(x_i)\approx\frac{u(x_{i + h_{i + 1}})-u(x_{i - h_i})}{h_{i + 1}+h_i}=u’(x_i)+\frac{u’‘(x_i)}{2}(h_{i + 1}-h_i)+O(h^2)
$$
对 $(p(x)u’(x))'$ 进行 Taylor 展开分析：
$$
\begin{align}
(p(x)u’(x))‘(x_i)&\approx\frac{2}{h_i+h_{i + 1}}[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}]\
&=(pu’)‘(x_i)+\frac{h_{i + 1}-h_i}{4}(pu’)‘’(x_i)+\frac{h_{i + 1}-h_i}{12}p(x_i)u’‘’(x_i)+O(h^2)
\end{align}
$$
综合截断误差：
$$
R_i(u)=(h_{i + 1}-h_i)[\frac{1}{4}(pu’)‘’(x_i)+\frac{1}{12}p(x)u’‘’(x_i)-\frac{1}{2}r(x_i)u’'(x_i)]+O(h^2)
$$
网格不均匀时：截断误差 $R_i(u)=O(h)$ 。这是因为在不均匀网格下，节点间距不一致，在 Taylor 展开分析导数近似时，由于步长差异，高阶项不能完全抵消，导致整体误差阶数为 $O(h)$ ，精度相对较低

网格均匀时：截断误差 $R_i(u)=O(h^2)$ 。当网格均匀，即 $h_i = h$ （$i = 1,\cdots,N - 1$ ）时，在对导数进行差分离散并通过 Taylor 展开分析误差过程中，一些含 $h$ 的高阶项相互抵消，使得误差阶数提高到 $O(h^2)$ ，差分格式具有更高精度

均匀网格下的差分格式形式

此时差分格式为
$$
L_hu_i = -\frac{1}{h^2}(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})+\frac{r_i}{2h}(u_{i + 1}-u_{i - 1})+q_iu_i = f_i
$$
当 $p(x)\equiv1$ 时，格式中
$$
-\frac{1}{h^2}(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})
$$
这部分就是常见的中心差分形式，用于逼近二阶导数，体现了均匀网格下差分格式的简洁性和规律性，也反映出在特定条件下格式与经典差分形式的联系