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

发布于：2024-04-24 ⋅ 阅读:(23) ⋅ 点赞:(0)

目录

一、研究目标

二、理论推导

2.1 三层显格式建立

2.2 三层显格式改进

三、算例实现

3.1 一阶显格式

3.2 二阶显格式

一、研究目标

        介绍完一阶双曲型偏微分方程的几种差分格式后，我们继续探讨二阶方程的差分格式。这里以非齐次二阶双曲型偏微分方程的初边值问题为研究对象（无法用解析解求解其精确解）：

\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是非零常数。公式（1）的第2式是初始条件，第3式是边界条件。若f(x,t)\equiv 0，则对应的初边值问题与一阶对流方程有类似属性，只不过公式（1）有两条特征线，精确解u(x,t)在点(x,t)处的值的依赖区域为[x-|a|t,x+|a|t]，从而也有相应的CFL条件。

二、理论推导

2.1 三层显格式建立

        差分步骤与之前的抛物型方程差分中介绍的相同。

        第一步：网格剖分。对矩形求解域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)

空间步长为h=1/m,\tau=T/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.

        第三步：偏导数用差商近似。用二阶中心差商近似公式（2）中的二阶偏导数，即

\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{j},t_{k})}\approx \frac{u(x_{j},t_{k+1})-2u(x_{j},t_{k})+u(x_{j},t_{k-1})}{\tau^{2}}

\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k})}\approx \frac{u(x_{j+1},t_{k})-2u(x_{j},t_{k})+u(x_{j-1},t_{k})}{h^{2}}

        初始条件中的一阶偏导数用向前差商近似，即

\frac{\partial u}{\partial t}|_{(x_{j},t_{0})}\approx\frac{u(x_{j},t_{1})-u(x_{j},t_{0})}{\tau}

        将上面三式代入公式（2）中，用数值解代替精确解并忽略高阶项，可得差分格式为

\left\{\begin{matrix} \frac{u^{k+1}_{j}-2u^{k}_{j}+u^{k-1}_{j}}{\tau^{2}}-a^{2}\frac{u^{k}_{j+1}-2u^{k}_{j}+u^{k}_{j-1}}{h^{2}}=f(x_{j},t_{k}),1\leqslant j\leqslant m-1,1\leqslant k\leqslant n-1,\\ u^{0}_{j}=\varphi(x_{j}),\frac{u^{1}_{j}-u^{0}_{j}}{\tau}=\Psi(x_{j}),0\leqslant j\leqslant m,\space\space(3)\\ 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}，可以将公式（3）整理为公式（4）所示的三层显格式：

\left\{\begin{matrix} u^{k+1}_{j}=ru^{k}_{j-1}+2(1-r)u^{k}_{j}+ru^{k}_{j+1}-u^{k-1}_{j}+\tau^{2}f(x_{j},t_{k}),1\leqslant j\leqslant m-1,1\leqslant k\geqslant n-1,\\ u^{0}_{j}=\varphi(x_{j}),\frac{u^{1}_{j}-u^{0}_{j}}{\tau}=\Psi(x_{j}),0\leqslant j\leqslant m,\space\space(4)\\ u^{k}_{0}=\alpha(t_{k}),u^{k}_{m}=\beta(t_{k}),1\leqslant k\leqslant n \end{matrix}\right.

公式（4）的局部截断误差为O(\tau+h^{2})。

2.2 三层显格式改进

        事实上，公式（4）的精度比较低，原因是在对初始条件进行一阶偏导数近似时，采用的是精度较低的向前差商。为了提高精度，可以考虑利用中心差商对时间的一阶偏导数进行近似，即

\frac{\partial u}{\partial t}|_{(x_{j},t_{0})}\approx\frac{u(x_{j},t_{1})-u(x_{j},t_{-1})}{2\tau}

        得到公式（1）的数值格式为

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

        这样处理引入虚拟越界点u^{-1}_{j}，如果认为公式（5）中第1式对k=0也成立，即

u^{1}_{j}=ru^{0}_{j-1}+2(1-r)u^{0}_{j}+ru^{k}_{j+1}-u^{-1}_{j}+\tau^{2}f(x_{j},t_{k}),1\leqslant j\leqslant m-1

得出    u^{-1}_{j}=ru^{0}_{j-1}+2(1-r)u^{0}_{j}+ru^{0}_{j+1}+\tau^{2}f(x_{j},t_{0}),1\leqslant j\leqslant m-1

于是可以在公式（5）第2式中消去u^{-1}_{j}，得到

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

从而得到改进的三层格式：

\left\{\begin{matrix} u^{k+1}_{j}=ru^{k}_{j-1}+2(1-r)u^{k}_{j}+ru^{k}_{j+1}-u^{k-1}_{j}+\tau^{2}f(x_{j},t_{k}),1\leqslant j\leqslant m-1,1\leqslant k\leqslant n-1,\\ u^{0}_{j}=\varphi(x_{j}),0\leqslant i\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.

        该格式的局部截断误差为O(\tau^{2}+h^{2})。

三、算例实现

        计算双曲型方程初边值问题：

\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/100和\tau_{1}=1/100,h_{1}=\pi/200，给出节点(\frac{i\pi}{10},\frac{4}{5}),i=1,\cdot\cdot\cdot,9处的数值解和误差。

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);

}

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

(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.

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

(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.

        计算结果显示是算法是一阶收敛，且精确解及数值解都关于x=\pi/2对称。 

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);

}

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

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.

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

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.

          计算结果显示是算法是二阶收敛，且精确解及数值解都关于x=\pi/2对称。

网站公告

今日签到

点亮在社区的每一天
去签到

热门文章

最新发布