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

一、示例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.