偏微分方程算法之九点紧差分法

一、研究目标

二、理论推导

三、算例实现

四、结论

一、研究目标

我们已经在专栏中介绍了椭圆型偏微分方程的五点菱形差分格式，这里我们继续以该方法为背景，探讨如何提高五点法的精度，即从二阶精度提升到四阶精度。

研究目标现继续以矩形区域 $\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(1)\\ u(x,y)=\varphi(x,y),(x,y)\in\partial \Omega=\Gamma \end{matrix}\right.$

二、理论推导

提高差分格式精度的一种有效方法是紧差分，即引入中间函数，对二阶偏导数进行更精细的逼近。

与五点菱形差分法类似，首先进行矩形区域的剖分，然后将原方程弱化，使其仅在离散节点处成立。在利用差商代替微商之前，引入中间函数，令

$\frac{\partial ^{2}u}{\partial x^{2}}=v(x,y),\frac{\partial^{2}u}{\partial y^{2}}=w(x,y)$

原方程在离散节点处满足

$-(v+w)|_{(x_{i},y_{j})}=f(x_{i},y_{j})\; (2)$

又有，当 $max_{(x,y)\in\Omega}[|\frac{\partial^{6}u(x,y)}{\partial x^{6}}|,|\frac{\partial^{6}u(x,y)}{\partial y^{6}}|]$ 有界时，

$\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}}-\frac{\Delta x^{2}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{i},y_{j})}+C_{1}(\Delta x)^{4}$

从而有

$v(x_{i},y_{j})=\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{i},y_{j})}=\frac{\delta ^{2}_{x}u(x_{i},y_{j})}{\Delta x^{2}}-\frac{\Delta x^{2}}{12}\frac{\partial^{2}v}{\partial x^{2}}|_{(x_{i},y_{j})}+C_{1}(\Delta x)^{4}$

$=\frac{\delta^{2}_{x}u(x_{i},y_{j})}{\Delta x^{2}}-\frac{\Delta x^{2}}{12}(\frac{v(x_{i-1},y_{j})-2v(x_{i},y_{j})+v(x_{i+1},y_{j})}{\Delta x^{2}})+C_{2}(\Delta x)^{4}$

也就是

$v(x_{i},y_{j})=\frac{\delta^{2}_{x}u(x_{i},y_{j})}{\Delta x^{2}}-\frac{1}{12}(v(x_{i-1},y_{j})+10v(x_{i},y_{j})+v(x_{i+1},y_{j}))+C_{2}(\Delta x)^{4}\; (3)$

由二维抛物型偏微分方程紧交替方向隐格式方法介绍中算子公式

$\varepsilon ^{2}_{x}v(x_{i},y,t)=\frac{1}{12}(v(x_{i-1},y,t)+10v(x_{i},y,t)+v(x_{i+1},y,t))$

可知公式（3）为

$\epsilon^{2}_{x}v(x_{i},y_{j})=\frac{\delta^{2}_{x}u(x_{i},y_{j})}{\Delta x^{2}}+C_{2}(\Delta x)^{4}\; (4)$

同理可得

$\epsilon^{2}_{y}w(x_{i},y_{j})=\frac{\delta^{2}_{y}u(x_{i},y_{j})}{\Delta y^{2}}+C_{3}(\Delta y)^{4}\; (5)$

将算子 $\epsilon^{2}_{x}\epsilon^{2}_{y}$ 作用到公式（2）上，有

$-\epsilon^{2}_{x}\epsilon^{2}_{y}(v+w)|_{(x_{i},y_{j})}=\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{i},y_{j})\; (6)$

然后将公式（4）、（5）一起代入公式（6），可得

$-\frac{\epsilon^{2}_{y}\delta^{2}_{x}u(x_{i},y_{j})}{\Delta x^{2}}-\frac{\epsilon^{2}_{x}\delta^{2}_{y}u(x_{i},y_{j})}{\Delta y^{2}}=\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{i},y_{j})+C_{2}(\Delta x)^{4}+C_{3}(\Delta y)^{4}$

用数值解代替精确解并忽略高阶项，可得紧差分格式：

$\left\{\begin{matrix} -\frac{\epsilon^{2}_{y}\delta^{2}_{x}u_{i,j}}{\Delta x^{2}}-\frac{\epsilon^{2}_{x}\delta^{2}_{y}u_{i,j}}{\Delta y^{2}}=\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{i},y_{j}),1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\; (7)\\ 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.$

该格式的局部截断误差为 $O(\Delta x^{4}+\Delta y^{4})$ ，将公式（7）中第1式整理化简，利用记号

$\frac{1}{\Delta x^{2}}=\beta,\frac{1}{\Delta y^{2}}=\gamma$

可得

$(\beta+\gamma)u_{i-1,j-1}+(10\gamma-2\beta)u_{i,j-1}+(\beta+\gamma)u_{i+1,j-1}+(10\beta-2\gamma)u_{i-1,j}-20(\beta+\gamma)u_{i,j}+(10\beta-2\gamma)u_{i+1,j}+(\beta+\gamma )u_{i-1,j+1}+(10\gamma-2\beta)u_{i,j+1}+(\beta+\gamma)u_{i+1,j+1}=-12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{i},y_{j})$

该紧格式每一步计算都涉及9个点，所以公式（7）称为九点紧差分格式。

该格式无法写成线性方程组 $\mathbf{A}x=\mathbf{b}$ 的形式，只能写成

$\begin{pmatrix} \eta_{2} & \xi & & & \\ \xi & \eta_{2} & \xi & 0 & \\ & & \ddots & & \\ & 0 & \xi & \eta_{2} & \xi \\ & & & \xi & \eta_{2} \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} -20\xi & \eta_{1} & & 0 & \\ \eta_{1} & -20\xi & \eta_{1} & & \\ & \ddots & \ddots & \ddots & \\ & & \eta_{1} & -20\xi & \eta_{1} \\ & 0 & & \eta_{1} & -20\xi \end{pmatrix}$ $\begin{pmatrix} u_{1,j}\\ u_{2,j}\\ \vdots\\ u_{m-2,j}\\ u_{m-1,j} \end{pmatrix}+$ $\begin{pmatrix} \eta_{2} & \xi & & & \\ \xi & \eta_{2} & \xi & 0 & \\ & & \ddots & & \\ & 0 & \xi & \eta_{2} & \xi \\ & & & \xi & \eta_{2} \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} -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{1},y_{j})-\xi (u_{0,j-1}+u_{0,j+1})-\eta_{1}u_{0,j}\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{2},y_{j})\\ \vdots\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{m-2},y_{j})\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{m-1},y_{j})-\xi(u_{m,j-1}+u_{m,j+1})-\eta_{1}u_{m,j} \end{pmatrix}$

其中， $\xi=\beta+\gamma,\eta_{1}=10\beta-2\gamma,\eta_{2}=10\gamma-2\beta$

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

则原数值格式可以简写为

$C(\mathbf{u}_{j-1}+\mathbf{u}_{j+1})+D\mathbf{u}_{j}=\mathbf{f}_{j},j=1,2,\cdot\cdot\cdot,n-1$

其中，

$C=\begin{pmatrix} \eta_{2} & \xi & & & \\ \xi & \eta_{2} & \xi & 0 & \\ & & \ddots & & \\ & 0 & \xi & \eta_{2} & \xi \\ & & & \xi & \eta_{2} \end{pmatrix}$ ， $D=\begin{pmatrix} -20\xi & \eta_{1} & & 0 & \\ \eta_{1} & -20\xi & \eta_{1} & & \\ & \ddots & \ddots & \ddots & \\ & & \eta_{1} & -20\xi & \eta_{1} \\ & 0 & & \eta_{1} & -20\xi \end{pmatrix}$

且 $\mathbf{f}_{j}=\begin{pmatrix} -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{1},y_{j})-\xi (u_{0,j-1}+u_{0,j+1})-\eta_{1}u_{0,j}\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{2},y_{j})\\ \vdots\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{m-2},y_{j})\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{m-1},y_{j})-\xi(u_{m,j-1}+u_{m,j+1})-\eta_{1}u_{m,j} \end{pmatrix}$

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

$\begin{pmatrix} D & C & & & \\ C & D & C & & \\ & \ddots & \ddots & \ddots & \\ & & C & D & C\\ & & & C & D \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}-C\mathbf{u}_{0}\\ \mathbf{f}_{2}\\ \vdots\\ \mathbf{f}_{n-2}\\ \mathbf{f}_{n-1}-C\mathbf{u}_{n} \end{pmatrix}$

上述方程组的高斯-赛德尔迭代公式为：

$u^{(k+1)}_{i,j}=-\frac{1}{20\xi}[-12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{i},y_{j})-\xi(u^{(k+1)}_{i-1,j-1}+u^{(k)}_{i-1,j+1}+u^{(k+1)}_{i+1,j-1}+u^{(k)}_{i+1,j+1})-\eta_{1}(u^{(k+1)}_{i-1,j}+u^{(k)}_{i+1,j})-\eta_{2}(u^{(k+1)}_{i,j-1}+u^{(k)}_{i,j+1})]\; (8)$

其中， $1\leqslant i \leqslant m-1, 1\leqslant j \leqslant n-1$ ，加括号的上标k表示迭代次数。

三、算例实现

采用九点紧差分格式求解椭圆型方程边值问题：

$\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/16$ 和 $\Delta x=\Delta y=1/32$ ，输出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>
#define pi 3.14159265359int 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,**g,**temp,kexi,eta1,eta2;double leftboundary(double y);double rightboundary(double y);double bottomboundary(double x);double topboundary(double x);double f(double x, double y);double **Gij(double *x, double *y, int m, int n);double exact(double x, double y);xa=0.0;xb=2.0;ya=0.0;yb=1.0;m=64;n=32;printf("m=%d, n=%d.\n",m,m);dx=(xb-xa)/m;dy=(yb-ya)/n;beta=1.0/(dx*dx);gamma=1.0/(dy*dy);kexi=beta+gamma;eta1=10*beta-2*gamma;eta2=10*gamma-2*beta;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;}g=Gij(x,y,m,n);for(i=0;i<=m;i++){for(j=0;j<=n;j++)temp[i][j]=u[i][j];}k=0;do{maxerr=0.0;for(i=1;i<m;i++){for(j=1;j<n;j++){temp[i][j]=(g[i][j]-kexi*(u[i-1][j-1]+temp[i-1][j+1]+u[i+1][j-1]+temp[i+1][j+1])-eta1*(u[i-1][j]+temp[i+1][j])-eta2*(u[i][j-1]+temp[i][j+1]))/(-20*kexi);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]));}for(i=0;i<=m;i++){free(u[i]);free(temp[i]);}free(u);free(temp);free(x);free(y);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 bottomboundary(double x)
{return 0.0;
}
double topboundary(double x)
{return 0.0;
}
double exact(double x, double y)
{return exp(x)*sin(pi*y);
}
double f(double x, double y)
{return (pi*pi-1)*exp(x)*sin(pi*y);
}
double **Gij(double *x, double *y, int m, int n)
{int i,j;double temp1,temp2,temp3,**ans;ans=(double**)malloc(sizeof(double*)*(m+1));for(i=0;i<=m;i++)ans[i]=(double*)malloc(sizeof(double)*(n+1));for(i=1;i<m;i++){for(j=1;j<n;j++){temp1=f(x[i-1],y[j-1])+10*f(x[i],y[j-1])+f(x[i+1],y[j-1]);temp2=f(x[i-1],y[j])+10*f(x[i],y[j])+f(x[i+1],y[j]);temp3=f(x[i-1],y[j+1])+10*f(x[i],y[j+1])+f(x[i+1],y[j+1]);ans[i][j]=-(temp1+temp3+10*temp2)/12.0;}}return ans;
}

相同的问题，与五点菱形差分法计算结果进行了对比，结果如下：

当 $\Delta x=\Delta y=1/16$ 时，计算结果如下：

+++++++++++++九点紧差分格式++++++++++++++
m=32, n=16.
k=747.
(0.50,0.25), y=1.165829, err=6.7140e-06.
(1.00,0.25), y=1.922127, err=1.1487e-05.
(1.50,0.25), y=3.169047, err=1.4319e-05.
(0.50,0.50), y=1.648731, err=9.4951e-06.
(1.00,0.50), y=2.718298, err=1.6245e-05.
(1.50,0.50), y=4.481709, err=2.0250e-05.

+++++++++++++五点菱形差分格式++++++++++++
m=32,n=16.
k=887
(0.50,0.25), y=1.169343, err=3.5207e-03.
(1.00,0.25), y=1.928138, err=6.0228e-03.
(1.50,0.25), y=3.176531, err=7.4986e-03.
(0.50,0.50), y=1.653700, err=4.9790e-03.
(1.00,0.50), y=2.726799, err=8.5175e-03.
(1.50,0.50), y=4.492294, err=1.0605e-02.

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

+++++++++++++九点紧差分格式++++++++++++++
m=64, n=32.
k=2790.
(0.50,0.25), y=1.165822, err=4.1552e-07.
(1.00,0.25), y=1.922116, err=7.1227e-07.
(1.50,0.25), y=3.169034, err=8.9066e-07.
(0.50,0.50), y=1.648722, err=5.8773e-07.
(1.00,0.50), y=2.718283, err=1.0074e-06.
(1.50,0.50), y=4.481690, err=1.2597e-06.

+++++++++++++五点菱形差分格式++++++++++++
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.