偏微分方程算法之椭圆型方程差分格式编程示例

一、示例1-五点菱形格式

1.1 C++代码

1.2 计算结果

二、示例2-九点紧差分格式

2.1 C++代码

2.2 计算结果

三、示例3-二阶混合边值

3.1 C++代码

3.2 计算结果

本专栏对椭圆型偏微分方程的三种主要差分方法进行了介绍，并给出相应格式的理论推导过程。为加深对差分格式的理解，分别对三种方法进行C++编程示例。

一、示例1-五点菱形格式

$\left\{\begin{matrix} -(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})=\frac{4y^{2}-2x^{2}}{(x^{2}+2y^{2})^{2}},1<x<2,0<y<3,\\ u(1,y)=ln(1+2y^{2}),u(2,y)=ln(4+2y^{2}),0\leqslant y\leqslant 1,\\ u(x,0)=2lnx,u(x,3)=ln(18+x^{2}),1<x<2 \end{matrix}\right.$

已知精确解为 $u(x,y)=ln(x^{2}+2y^{2})$ 。分别取两种剖分数：m=20，n=30和m=40，n=60，输出10个节点 $(1.25,0.5i)$ 和 $(1.75,0.5i),i=1,2,3,4,5$ 处的数值解，并给出误差。要求在各个节点处最大误差的迭代误差限为 $0.5\times10^{-10}$ 。

1.1 C++代码

#include <cmath>
#include <stdlib.h>
#include <stdio.h>int main(int argc, char* argv[])
{int m,n,i,j,k,num;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 bottomboundary(double x);double topboundary(double x);double f(double x, double y);double exact(double x, double y);xa=1.0;xb=2.0;ya=0.0;yb=3.0;m=20;n=30;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.0*(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];}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=fabs(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=n/6;num=m/4;for(j=k;j<n;j=j+k){printf("(1.25,%.3f), y=%f, err=%.4e.\n",y[j],u[num][j],fabs(exact(x[num],y[j])-u[num][j]));}num=3*m/4;for(j=k;j<n;j=j+k){printf("(1.75,%.3f), y=%f, err=%.4e.\n",y[j],u[num][j],fabs(exact(x[num],y[j])-u[num][j]));}for(i=0;i<=m;i++){free(u[i]);free(temp[i]);}free(x);free(y);return 0;
}double leftboundary(double y)
{return log(1.0+2*y*y);
}
double rightboundary(double y)
{return log(4.0+2*y*y);
}
double bottomboundary(double x)
{return 2*log(x);
}
double topboundary(double x)
{return log(18.0+x*x);
}
double f(double x, double y)
{double temp1,temp2,z;temp1=x*x; temp2=y*y;z=temp1+2*temp2;return (4*temp2-2*temp1)/(z*z);
}
double exact(double x, double y)
{return log(x*x+2*y*y);
}

1.2 计算结果

当m=20，n=30时，计算结果为：

m=20,n=30.
k=959.
(1.25,0.500), y=0.724037, err=1.1827e-04.
(1.25,1.000), y=1.270654, err=1.9108e-04.
(1.25,1.500), y=1.802202, err=7.9937e-05.
(1.25,2.000), y=2.257872, err=2.3280e-05.
(1.25,2.500), y=2.643516, err=4.1352e-06.
(1.75,0.500), y=1.270488, err=2.5584e-05.
(1.75,1.000), y=1.621992, err=1.3181e-04.
(1.75,1.500), y=2.023279, err=7.7668e-05.
(1.75,2.000), y=2.403589, err=2.7872e-05.
(1.75,2.500), y=2.744871, err=6.9853e-06.

当m=40，n=60时，计算结果为：

m=40,n=60.
k=3582.
(1.25,0.500), y=0.723948, err=2.9304e-05.
(1.25,1.000), y=1.270510, err=4.7781e-05.
(1.25,1.500), y=1.802142, err=1.9972e-05.
(1.25,2.000), y=2.257855, err=5.8033e-06.
(1.25,2.500), y=2.643513, err=1.0237e-06.
(1.75,0.500), y=1.270469, err=6.1963e-06.
(1.75,1.000), y=1.621893, err=3.2942e-05.
(1.75,1.500), y=2.023221, err=1.9426e-05.
(1.75,2.000), y=2.403568, err=6.9568e-06.
(1.75,2.500), y=2.744866, err=1.7374e-06.

二、示例2-九点紧差分格式

$\left\{\begin{matrix} -(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})=\frac{4y^{2}-2x^{2}}{(x^{2}+2y^{2})^{2}},1<x<2,0<y<3,\\ u(1,y)=ln(1+2y^{2}),u(2,y)=ln(4+2y^{2}),0\leqslant y\leqslant 1,\\ u(x,0)=2lnx,u(x,3)=ln(18+x^{2}),1<x<2 \end{matrix}\right.$

2.1 C++代码

#include <cmath>
#include <stdlib.h>
#include <stdio.h>int main(int argc, char*argv[])
{int m,n,i,j,k,num;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=1.0;xb=2.0;ya=0.0;yb=3.0;m=20;n=30;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);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=fabs(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=n/6;num=m/4;for(j=k;j<n;j=j+k){printf("(1.25,%.3f), y=%f, err=%.4e.\n",y[j],u[num][j],fabs(exact(x[num],y[j])-u[num][j]));}num=3*m/4;for(j=k;j<n;j=j+k){printf("(1.75,%.3f), y=%f, err=%.4e.\n",y[j],u[num][j],fabs(exact(x[num],y[j])-u[num][j]));}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 log(1.0+2*y*y);
}
double rightboundary(double y)
{return log(4.0+2*y*y);
}
double bottomboundary(double x)
{return 2*log(x);
}
double topboundary(double x)
{return log(18+x*x);
}
double f(double x, double y)
{double temp1, temp2, z;temp1=x*x;temp2=y*y;z=temp1+2*temp2;return (4*temp2-2*temp1)/(z*z);
}
double exact(double x, double y)
{return log(x*x+2*y*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=0;i<=m;i++){for(j=0;j<=n;j++)ans[i][j]=0.0;}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;
}

2.2 计算结果

当m=20，n=30时，计算结果为：

m=20,n=30.
k=805.
(1.25,0.500), y=0.723921, err=2.5068e-06.
(1.25,1.000), y=1.270463, err=4.0234e-07.
(1.25,1.500), y=1.802122, err=8.8970e-08.
(1.25,2.000), y=2.257849, err=6.0205e-08.
(1.25,2.500), y=2.643512, err=2.1371e-08.
(1.75,0.500), y=1.270463, err=8.8774e-07.
(1.75,1.000), y=1.621861, err=5.0648e-07.
(1.75,1.500), y=2.023202, err=1.3736e-10.
(1.75,2.000), y=2.403561, err=4.9714e-08.
(1.75,2.500), y=2.744864, err=2.2523e-08.

当m=40，n=60时，计算结果为：

m=40,n=60.
k=3012.
(1.25,0.500), y=0.723919, err=1.5248e-07.
(1.25,1.000), y=1.270463, err=2.0549e-08.
(1.25,1.500), y=1.802122, err=1.0963e-08.
(1.25,2.000), y=2.257849, err=8.4329e-09.
(1.25,2.500), y=2.643512, err=4.0188e-09.
(1.75,0.500), y=1.270463, err=5.2372e-08.
(1.75,1.000), y=1.621860, err=2.7195e-08.
(1.75,1.500), y=2.023202, err=5.0463e-09.
(1.75,2.000), y=2.403561, err=7.4797e-09.
(1.75,2.500), y=2.744864, err=3.9218e-09.

三、示例3-二阶混合边值

$\left\{\begin{matrix} -(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})=\frac{4y^{2}-2x^{2}}{(x^{2}+2y^{2})^{2}},1<x<2,0<y<3,\\ (\frac{\partial u(x,y)}{\partial x}-u)|_{(1,y)}=\frac{2}{1+2y^{2}}-ln(1+2y^{2}),0\leqslant y\leqslant 3,\\ (\frac{\partial u(x,y)}{\partial x}+u)|_{(2,y)}=\frac{2}{2+y^{2}}+ln(4+2y^{2}),0\leqslant y\leqslant 3,\\ (\frac{\partial u(x,y)}{\partial y}-u)|_{(x,0)} = \-2lnx,1\leqslant x\leqslant 2,\\ (\frac{\partial u(x,y)}{\partial y})|_{(x,3)}=\frac{12}{18+x^{2}}+ln(18+x^{2}),1\leqslant x\leqslant 2 \end{matrix}\right.$

3.1 C++代码

#include <cmath>
#include <stdlib.h>
#include <stdio.h>int main(int argc, char*argv[])
{int m, n, i, j, k, num;double xa, xb, ya, yb, dx, dy, alpha, beta, gamma, maxerr;double *x, *y, **u, **v, **lambda, kexi, eta, *d, temp;double f(double x, double y);double lambda_function(double x, double y);double phi1(double y);double phi2(double y);double psi1(double x);double psi2(double x);double exact(double x, double y);xa=1.0;xb=2.0;ya=0.0;yb=3.0;m=20;n=30;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);kexi=2.0/dx;eta=2.0/dy;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));v=(double**)malloc(sizeof(double*)*(m+1));lambda=(double**)malloc(sizeof(double*)*(m+1));for(i=0;i<=m;i++){u[i]=(double*)malloc(sizeof(double)*(n+1));v[i]=(double*)malloc(sizeof(double)*(n+1));lambda[i]=(double*)malloc(sizeof(double)*(n+1));}for(i=0;i<=m;i++){for(j=0;j<=n;j++){u[i][j]=0.0;v[i][j]=0.0;lambda[i][j]=lambda_function(x[i], y[j]);}}d=(double*)malloc(sizeof(double)*(m+1));k=0;do{maxerr=0.0;for(i=0;i<=m;i++)d[i]=f(x[i],y[0])-eta*psi1(x[i]);d[0]=d[0]-kexi*phi1(y[0]);d[m]=d[m]+kexi*phi2(y[0]);v[0][0]=(d[0]+2*gamma*u[0][1]+2*beta*u[1][0])/(alpha+(kexi+eta)*lambda[0][0]);for(i=1;i<m;i++)v[i][0]=(d[i]+2*gamma*u[i][1]+beta*(v[i-1][0]+u[i+1][0]))/(alpha+eta*lambda[i][0]);v[m][0]=(d[m]+2*gamma*u[m][1]+2*beta*v[m-1][0])/(alpha+(kexi+eta)*lambda[m][0]);for(j=1;j<n;j++){for(i=0;i<=m;i++)d[i]=f(x[i],y[j]);d[0]=d[0]-kexi*phi1(y[j]);d[m]=d[m]+kexi*phi2(y[j]);v[0][j]=(d[0]+gamma*(u[0][j+1]+v[0][j-1])+2*beta*u[1][j])/(alpha+kexi*lambda[0][j]);for(i=1;i<m;i++)v[i][j]=(d[i]+gamma*(v[i][j-1]+u[i][j+1])+beta*(v[i-1][j]+u[i+1][j]))/alpha;v[m][j]=(d[m]+gamma*(v[m][j-1]+u[m][j+1])+2*beta*v[m-1][j])/(alpha+kexi*lambda[m][j]);}for(i=0;i<=m;i++)d[i]=f(x[i],y[n])+eta*psi2(x[i]);d[0]=d[0]-kexi*phi1(y[n]);d[m]=d[m]+kexi*phi2(y[n]);v[0][n]=(d[0]+2*beta*u[1][n]+2*gamma*v[0][n-1])/(alpha+(kexi+eta)*lambda[0][n]);for(i=1;i<m;i++)v[i][n]=(d[i]+beta*(v[i-1][n]+u[i+1][n])+2*gamma*v[i][n-1])/(alpha+eta*lambda[i][n]);v[m][n]=(d[m]+2*beta*v[m-1][n]+2*gamma*v[m][n-1])/(alpha+(kexi+eta)*lambda[m][n]);for(i=0;i<=m;i++){for(j=0;j<=n;j++){temp=fabs(u[i][j]-v[i][j]);if(temp>maxerr)maxerr=temp;u[i][j]=v[i][j];}}k=k+1;}while((maxerr>0.5*1e-10)&&(k<=1e+8));printf("k=%d\n", k);k=n/6;num=m/4;for(j=k;j<n;j=j+k){printf("(1.25,%.3f), y=%f, err=%.4e.\n",y[j],u[num][j],fabs(exact(x[num],y[j])-u[num][j]));}num=3*m/4;for(j=k;j<n;j=j+k){printf("(1.75,%.3f), y=%f, err=%.4e.\n",y[j],u[num][j],fabs(exact(x[num],y[j])-u[num][j]));}for(i=0;i<=m;i++){free(u[i]);free(v[i]);free(lambda[i]);}free(u);free(v);free(lambda);free(x);free(y);free(d);return 0;
}double f(double x, double y)
{double temp1, temp2, z;temp1=x*x;temp2=y*y;z=temp1+2*temp2;return (4*temp2-2*temp1)/(z*z);
}
double lambda_function(double x, double y)
{return 1.0;
}
double phi1(double y)
{double z;z=1.0+2*y*y;return 2.0/z-log(z);
}
double phi2(double y)
{double z;z=2+y*y;return 2.0/z+log(2*z);
}
double psi1(double x)
{return -2*log(x);
}
double psi2(double x)
{double z;z=x*x+18.0;return 12.0/z+log(z);
}
double exact(double x, double y)
{return log(x*x+2*y*y);
}

3.2 计算结果

当m=20，n=30时，计算结果为：

m=20, n=30
k=4470
(1.25,0.500), y=0.723996, err=7.7043e-05.
(1.25,1.000), y=1.270860, err=3.9760e-04.
(1.25,1.500), y=1.802391, err=2.6918e-04.
(1.25,2.000), y=2.257989, err=1.3972e-04.
(1.25,2.500), y=2.643565, err=5.3582e-05.
(1.75,0.500), y=1.270387, err=7.5935e-05.
(1.75,1.000), y=1.622151, err=2.9080e-04.
(1.75,1.500), y=2.023479, err=2.7756e-04.
(1.75,2.000), y=2.403726, err=1.6475e-04.
(1.75,2.500), y=2.744937, err=7.3239e-05.

当m=40，n=60时，计算结果为：

m=40, n=60
k=16565
(1.25,0.500), y=0.723937, err=1.8621e-05.
(1.25,1.000), y=1.270562, err=9.9132e-05.
(1.25,1.500), y=1.802189, err=6.7202e-05.
(1.25,2.000), y=2.257884, err=3.4879e-05.
(1.25,2.500), y=2.643525, err=1.3353e-05.
(1.75,0.500), y=1.270443, err=1.9346e-05.
(1.75,1.000), y=1.621933, err=7.2431e-05.
(1.75,1.500), y=2.023271, err=6.9315e-05.
(1.75,2.000), y=2.403602, err=4.1144e-05.
(1.75,2.500), y=2.744882, err=1.8266e-05.