目录
一、研究目标
二、理论推导
三、算例实现
四、结论
一、研究目标
我们在前几节中介绍了Poisson方程的边值问题,接下来对椭圆型偏微分方程的混合边值问题进行探讨,研究对象为:
其中,为矩形区域,为上的连续函数,是上的单位法向量,从而表示方向导数,为的连续函数且非负。对于矩形区域而言,其边界上的法向量没有统一的表达式,需要对四条边界线段分别讨论。可见:
在左边界,边界条件为;
在右边界,边界条件为;
在下边界,边界条件为;
在上边界,边界条件为。
于是,公式(1)可以具体写为:
二、理论推导
首先进行矩形区域的等距剖分,得到各网格节点,且,;,。然后弱化原方程,使其仅在离散节点上成立,从而有
将上式中的一阶、二阶偏导数分别用关于一阶导数的中心差商和关于二阶导数的中心差商来近似,得
用数值解代替精确解并忽略高阶项,得到以下数值格式:
其中,。该格式的局部截断误差为,同时需要处理其中的下标越界问题。由于函数的连续性,我们认为公式(2)中的第1式对i=0也成立,即有
再成公式(2)的第2式中解出代入上式,得
同样,设公式(2)中的第1式分别对成立,再分别从公式(2)的第3、4、5式中解出代入前面刚得到的方程,就可以处理掉越界下标,得到以下格式:
至此我们共有(m+1)x(n+1)个待求量,而现有(m-1)x(n-1)个关于内节点的方程,2(n-1)个关于左、右边界上的节点(不含端点)的方程及2(m-1)个关于上、下边界上的节点(也不含端点)的方程,还需要补充
个方程,也就是关于矩形区域的4个顶点的方程。为此,设公式(2)中第1式对成立,即
然后再从公式(2)的第2式和第4式中分别解出和,代入公式(3)中,得到左下顶点处的方程:
上式可以简化为
同样,设公式(2)中第1式分别对、和成立,然后再从公式(2)中第3式和第4式解出分别代入刚才得到的3个方程,就得到的右下顶点、左上顶点、和右上顶点处的方程,这样,我们就有了完整的处理带导数边界条件的椭圆型方程的数值格式:
记,有
为计算分别,可将上述方程组写成矩阵形式。
首先,将公式(4)中第4、6、7式写成:
上面的式子可以简记为
其中,,为m+1阶单位矩阵,且C为公式(5)最左端的三对角矩阵,为公式(5)右端的向量,。接着,公式(4)中的第1,2,3式可以写成
,
该式可以简记为
其中,为公式(6)中的三对角矩阵,为公式(6)右端的向量。最后,公式(4)的第5、8、9式可以写成
该式可以简记为
其中,D为公式(7)中的三对角矩阵,为公式(7)右端的向量。于是,由公式(5)、(6)、(7)可得到公式(4)写成块三对角矩阵为
采用Gauss-Seidel迭代法求解公式(8)。
三、算例实现
用差分格式-公式(8)求解椭圆型方程混合边值问题:
已知精确解为。分别取步长为和,输出6个节点和处的数值解,并给出误差,要求在各节点处最大误差的迭代误差限为。
代码如下:
#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,maxerr;double *x,*y,**u,**v,**lambda,*d,kexi,eta,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=0.0;xb=2.0;ya=0.0;yb=1.0;m=64;n=32;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*gamma*v[0][n-1]+2*beta*u[1][n])/(alpha+(kexi+eta)*lambda[0][n]);for(i=1;i<m;i++)v[i][n]=(d[i]+2*gamma*v[i][n-1]+beta*(v[i-1][n]+u[i+1][n]))/(alpha+eta*lambda[i][n]);v[m][n]=(d[m]+2*gamma*v[m][n-1]+2*beta*v[m-1][n])/(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=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(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)
{return (pi*pi-1)*exp(x)*sin(pi*y);
}
double lambda_function(double x, double y)
{return x*x+y*y;
}
double phi1(double y)
{return (1-y*y)*sin(pi*y);
}
double phi2(double y)
{return (5.0+y*y)*exp(2.0)*sin(pi*y);
}
double psi1(double x)
{return pi*exp(x);
}
double psi2(double x)
{return -pi*exp(x);
}
double exact(double x, double y)
{return exp(x)*sin(pi*y);
}
当时,计算结果如下:
m=32,n=16.
k=4323
(0.50,0.25), y=1.152179, err=1.3643e-02.
(1.00,0.25), y=1.911016, err=1.1100e-02.
(1.50,0.25), y=3.162159, err=6.8738e-03.
(0.50,0.50), y=1.638607, err=1.0115e-02.
(1.00,0.50), y=2.711255, err=7.0265e-03.
(1.50,0.50), y=4.479936, err=1.7526e-03.
当时,计算结果如下:
m=64,n=32.
k=16007
(0.50,0.25), y=1.162412, err=3.4097e-03.
(1.00,0.25), y=1.919341, err=2.7745e-03.
(1.50,0.25), y=3.167313, err=1.7193e-03.
(0.50,0.50), y=1.646193, err=2.5286e-03.
(1.00,0.50), y=2.716524, err=1.7575e-03.
(1.50,0.50), y=4.481249, err=4.3972e-04.
四、结论
从计算结果可知,当步长减小为1/2时,误差减小为原来的1/4,可见该数值格式是二阶收敛的。