目录
一、研究目标
二、理论推导
2.1 三层显格式建立
2.2 三层显格式改进
三、算例实现
3.1 一阶显格式
3.2 二阶显格式
一、研究目标
介绍完一阶双曲型偏微分方程的几种差分格式后,我们继续探讨二阶方程的差分格式。这里以非齐次二阶双曲型偏微分方程的初边值问题为研究对象(无法用解析解求解其精确解):
公式(1)中u表示一个与时间t和位置x有关的待求波函数,及方程右端项函数都是已知函数,是非零常数。公式(1)的第2式是初始条件,第3式是边界条件。若,则对应的初边值问题与一阶对流方程有类似属性,只不过公式(1)有两条特征线,精确解在点处的值的依赖区域为,从而也有相应的CFL条件。
二、理论推导
2.1 三层显格式建立
差分步骤与之前的抛物型方程差分中介绍的相同。
第一步:网格剖分。对矩形求解域进行等距剖分,即
空间步长为。
第二步:弱化原方程。将原来的连续方程离散到网格节点上成立,得到离散方程:
第三步:偏导数用差商近似。用二阶中心差商近似公式(2)中的二阶偏导数,即
初始条件中的一阶偏导数用向前差商近似,即
将上面三式代入公式(2)中,用数值解代替精确解并忽略高阶项,可得差分格式为
记,可以将公式(3)整理为公式(4)所示的三层显格式:
公式(4)的局部截断误差为。
2.2 三层显格式改进
事实上,公式(4)的精度比较低,原因是在对初始条件进行一阶偏导数近似时,采用的是精度较低的向前差商。为了提高精度,可以考虑利用中心差商对时间的一阶偏导数进行近似,即
得到公式(1)的数值格式为
这样处理引入虚拟越界点,如果认为公式(5)中第1式对k=0也成立,即
得出
于是可以在公式(5)第2式中消去,得到
从而得到改进的三层格式:
该格式的局部截断误差为。
三、算例实现
计算双曲型方程初边值问题:
已知其精确解为。分别取步长为和,给出节点处的数值解和误差。
3.1 一阶显格式
代码如下:
#include <cmath>
#include <stdlib.h>
#include <stdio.h>
int main(int argc, char* argv[])
{
int i,j,k,m,n;
double a,h,tau,r,pi;
double *x,*t,**u;
double phi(double x);
double psi(double x);
double alpha(double t);
double beta(double t);
double f(double x, double t);
double exact(double x, double t);
pi=3.14159265359;
m=100;
n=50;
a=1.0;
h=pi/m;
tau=1.0/n;
r=a*tau/h;
r=r*r;
x=(double *)malloc(sizeof(double)*(m+1));
for(i=0;i<=m;i++)
x[i]=i*h;
t=(double *)malloc(sizeof(double)*(n+1));
for(k=0;k<=n;k++)
t[k]=k*tau;
u=(double **)malloc(sizeof(double *)*(m+1));
for(i=0;i<=m;i++)
u[i]=(double *)malloc(sizeof(double)*(n+1));
for(i=0;i<=m;i++)
{
u[i][0]=phi(x[i]);
u[i][1]=phi(x[i])+tau*psi(x[i]);
}
for(k=1;k<=n;k++)
{
u[0][k]=alpha(t[k]);
u[m][k]=beta(t[k]);
}
for(k=1;k<n;k++)
{
for(i=1;i<m;i++)
u[i][k+1]=r*u[i-1][k]+2*(1-r)*u[i][k]+r*u[i+1][k]-u[i][k-1]+tau*tau*f(x[i],t[k]);
}
k=4*n/5;
j=m/10;
for(i=j;i<m;i=i+j)
printf("(x,t)=(%.2f,%.2f),y=%f,error=%.4e.\n",x[i],t[k],u[i][k],fabs(u[i][k]-exact(x[i],t[k])));
free(x);free(t)
for(j=0;j<=m;j++)
free(u[j]);
free(u);
return 0;
}
double phi(double x)
{
return sin(x);
}
double psi(double x)
{
return sin(x);
}
double alpha(double t)
{
return 0.0;
}
double beta(double t)
{
return 0.0;
}
double f(double x, double t)
{
return 2*sin(x)*exp(t);
}
double exact(double x, double t)
{
return sin(x)*exp(t);
}
当时,计算结果如下:
(x,t)=(0.31,0.80),y=0.685504,error=2.2257e-03.
(x,t)=(0.63,0.80),y=1.303907,error=4.2336e-03.
(x,t)=(0.94,0.80),y=1.794673,error=5.8271e-03.
(x,t)=(1.26,0.80),y=2.109765,error=6.8501e-03.
(x,t)=(1.57,0.80),y=2.218338,error=7.2027e-03.
(x,t)=(1.88,0.80),y=2.109765,error=6.8501e-03.
(x,t)=(2.20,0.80),y=1.794673,error=5.8271e-03.
(x,t)=(2.51,0.80),y=1.303907,error=4.2336e-03.
(x,t)=(2.83,0.80),y=0.685504,error=2.2257e-03.
当时,计算结果如下:
(x,t)=(0.31,0.80),y=0.686619,error=1.1106e-03.
(x,t)=(0.63,0.80),y=1.306028,error=2.1124e-03.
(x,t)=(0.94,0.80),y=1.797593,error=2.9075e-03.
(x,t)=(1.26,0.80),y=2.113197,error=3.4180e-03.
(x,t)=(1.57,0.80),y=2.221947,error=3.5939e-03.
(x,t)=(1.88,0.80),y=2.113197,error=3.4180e-03.
(x,t)=(2.20,0.80),y=1.797593,error=2.9075e-03.
(x,t)=(2.51,0.80),y=1.306028,error=2.1124e-03.
(x,t)=(2.83,0.80),y=0.686619,error=1.1106e-03.
计算结果显示是算法是一阶收敛,且精确解及数值解都关于对称。
3.2 二阶显格式
代码如下:
#include <cmath>
#include <stdlib.h>
#include <stdio.h>
int main(int argc, char* argv[])
{
int i,j,k,m,n;
double a,h,tau,r,pi;
double *x,*t,**u;
double phi(double x);
double psi(double x);
double alpha(double t);
double beta(double t);
double f(double x, double t);
double exact(double x, double t);
pi=3.14159265359;
m=100;
n=50;
a=1.0;
h=pi/m;
tau=1.0/n;
r=a*tau/h;
r=r*r;
printf("r=%.4f.\n",r);
x=(double *)malloc(sizeof(double)*(m+1));
for(i=0;i<=m;i++)
x[i]=i*h;
t=(double *)malloc(sizeof(double)*(n+1));
for(k=0;k<=n;k++)
t[k]=k*tau;
u=(double **)malloc(sizeof(double *)*(m+1));
for(i=0;i<=m;i++)
u[i]=(double *)malloc(sizeof(double)*(n+1));
for(i=0;i<=m;i++)
u[i][0]=phi(x[i]);
for(k=1;k<=n;k++)
{
u[0][k]=alpha(t[k]);
u[m][k]=beta(t[k]);
}
for(i=1;i<m;i++)
u[i][1]=(r*u[i-1][0]+2*(1-r)*u[i][0]+r*u[i+1][0]+tau*tau*f(x[i],t[0])+2*tau*psi(x[i]))/2.0;
for(k=1;k<n;k++)
{
for(i=1;i<m;i++)
u[i][k+1]=r*u[i-1][k]+2*(1-r)*u[i][k]+r*u[i+1][k]-u[i][k-1]+tau*tau*f(x[i],t[k]);
}
k=4*n/5;
j=m/10;
for(i=j;i<m;i=i+j)
printf("(x,t)=(%.2f,%.2f),y=%f,error=%.4e.\n",x[i],t[k],u[i][k],fabs(u[i][k]-exact(x[i],t[k])));
free(x);free(t);
for(j=0;j<=m;j++)
free(u[j]);
free(u);
return 0;
}
double phi(double x)
{
return sin(x);
}
double psi(double x)
{
return sin(x);
}
double alpha(double t)
{
return 0.0;
}
double beta(double t)
{
return 0.0;
}
double f(double x, double t)
{
return 2*sin(x)*exp(t);
}
double exact(double x, double t)
{
return sin(x)*exp(t);
}
当时,计算结果如下:
r=0.4053.
(x,t)=(0.31,0.80),y=0.687721,error=8.6480e-06.
(x,t)=(0.63,0.80),y=1.308124,error=1.6450e-05.
(x,t)=(0.94,0.80),y=1.800478,error=2.2641e-05.
(x,t)=(1.26,0.80),y=2.116589,error=2.6616e-05.
(x,t)=(1.57,0.80),y=2.225513,error=2.7986e-05.
(x,t)=(1.88,0.80),y=2.116589,error=2.6616e-05.
(x,t)=(2.20,0.80),y=1.800478,error=2.2641e-05.
(x,t)=(2.51,0.80),y=1.308124,error=1.6450e-05.
(x,t)=(2.83,0.80),y=0.687721,error=8.6480e-06.
当时,计算结果如下:
r=0.4053.
(x,t)=(0.31,0.80),y=0.687728,error=2.1615e-06.
(x,t)=(0.63,0.80),y=1.308136,error=4.1115e-06.
(x,t)=(0.94,0.80),y=1.800495,error=5.6590e-06.
(x,t)=(1.26,0.80),y=2.116609,error=6.6526e-06.
(x,t)=(1.57,0.80),y=2.225534,error=6.9949e-06.
(x,t)=(1.88,0.80),y=2.116609,error=6.6526e-06.
(x,t)=(2.20,0.80),y=1.800495,error=5.6590e-06.
(x,t)=(2.51,0.80),y=1.308136,error=4.1115e-06.
(x,t)=(2.83,0.80),y=0.687728,error=2.1615e-06.
计算结果显示是算法是二阶收敛,且精确解及数值解都关于对称。