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

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

对于椭圆型方程

$\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](x_i)\\ &=-\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})$

这部分就是常见的中心差分形式，用于逼近二阶导数，体现了均匀网格下差分格式的简洁性和规律性，也反映出在特定条件下格式与经典差分形式的联系