偏微分方程算法之五点菱形差分法

一、研究目标

上个专栏我们介绍了双曲型偏微分方程的主要算法及实现。从今天开始，我们在新的专栏介绍另一种形式偏微分方程-椭圆型的解法。

研究目标选取经典的二维椭圆型方程（也称泊松Poisson方程）：

$-\Delta u=-(\frac{\partial^{2}u(x,y)}{\partial x^{2}}+\frac{\partial^{2}u(x,y)}{\partial y^{2}})=f(x,y)$

当f=0时，就是著名的拉普拉斯（Laplace）方程。椭圆型方程在流体力学、弹性力学、电磁学、几何学和变分法中都广泛应用。现假设所要讨论的为矩形区域 $\Omega=[(x,y)|a\leqslant x\leqslant b,c\leqslant y\leqslant d]$ ，考虑以下Poisson方程的边值问题：

$\left\{\begin{matrix} -(\frac{\partial^{2}u(x,y)}{\partial x^{2}}+\frac{\partial^{2}u(x,y)}{\partial y^{2}})=f(x,y),(x,y)\in\Omega,\space\space(2)\\ u(x,y)=\varphi(x,y),(x,y)\in\partial \Omega=\Gamma \end{matrix}\right.$

固定边界的无厚薄膜，受外力作用后达到平衡状态时的位移函数u满足上述方程。一般情况下，公式（2）是很难直接用解析的方式计算精确解的，所以需要利用数值方法求解。

二、理论推导

首先介绍五点菱形差分格式，推导过程如下：

第一步：网格剖分。对矩形区域进行剖分，即在x方向对[a,b]进行步长为 $\Delta x$ 的等距剖分，分成m份，得到m+1个节点 $x_{i}=a+i\cdot\Delta x,i=0,1,\cdot\cdot\cdot,m$ ，其中 $\Delta x=(b-a)/m$ 。在y方向对[c,d]进行步长为 $\Delta y$ 的等距剖分，分成n份，得到n+1个节点 $y_{j}=c+j\cdot\Delta y,j=0,1,\cdot\cdot\cdot,n$ ，其中 $\Delta y=(c-d)/n$ 。然后用两族平行线 $x=x_{i},y=y_{j}$ 将区域 $\Omega$ 分成mn个小矩形，得到节点 $(x_{i},y_{j})$ ，如图所示。X表示边界节点，其余为内部节点。

第二步：弱化原方程，使得在离散点处成立，即

$\left\{\begin{matrix} -(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})|_{(x_{i},y_{j})}=f(x_{i},y_{j}),(x_{i},y_{j})\in\Omega,\\ u(x_{s},y_{t})=\varphi(x_{s},y_{t}),(x_{s},y_{t})\in\Gamma \end{matrix}\right.$

其中， $1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,s=0$ 或 $m$ 且 $0\leqslant t\leqslant n,t=0$ 或 $n$ 且 $0\leqslant s\leqslant m$ 。也就是 $(x_{i},y_{j})$ 为内节点， $(x_{s},y_{t})$ 为边界节点。

第三步：用差商近似微商，建立数值格式，即

$\frac{\partial^{2}u}{\partial x^{2}}(x_{i},y_{j})=\frac{u(x_{i-1},y_{j})-2u(x_{i},y_{j})+u(x_{i+1},y_{j})}{\Delta x^{2}}+O(\Delta x^{2})$

$\frac{\partial^{2}u}{\partial y^{2}}(x_{i},y_{j})=\frac{u(x_{i},y_{j-1})-2u(x_{i},y_{j})+u(x_{i},y_{j+1})}{\Delta y^{2}}+O(\Delta y^{2})$

将上面两式代入离散节点处的方程，可得

$\left\{\begin{matrix} -(\frac{u(x_{i-1},y_{j})-2u(x_{i},y_{j})+u(x_{i+1},y_{j})}{\Delta x^{2}}+\frac{u(x_{i},y_{j-1})-2u(x_{i},y_{j})+u(x_{i},y_{j+1})}{\Delta y^{2}})=\\ f(x_{i},y_{j})+C_{1}(\Delta x)^{2}+C_{2}(\Delta y)^{2},1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ u_{s,t}=\varphi(x_s),y_{t},s=0,m,0\leqslant t\leqslant n;t=0,n,0\leqslant s\leqslant m \end{matrix}\right.$

用数值解代替精确解并忽略高阶项，可得数值格式为

$\left\{\begin{matrix} -(\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{(\Delta x)^{2}}+\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{(\Delta y)^{2}})=f(x_{i},y_{j}),1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ u_{s,t}=\varphi(x_{s},y_{t}), s=0,m,0\leqslant t\leqslant n;t=0,n,0\leqslant s\leqslant m \end{matrix}\right.$

整理上式可得

$\left\{\begin{matrix} -\frac{u_{i-1,j}}{\Delta x^{2}}-\frac{u_{i+1,j}}{\Delta x^{2}}+2(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}})u_{i,j}-\frac{u_{i,j-1}}{\Delta y^{2}}-\frac{u_{i,j+1}}{\Delta y^{2}}=f(x_{i},y_{j}),1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ u_{s,t}=\varphi(x_{s},y_{t}), s=0,m,0\leqslant t\leqslant n;t=0,n,0\leqslant s\leqslant m \space\space(3) \end{matrix}\right.$

公式（3）每一步计算要涉及5个点，除中心点外其余4个点正好位于一个菱形的4个顶点，所以这个格式称为“五点菱形差分格式”，简称“五点格式”。

第四步：差分格式求解。公式（3）无法写成线性方程组 $Ax=b$ 的简单形式，只能写成

$-\frac{1}{\Delta y^{2}}\begin{pmatrix} 1 & & & & \\ & 1 & & 0 & \\ & & \ddots & & \\ & 0 & & 1 & \\ & & & & 1 \end{pmatrix}$ $\begin{pmatrix} u_{1,j-1}\\ u_{2,j-1}\\ \vdots\\ u_{m-2,j-1}\\ u_{m-1,j-1} \end{pmatrix}$ $-\frac{1}{\Delta y^{2}}\begin{pmatrix} 1 & & & & \\ & 1 & & 0 & \\ & & \ddots & & \\ & 0 & & 1 & \\ & & & & 1 \end{pmatrix}$ $\begin{pmatrix} u_{1,j+1}\\ u_{2,j+1}\\ \vdots\\ u_{m-2,j+1}\\ u_{m-1,j+1} \end{pmatrix}$ $\begin{pmatrix} 2(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}}) & -\frac{1}{\Delta x^{2}} & & 0 & \\ -\frac{1}{\Delta x^{2}} & 2(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}}) & -\frac{1}{\Delta x^{2}}& & \\ & \ddots & \ddots& \ddots & \\ & & -\frac{1}{\Delta x^{2}} & 2(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}}) & -\frac{1}{\Delta x^{2}}\\ & 0 & & -\frac{1}{\Delta x^{2}} & 2(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}}) \end{pmatrix}$ $\begin{pmatrix} u_{1,j}\\ u_{2,j}\\ \vdots\\ u_{m-2,j}\\ u_{m-1,j} \end{pmatrix}$ $=\begin{pmatrix} f(x_{1},y_{j})+\frac{1}{\Delta x^{2}}u_{0,j}\\ f(x_{2},y_{j})\\ \vdots\\ f(x_{m-2},y_{j})\\ f(x_{m-1},y_{j})+\frac{1}{\Delta x^{2}}u_{m,j} \end{pmatrix}$

记 $\mathbf{u}_{j}=(u_{1,j},u_{2,j},\cdot\cdot\cdot,u_{m-1,j})^{T},0\leqslant j\leqslant n$

且设 $2(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}})=\alpha,\frac{1}{\Delta x^{2}}=\beta,\frac{1}{\Delta y^{2}}=\gamma$

则数值格式可写为

$\begin{pmatrix} -\gamma & & & \\ & -\gamma & & \\ & & \ddots & \\ & & & -\gamma \end{pmatrix}$ $\begin{pmatrix} u_{1,j-1}\\ u_{2,j-1}\\ \vdots\\ u_{m-1,j-1} \end{pmatrix}$ $+\begin{pmatrix} \alpha & -\beta & & \\ -\beta & \alpha & -\beta & \\ & & \ddots & \\ & & -\beta & \alpha \end{pmatrix}$ $\begin{pmatrix} u_{1,j}\\ u_{2,j}\\ \vdots\\ u_{m-1,j} \end{pmatrix}$ $+\begin{pmatrix} -\gamma & & & \\ & -\gamma & & \\ & & \ddots & \\ & & & -\gamma \end{pmatrix}$ $\begin{pmatrix} u_{1,j+1}\\ u_{2,j+1}\\ \vdots\\ u_{m-1,j+1} \end{pmatrix}$ $=\begin{pmatrix} f(x_{1},y_{j})+\beta u_{0,j}\\ f(x_{2},y_{j})\\ \vdots\\ f(x_{m-2},y_{j})\\ f(x_{m-1},y_{j})+\beta u_{m,j} \end{pmatrix},1\leqslant j\leqslant n-1$

上式可简写为 $A(\mathbf{u}_{i-1}+\mathbf{u}_{i+1})+B\mathbf{u}_{j}=\mathbf{f}_{j},j=1,2,\cdot\cdot\cdot,n-1$ 。其中， $A=-\gamma I$ ，且 $I$ 为m-1阶单位矩阵：

$B=\begin{pmatrix} \alpha & -\beta & & 0 & \\ -\beta & \alpha & -\beta & & \\ & \ddots & \ddots & \ddots & \\ & & -\beta & \alpha & -\beta \\ & 0 & & -\beta & \alpha \\ \end{pmatrix}$ 且 $\mathbf{f}_{j}=\begin{pmatrix} f(x_{1},y_{j})+\beta u_{0,j}\\ f(x_{2},y_{j})\\ \vdots\\ f(x_{m-2},y_{j})\\ f(x_{m-1},y_{j})+\beta u_{m,j} \end{pmatrix}$

为解出此方程组，将未知量 $\mathbf{u}_{j}$ 按下标拉长为一个列向量，并写成块矩阵形式，有

$\begin{pmatrix} B & A & & & \\ A & B & A & & \\ & \ddots & \ddots & \ddots & \\ & & A & B & A\\ & & & A &B \end{pmatrix}$ $\begin{pmatrix} \mathbf{u}_{1}\\ \mathbf{u}_{2}\\ \vdots\\ \mathbf{u}_{n-2}\\ \mathbf{u}_{n-1} \end{pmatrix}$ $=\begin{pmatrix} \mathbf{f}_{1}-A\mathbf{u}_{0}\\ \mathbf{f}_{2}\\ \vdots\\ \mathbf{f}_{n-2}\\ \mathbf{f}_{n-1}-A\mathbf{u}_{n} \end{pmatrix} \space\space(4)$

公式（4）的特点是：系数矩阵对称、正定，且绝大多数都是零元素，每一行最多只有5个非零元素，为稀疏矩阵。对于阶数不高的线性方程组的求解，直接法非常有效，而对于阶数高、系数矩阵稀疏的线性方程组，若采用直接法求解，就需要存储大量零元素。为减少运算律、节省内存，通常采用迭代法进行求解。在二维抛物型、双曲型方程的初边值问题中都曾遇到过这一类方程组，因为存在求解上的困难，后来就直接借助新的思路用交替方向隐式方法去处理数值逼近，从而避免了上述问题的求解。但事实上，公式（4）还是可以通过迭代法处理的，相比二维抛物型、双曲型方程初边值问题，由于不存在时间变量，处理起来会简单许多。具体的迭代法以及相应理论推导在下节中介绍（包括Jacobi迭代、Gauss-Seidel迭代、SOR迭代）。

三、算例实现

用五点菱形格式求解椭圆型方程边值问题：

$\left\{\begin{matrix} -(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})=(\pi^{2}-1)e^{x}sin(\pi y),0<x<2,0<y<1,\\ u(0,y)=sin(\pi y),u(2,y)=e^{2}sin(\pi y),0\leqslant y\leqslant 1,\\ u(x,0)=u(x,1)=0,0<x<2 \end{matrix}\right.$

已知该问题精确解为 $u(x,y)=e^{x}sin(\pi y)$ 。分别取步长 $\Delta x=\Delta y=1/32$ 和 $\Delta x=\Delta y=1/64$ ，输出6个节点 $(0.5i,0.25)$ 和 $(0.5i,0.5),i=1,2,3$ 处的数值解和误差。要求在各节点处最大误差的迭代误差限为 $0.5\times10^{-10}$ 。

代码如下：（采用Gauss-Seidel迭代）

#include <cmath>
#include <stdlib.h>
#include <stdio.h>
double pi=3.14159265359;

int main(int argc, char* argv[])
{
        int m, n, i, j, k;
        double xa, xb, ya, yb, dx, dy, alpha, beta, gamma, err, maxerr;
        double *x, *y, **u, **temp;
        double leftboundary(double y);
        double rightboundary(double y);
        double topboundary(double x);
        double bottomboundary(double x);
        double f(double x, double y);
        double exact(double x, double y);

        xa=0.0;
        xb=2.0;
        ya=0.0;
        yb=1.0;
        m=128;
        n=64;
        printf("m=%d,n=%d.\n",m,n);

        dx=(xb-xa)/m;
        dy=(yb-ya)/n;
        beta=1.0/(dx*dx);
        gamma=1.0/(dy*dy);
        alpha=2*(beta+gamma);

        x=(double*)malloc(sizeof(double)*(m+1));
        for(i=0;i<=m;i++)
                x[i]=xa+i*dx;

        y=(double*)malloc(sizeof(double)*(n+1));
        for(j=0;j<=n;j++)
                y[j]=ya+j*dy;

        u=(double**)malloc(sizeof(double *)*(m+1));
        temp=(double**)malloc(sizeof(double *)*(m+1));
        for(i=0;i<=m;i++)
        {
                u[i]=(double*)malloc(sizeof(double)*(n+1));
                temp[i]=(double*)malloc(sizeof(double)*(n+1));
        }

        for(j=0;j<=n;j++)
        {
                u[0][j]=leftboundary(y[j]);
                u[m][j]=rightboundary(y[j]);

        }

        for(i=1;i<m;i++)
        {
                u[i][0]=bottomboundary(x[i]);
                u[i][n]=topboundary(x[i]);
        }

        for(i=1;i<m;i++)
                for(j=1;j<n;j++)
                        u[i][j]=0.0;

        for(i=0;i<=m;i++)
                for(j=0;j<=n;j++)
                        temp[i][j]=u[i][j];

         //Gauss-Seidel迭代
        k=0;
        do
        {
                maxerr=0.0;
                for(i=1;i<m;i++)
                {
                        for(j=1;j<n;j++)
                        {
                                temp[i][j]=(f(x[i],y[j])+beta*(u[i-1][j]+temp[i+1][j])+gamma*(u[i][j-1]+temp[i][j+1]))/alpha;
                                err=temp[i][j]-u[i][j];
                                if(err>maxerr)
                                        maxerr=err;
                                u[i][j]=temp[i][j];
                        }
                }
                k=k+1;
        }while(maxerr>0.5*1e-10);
        printf("k=%d\n",k);

        k=m/4;
        for(i=k;i<m;i=i+k)
        {
                printf("(%.2f,0.25), y=%f, err=%.4e.\n",x[i],u[i][n/4],fabs(exact(x[i],y[n/4])-u[i][n/4]));
        }

        k=m/4;
        for(i=k;i<m;i=i+k)
        {
                printf("(%.2f,0.50), y=%f, err=%.4e.\n",x[i],u[i][n/2],fabs(exact(x[i],y[n/2])-u[i][n/2]));
        }

        return 0;
}


double leftboundary(double y)
{
        return sin(pi*y);
}
double rightboundary(double y)
{
        return exp(1.0)*exp(1.0)*sin(pi*y);
}
double topboundary(double x)
{
        return 0.0;
}
double bottomboundary(double x)
{
        return 0.0;
}
double f(double x, double y)
{
        return (pi*pi - 1)*exp(x)*sin(pi*y);
}
double exact(double x, double y)
{
        return exp(x)*sin(pi*y);
}

步长为 $\Delta x=\Delta y=1/32$ 时，计算结果如下：

m=64,n=32.
k=3315
(0.50,0.25), y=1.166702, err=8.7958e-04.
(1.00,0.25), y=1.923620, err=1.5048e-03.
(1.50,0.25), y=3.170908, err=1.8751e-03.
(0.50,0.50), y=1.649965, err=1.2439e-03.
(1.00,0.50), y=2.720410, err=2.1281e-03.
(1.50,0.50), y=4.484341, err=2.6518e-03.

步长为 $\Delta x=\Delta y=1/64$ 时，计算结果如下：

m=128,n=64.
k=12332
(0.50,0.25), y=1.166042, err=2.1984e-04.
(1.00,0.25), y=1.922492, err=3.7612e-04.
(1.50,0.25), y=3.169502, err=4.6879e-04.
(0.50,0.50), y=1.649032, err=3.1090e-04.
(1.00,0.50), y=2.718814, err=5.3191e-04.
(1.50,0.50), y=4.482352, err=6.6297e-04.