有限差分法-二维泊松方程及其Matlab程序实现

2.2 偏微分方程的差分解法

2.2.1 二维泊松方程

考虑区域 $\Omega$ 上的二维泊松问题:

$$\{\begin{array}{l}-\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)u=f(x,y),(x,y)\in \Omega \\ {u|}_{\partial \Omega }=\phi (x,y),\end{array}\text{\hspace{1em}(2-16)}$$
其中, $\partial \Omega$ 为区域 $\Omega$ 的边界, $f(x,y)$ 和 $\phi (x,y)$ 为已知函数。 $u{\mid }_{\partial \Omega }=\phi (x,y)$ 描述了函数 $u$ 在 边界上的取值, 即所谓的边界条件。为简单起见, 这里仅讨论 $\Omega$ 为矩形的情况, 即 $a<x<b$, $c<y<d$ 。首先, 用间隔 $h_1=(b-a)/N、h_2=(d-c)/M$ 分别将 $x$ 轴上的区间 $[a,b]、y$ 轴上的区间 $[c,d]$ 划分为 $N、M$ 等分, 得 $x_i=a+i{h}_1、0⩽i⩽N$ 和 $y_j=c+j{h}_2、0⩽j⩽M$ 。再根据 $x_i$ 和 $y_j$ 对矩形 区域进行网格剖分, 如下图所示。横线与坚线的交点就是网格结点, 称位于边界 $\partial \Omega$ 上的 结点为边界结点, 其余结点为内结点。

在内结点处考虑泊松问题:
$$-\left[\frac{{\partial }^{2}u\left(x_i,y_j\right)}{\partial x^2}+\frac{{\partial }^{2}u\left(x_i,y_j\right)}{\partial y^2}\right]=f\left(x_i,y_j\right),1⩽i⩽N-1,1⩽j⩽M-1\text{\hspace{1em}(2-17)}$$
由泰勒公式, 有:
$$\begin{array}{c}\frac{{\partial }^{2}u\left(x_i,y_j\right)}{\partial x^2}=\frac{1}{h_1^2}\left[u\left(x_{i-1},y_j\right)-2u\left(x_i,y_j\right)+u\left(x_{i+1},y_j\right)\right]-\frac{h_1^2}{12}\frac{{\partial }^{4}u\left({\xi }_{ij},y_j\right)}{\partial x^4},\\ x_{i-1}<{\xi }_{ij}<x_{i+1}\end{array}\text{\hspace{1em}(2-18)}$$
$$\begin{array}{c}\frac{{\partial }^{2}u\left(x_i,y_j\right)}{\partial y^2}=\frac{1}{h_2^2}\left[u\left(x_i,y_{j-1}\right)-2u\left(x_i,y_j\right)+u\left(x_i,y_{j+1}\right)\right]-\frac{h_2^2}{12}\frac{{\partial }^{4}u\left(x_i,{\eta }_{ij}\right)}{\partial y^4},\\ y_{j-1}<{\eta }_{ij}<y_{j+1}\end{array}\text{\hspace{1em}(2-19)}$$
将以上两式代入式 $(2-17)$, 忽略小量 $O\left(h_1{}^2+h_2{}^2\right)$, 并使用简写 $u_{ij}=u\left(x_i,y_j\right)$, 得到差分 格式 (2-20)。所谓差分格式, 就是用几个相邻数值点的差商来替代方程中的导数或偏导数 的近似算法。
$$\begin{array}{c}-\frac{1}{h_2^2}{u}_{i,j-1}-\frac{1}{h_1^2}{u}_{i-1,j}+2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right)u_{i,j}-\frac{1}{h_1^2}{u}_{i+1,j}-\frac{1}{h_2^2}{u}_{i,j+1}=f\left(x_i,y_j\right),\\ 1⩽i⩽N-1,1⩽j⩽M-1\end{array}\text{\hspace{1em}(2-20)}$$
先定义向量 ${\mathbf{u}}_j:$
$${\mathbf{u}}_j={\left(u_{1j},u_{2j},\dots ,u_{N-1,j}\right)}^{\mathrm{T}},0⩽j⩽M\text{\hspace{1em}(2-21)}$$
再将差分格式 (2-20) 写为矩阵形式:
$$\mathbf{D}{\mathbf{u}}_{j-1}+\mathbf{C}{\mathbf{u}}_j+\mathbf{D}{\mathbf{u}}_{j+1}={\mathbf{f}}_j,1⩽j⩽M-1\text{\hspace{1em}(2-22)}$$
其中, 矩阵 $\mathbf{C}、\mathbf{D}、$ 向量 ${\mathbf{f}}_j$ 的定义如下, 注意向量 ${\mathbf{f}}_j$ 的首尾元素已包含了 $x=a$ 和 $x=b$ 处的边界条件。

$$\mathbf{C}=\left(\begin{array}{ccccc}2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right)& -\frac{1}{h_1^2}& & & \\ -\frac{1}{h_1^2}& 2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right)& -\frac{1}{h_1^2}& & \\ & \ddots & \ddots & \ddots & \\ & & -\frac{1}{h_1^2}& 2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right)& -\frac{1}{h_1^2}\\ & & & -\frac{1}{h_1^2}& 2\left(\frac{1}{h_1^2}+\frac{1}{h_2^2}\right)\end{array}\right)\text{\hspace{1em}(2-23)}$$
$$\mathbf{D}=\left(\begin{array}{cccc}-\frac{1}{h_2^2}& & & \\ & -\frac{1}{h_2^2}& & \\ & & \ddots & \\ & & -\frac{1}{h_2^2}& \\ & & & -\frac{1}{h_2^2}\end{array}\right),{\mathbf{f}}_j=\left(\begin{array}{c}f\left(x_1,y_j\right)+\frac{1}{h_1^2}\phi \left(x_0,y_j\right)\\ f\left(x_2,y_j\right)\\ ⋮\\ f\left(x_{N-2},y_j\right)\\ f\left(x_{N-1},y_j\right)+\frac{1}{h_1^2}\phi \left(x_N,y_j\right)\end{array}\right)\text{\hspace{1em}(2-24)}$$
式 (2-22) 还可以进一步写成如下的矩阵形式, 注意等号右边向量的首尾元素加入了 $y=c$ 和 $y=d$ 处的边界条件。

$$\left(\begin{array}{ccccc}C& D& & & \\ D& C& D& & \\ & \ddots & \ddots & \ddots & \\ & & D& C& D\\ & & & D& C\end{array}\right)\left(\begin{array}{c}{\mathbf{u}}_1\\ {\mathbf{u}}_2\\ ⋮\\ {\mathbf{u}}_{M-2}\\ {\mathbf{u}}_{M-1}\end{array}\right)=\left(\begin{array}{c}{\mathbf{f}}_1-\mathbf{D}{\mathbf{u}}_0\\ {\mathbf{f}}_2\\ ⋮\\ {\mathbf{f}}_{M-2}\\ {\mathbf{f}}_{M-1}-\mathbf{D}{\mathbf{u}}_M\end{array}\right)\text{\hspace{1em}(2-25)}$$
因此, 矩形 $\Omega$ 上的二维泊松问题就转化为上面形如 $\mathbf{A}\mathbf{u}=\mathbf{f}$ 的矩阵问题。对于这类问题, 最直接的解法就是先求 $\mathbf{A}$ 的逆矩阵 ${\mathbf{A}}^{-1}$, 然后 $\mathbf{u}={\mathbf{A}}^{-1}\mathbf{f}$ 即可, 这在 Matlab 中可用左除 \ 实现。注意式 (2-22）只针对内结点, 所以 $\mathbf{C}、\mathbf{D}$ 均是 $N-1$ 阶方阵, $\mathbf{A}$ 是 $(N-1)(M-1)$ 阶方 阵。下面给出一个实际例子。
$$\{\begin{array}{lll}-\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)u=\left({\pi }^2-1\right){\mathrm{e}}^x\sin (\pi y),& 0<x<2,& 0<y<1\\ u(0,y)=\sin (\pi y),& u(2,y)={\mathrm{e}}^2\sin (\pi y),& 0⩽y⩽1\\ u(x,0)=0,& u(x,1)=0,& 0<x<2\end{array}\text{\hspace{1em}(2-26)}$$
根据式 (2-25) 求解这个泊松问题, 代码如下:

主程序代码如下

clear; close all;
%生成网格上的坐标
h=0.1; x=[0:h:2]'; y=[0:h:1]';
N=length(x)-1; M=length(y)-1;
[X,Y]=meshgrid(x,y);
%解析解
u_analytical=exp(X).*sin(pi*Y);
X=X(2:M,2:N); Y=Y(2:M,2:N);
%生成f(x,y)的矩阵
f=(pi^2-1)*exp(X).*sin(pi*Y);
f(:,1)=f(:,1)+sin(pi*y(2:M))/h^2;
f(:,end)=f(:,end)+exp(2)*sin(pi*y(2:M))/h^2;
%构造矩阵D、C、A
e=ones(N-1,1);
C=1/h^2*spdiags([-e 4*e -e],[-1 0 1],N-1,N-1);
D=-1/h^2*eye(N-1);
e=ones(M-1,1);
A=kron(eye(M-1),C)+kron(spdiags([e e],[-1 1],M-1,M-1),D);
%左除求解
f=f'; u=zeros(M+1,N+1);
u(2:M,2:N)=reshape(A\f(:),N-1,M-1)';
u(:,1)=sin(pi*y); u(:,end)=exp(2)*sin(pi*y);
%画图
figure(1), spy(A,5)
figure(2), subplot(2,2,1), mesh(x,y,u)
xlabel x, ylabel y, zlabel u
subplot(2,2,2), mesh(x,y,u-u_analytical)
axis([0 2 0 1 0 0.04])
xlabel x, ylabel y, zlabel Error

程序输出结果如图所示, 左图为 $h_1=h_2=0.1$ 时泊松问题 (2-26) 的数值解, 右图为数值解与解析解 $u={\mathrm{e}}^x\sin (\pi y)$ 之间的误差。当然, 进一步缩小 $h_1$ 和 $h_2$ 会以增加运算量为代价降低这个误差。

为了使读者产生直观的理解, 下图显示了 $h_1=h_2=0.2$ 时方阵 $\mathbf{A}$ 的非零元素分布情况。