偏微分方程算法之二阶双曲型方程紧差分方法

一、研究目标

二、理论推导

三、算例实现

一、研究目标

前面我们已经介绍了二阶双曲型方程显式、隐式差分格式，可否像抛物型方程一样，构建更高精度的差分格式。接下来我们介绍紧差分格式。这里继续以非齐次二阶双曲型偏微分方程的初边值问题为研究对象：

$\left\{\begin{matrix} \frac{\partial^{2}u(x,t)}{\partial t^{2}}-a^{2}\frac{\partial^{2}u(x,t)}{\partial x^{2}}=f(x,t),0<x<1,0<t\leqslant T,\\ u(x,0)=\varphi(x),\frac{\partial u}{\partial t}(x,0)=\Psi(x),0\leqslant x\leqslant 1,\space\space(1)\\ u(0,t)=\alpha(t),u(1,t)=\beta(t),0<t\leqslant T \end{matrix}\right.$

公式（1）中u表示一个与时间t和位置x有关的待求波函数， $\varphi(x),\Psi(x),\alpha(t),\beta(t)$ 及方程右端项函数 $f(x,t)$ 都是已知函数， $a,T$ 是非零常数。

二、理论推导

第一步：网格剖分。对矩形求解域 $0\leqslant x\leqslant 1,0\leqslant t\leqslant T$ 进行等距剖分，即

$x_{j}=jh(j=0,1,\cdot\cdot\cdot,m),t_{k}=k\tau(k=0,1,\cdot\cdot\cdot,n)$

第二步：弱化原方程。将原来的连续方程离散到网格节点上成立，得到离散方程：

$\left\{\begin{matrix} \frac{\partial^{2}u}{\partial t^{2}}|_{(x_{j},t_{k})}-a^{2}\frac{\partial ^{2}u}{\partial x^{2}}|_{(x_{j},t_{k})}=f(x_{j},t_{k}),0<j<m,0<k\leqslant n,\\ u(x_{j},t_{0})=\varphi(x_{j}),\frac{\partial u}{\partial t}(x_{j},t_{0})=\Psi(x_{j}),0\leqslant j\leqslant m,\space\space(2)\\ u(x_{0},t_{k})=\alpha(t_{k}),u(x_{m},t_{k})=\beta(t_{k}),1\leqslant k\leqslant n \end{matrix}\right.$

第三步：对偏导数进行更高精度近似。由泰勒公式（固定时间t不变）：

$u(x_{j-1},t_{k})=(u-h\frac{\partial u}{\partial x}+\frac{h^{2}}{2}\frac{\partial ^{2}u}{\partial x^{2}}-\frac{h^{3}}{6}\frac{\partial^{3}u}{\partial x^{3}}+\frac{h^{4}}{24}\frac{\partial^{4}u}{\partial x^{4}}-\frac{h^{5}}{120}\frac{\partial^{5}u}{\partial x^{5}})|_{(x_{j},t_{k})}+O(h^{6})$

$u(x_{j+1},t_{k})=(u+h\frac{\partial u}{\partial x}+\frac{h^{2}}{2}\frac{\partial ^{2}u}{\partial x^{2}}+\frac{h^{3}}{6}\frac{\partial^{3}u}{\partial x^{3}}+\frac{h^{4}}{24}\frac{\partial^{4}u}{\partial x^{4}}+\frac{h^{5}}{120}\frac{\partial^{5}u}{\partial x^{5}})|_{(x_{j},t_{k})}+O(h^{6})$

将上面两式相加可得

$u(x_{j+1},t_{k})-2u(x_{j},t_{k})+u(x_{j-1},t_{k})=h^{2}\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k})}+\frac{h^{4}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k})}+O(h^{6})$

有 $\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k})}=\frac{u(x_{j+1},t_{k})-2u(x_{j},t_{k})+u(x_{j-1},t_{k})}{h^{2}}-\frac{h^{2}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k})}+O(h^{4})$

类似的有

$\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k-1})}=\frac{u(x_{j+1},t_{k-1})-2u(x_{j},t_{k-1})+u(x_{j-1},t_{k-1})}{h^{2}}-\frac{h^{2}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k-1})}+O(h^{4})$

$\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k+1})}=\frac{u(x_{j+1},t_{k+1})-2u(x_{j},t_{k+1})+u(x_{j-1},t_{k+1})}{h^{2}}-\frac{h^{2}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k+1})}+O(h^{4})$

将上面两式相加后除以2得

$\frac{1}{2}(\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k-1})}+\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k+1})})=\frac{u(x_{j+1},t_{k-1})-2u(x_{j},t_{k-1})+u(x_{j-1},t_{k+1})}{2h^{2}}+\frac{u(x_{j+1},t_{k+1})-2u(x_{j},t_{k+1})+u(x_{j-1},t_{k+1})}{2h{^{2}}}-\frac{h^{2}}{12}\frac{1}{2}(\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k-1})}+\frac{\partial ^{4}u}{\partial x^{4}}|_{(x_{j},t_{k+1})})+O(h^{4})$

从而 $\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k})}+O(\tau^{2})=\frac{u(x_{j+1},t_{k-1})-2u(x_{j},t_{k-1})+u(x_{j-1},t_{k-1})}{2h^{2}}+\frac{u(x_{j+1},t_{k+1})-2u(x_{j},t_{_{k+1}})+u(x_{j-1},t_{k+1})}{2h^{2}}-\frac{h^{2}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k})}+O(h^{4}+\tau^{2}h^{2}) \space\space(3)$

现令

$D=u(x_{j+1},t_{k-1})-2u(x_{j},t_{k-1})+u(x_{j-1},t_{k-1})+u(x_{j+1},t_{k+1})-2u(x_{j},t_{k+1})+u(x_{j-1},t_{k+1})$ 及 $\frac{\partial^{2}u}{\partial x^{2}}=v(x,t)$

则公式（3）可写作

$v(x_{j},t_{k})=\frac{D}{2h^{2}}-\frac{h^{2}}{12}\frac{\partial^{2}v}{\partial x^{2}}|_{(x_{j},t_{k})}+O(\tau^{2}h^{2}+h^{4}+\tau^{2})$

$=\frac{D}{2h^{2}}-\frac{h^{2}}{12}\frac{v(x_{j+1},t_{k})-2v(x_{j},t_{k})+v(x_{j-1},t_{k})}{h^{2}}+O(\tau^{2}h^{2}+h^{4}+\tau^{2})$

整理可得

$\frac{v(x_{j+1},t_{k})+10v(x_{j},t_{k})+v(x_{j-1},t_{k})}{12}=\frac{D}{2h^{2}}+O(\tau^{2}h^{2}+h^{4}+\tau^{2})\space\space\space\space(4)$

由于原双曲型方程为 $\frac{\partial^{2}u}{\partial t^{2}}-a^{2}v(x,t)=f(x,t)$ ，也即 $v(x,t)=\frac{1}{a^{2}}(\frac{\partial^{2}u}{\partial t^{2}}-f(x,t))$

公式（4）可改写为

$\frac{1}{12a^{2}}(\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{j-1},t_{k})}-f(x_{j-1},t_{k})+10\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{j},t_{k})}-10f(x_{j},t_{k})+\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{j+1},t_{k})}-f(x_{j+1},t_{k}))=\frac{D}{2h^{2}}+O(\tau^{2}h^{2}+h^{4}+\tau^{2})$

再利用中心差商近似

$\frac{1}{12a^{2}}(\frac{u(x_{j-1},t_{k+1})-2u(x_{j-1},t_{k})+u(x_{j-1},t_{k-1})}{\tau^{2}}-f(x_{j-1},t_{k})+10\frac{u(x_{j},t_{k+1})-2u(x_{j},t_{k})+u(x_{j},t_{k-1})}{\tau^{2}}-10f(x_{j},t_{k})+\frac{u(x_{j+1},t_{k+1})-2u(x_{j+1},t_{k}+u(x_{j+1},t_{k-1}))}{\tau^{2}}-f(x_{j+1},t_{k}))=\frac{D}{2h^{2}}+O(\tau^{2}h^{2}+h^{4}+\tau^{2})$

上式中用数值解代替精确解并忽略高阶项，可得

$\frac{1}{12a^{2}\tau^{2}}(u^{k+1}_{j-1}-2u^{k}_{j-1}+u^{k+1}_{j-1}+10(u^{k+1}_{j}-2u^{k}_{j}+u^{k-1}_{j})+u^{k+1}_{j+1}-2u^{k}_{j+1}+u^{k-1}_{j+1})=\frac{1}{2}h^{2}(u^{k-1}_{j+1}-2u^{k-1}_{j}+u^{k-1}_{j-1}+u^{k+1}_{j-1}-2u^{k+1}_{j}+u^{k+1}_{j-1})+\frac{1}{12a^{2}}(f^{k}_{j+1}+10f^{k}_{j}+f^{k}_{j-1})$

其中， $f^{k}_{l}=f(x_{l},t_{k}),l=j-1,j,j+1$ 。

联合初边值条件，可得到以下紧差分格式：

$\left\{\begin{matrix} (6r-1)u^{k+1}_{j-1}+(10+12r)u^{k+1}_{j}+(1-6r)u^{k+1}_{j+1}=(6r-1)u^{k-1}_{j-1}-(12r+10)u^{k-1}_{j}+(6r-1)u^{k-1}_{j+1}+\\2(u^{k}_{j-1}+10u^{k}_{j}+u^{k}_{j+1})+\tau^{2}(f^{k}_{j-1}+10f^{k}_{j}+f^{k}_{j+1}),1\leqslant i\leqslant m-1,1\leqslant k\leqslant n-1,\\ u^{0}_{j}=\varphi(_{j}),0\leqslant j\leqslant m,\\ u^{1}_{j}=(ru^{0}_{j-1}+2(1-r)u^{0}_{j}+ru^{0}_{j+1}+\tau^{2}f(x_{j},t_{0})+2\tau\Psi(x_{j}))/2,1\leqslant j\leqslant m-1,\\ u^{k}_{0}=\alpha(t_{k}),u^{k}_{m}=\beta(t_{k}),1\leqslant k\leqslant n \end{matrix}\right.$

其中 $r=a^{2}\tau^{2}/h^{2}>0$ ，局部截断误差为 $O(\tau^{2}+h^{4})$ ，关于时间t是二阶的，关于空间x是四阶的。

三、算例实现

紧差分格式计算双曲型偏微分方程初边值问题：

$\left\{\begin{matrix} \frac{\partial^{2}u(x,t)}{\partial t^{2}}-\frac{\partial^{2}u(x,t)}{\partial x^{2}}=2e^{t}sinx,0<x<\pi,0<t\leqslant 1,\\ u(x,0)=sinx,\frac{\partial u}{\partial t}(x,0)=sinx,0\leqslant x\leqslant \pi,\\ u(0,t)=0,u(\pi,t)=0,0<t\leqslant 1 \end{matrix}\right.$

已知其精确解为 $u(x,t)=e^{t}sinx$ 。分别取步长为 $\tau_{1}=1/50,h_{1}=\pi/200$ 和 $\tau_{1}=1/100,h_{1}=\pi/400$ ，给出节点 $(\frac{i\pi}{10},\frac{4}{5}),i=1,\cdot\cdot\cdot,5$ 处的数值解和误差。

代码如下：

#include <cmath>
#include <stdlib.h>
#include <stdio.h>


int main(int argc, char* argv[])
{
        int i,j,k,m,n;
        double a,h,r,tau,pi,c1,c2;
        double *x,*t,**u,*a1,*b,*c,*d,*ans;
        double phi(double x);
        double ddphi(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);
        double *chase_algorithm(double *a, double *b, double *c, double *d, int n);

        m=200;
        n=50;
        a=1.0;
        pi=3.14159265359;
        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;

        a1=(double*)malloc(sizeof(double)*(m-1));
        b=(double*)malloc(sizeof(double)*(m-1));
        c=(double*)malloc(sizeof(double)*(m-1));
        d=(double*)malloc(sizeof(double)*(m-1));
        ans=(double*)malloc(sizeof(double)*(m-1));
        c1=1.0-6*r;
        c2=10.0+12*r;

        for(k=1;k<n;k++)
        {
                for(i=1;i<m;i++)
                {
                        d[i-1]=(-c1)*(u[i-1][k-1]+u[i+1][k-1])-c2*u[i][k-1]+2*(u[i-1][k]+10*u[i][k]+u[i+1][k])+tau*tau*(f(x[i-1],t[k])+10*f(x[i],t[k])+f(x[i+1],t[k]));
                        a1[i-1]=c1;
                        b[i-1]=c2;
                        c[i-1]=a1[i-1];
                }
                d[0]=d[0]-c1*u[0][k+1];
                d[m-2]=d[m-2]-c1*u[m][k+1];
                ans=chase_algorithm(a1,b,c,d,m-1);
                for(i=0;i<m-1;i++)
                        u[i+1][k+1]=ans[i];
        }

        free(ans);
        k=4*n/5;
        j=m/10;
        for(i=j;i<=m/2;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(a1);free(b);free(c);free(d);free(x);free(t);

        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);
}
double *chase_algorithm(double *a, double *b, double *c, double *d, int n)
{
        int i;
        double *ans,*g,*w,p;

        ans=(double*)malloc(sizeof(double)*n);
        g=(double*)malloc(sizeof(double)*n);
        w=(double*)malloc(sizeof(double)*n);;
        g[0]=d[0]/b[0];
        w[0]=c[0]/b[0];

        for(i=1;i<n;i++)
        {
                p=b[i]-a[i]*w[i-1];
                g[i]=(d[i]-a[i]*g[i-1])/p;
                w[i]=c[i]/p;
        }

        ans[n-1]=g[n-1];
        i=n-2;
        do
        {
                ans[i]=g[i]-w[i]*ans[i+1];
                i=i-1;
        }while(i>=0);

        free(g);free(w);
        return ans;
}

当 $\tau_{1}=1/50,h_{1}=\pi/200$ 时，计算结果如下：

r=1.6211.
(x,t)=(0.31,0.80),y=0.687686,error=4.3538e-05.
(x,t)=(0.63,0.80),y=1.308057,error=8.2815e-05.
(x,t)=(0.94,0.80),y=1.800386,error=1.1398e-04.
(x,t)=(1.26,0.80),y=2.116481,error=1.3400e-04.
(x,t)=(1.57,0.80),y=2.225400,error=1.4089e-04.

当 $\tau_{1}=1/100,h_{1}=\pi/400$ 时，计算结果如下：

r=0.4053.
(x,t)=(0.31,0.80),y=0.687727,error=2.7423e-06.
(x,t)=(0.63,0.80),y=1.308135,error=5.2162e-06.
(x,t)=(0.94,0.80),y=1.800493,error=7.1794e-06.
(x,t)=(1.26,0.80),y=2.116607,error=8.4399e-06.
(x,t)=(1.57,0.80),y=2.225532,error=8.8743e-06.

从计算结果可知，当时间步长减小为原来的1/4、空间步长减小为原来的1/2时，误差减小为原来的1/16。