紧致有限差分方法求解全离散波动方程

Compact Finite Difference Method for Solving Fully Discrete Wave Equations

摘要: 本文针对整数阶波动方程，给出了一种基于紧致有限差分方法的隐式全离散格式。该格式在时间方向采用中心差分格式来离散，在空间方向采用紧致中心差商的权平均来离散。离散格式的稳定性分析及误差估计表明，该离散格式在时间方向达到二阶收敛，空间方向达到四阶收敛。并且通过数值实验证明该离散格式的收敛阶为O(τ2h4)。

Abstract: This article proposes an implicit fully discrete scheme based on the compact finite difference method for integer order wave equations. This discretization scheme uses a central difference scheme for discretization in the temporal direction and a compact central difference quotient weighted average for discretization in the spatial direction. The stability analysis and error estimation of the discrete format indicate that it achieves second-order convergence in the temporal direction and fourth-order convergence in the spatial direction. And numerical experiments have shown that the convergence order of the discrete format isO(τ2h4).

1. 引言

在过去的许多年中，偏微分方程发数值解法一直备受学者们的关注。波动方程作为偏微分方程中的一类，其精确解很难直接求出，通过各种数值方法求解波动方程的精确解备受关注。因此，求解波动方程的初边值问题是偏微分方程数值解研究方向的热门问题。例如，Du等人为了研究分数阶扩散波动方程，提出了一种高阶的差分方法 [1] 。安文静等人用有限差分方法给出了波动方程的一种离散格式，并且用Fourier方法分析了该格式的无条件稳定性，用能量不等式法分析了误差估计，理论分析表明他们提出的格式收敛阶为 $O\left({\tau }^2+h^4\right)$ [2] 。廖洪林利用离散H²能量法给出三层显式格式和交替方向隐式格式的最大模误差估计，并且得到时间和空间方向的收敛阶都为二阶精度 [3] 。朱爱玲等人针对波动方程推导出了一种显式差分外推格式，将收敛精度由原来的二阶提高为四阶 [4] 。由于传统的有限差分方法不能达到高阶的差分近似，紧致有限差分方法被许多学者广泛的关注。紧致差分方法之所以被许多学者们关注，是因为它可以使用比较少的网格点来达到高精度的收敛阶 [5] 。关于紧致有限差分格式的详细介绍可以参见 [6] 。例如，Gao和Sun对于满足Neumann边界条件的分数阶次扩散方程和热方程提出了一种紧致有限差分格式 [7] [8] 。Ren等人提出了一种紧致有限差分格式去求解具有Neumann边界条件的时间分数阶次扩散方程， [9] 。Li等人为了解决地下水污染问题，针对空间二维对流扩散方程研究了一种四阶紧致有限差分方法 [10] 。Liao等人提出了一种紧致算法去解决非线性反应扩散方程，并且时间和空间方向的收敛精度都达到了四阶 [11] 。

本文所研究的方程如下：

$\{\begin{array}{l}\frac{{\partial }^{2}u}{\partial t^2}=a^2\frac{{\partial }^{2}u}{\partial x^2}+f\left(x,t\right),x\in \left[0,L\right],t\in \left[0,T\right],\\ u\left(x,0\right)=\phi \left(x\right),\frac{\partial u}{\partial t}\left(x,0\right)=\psi \left(x\right),x\in \left[0,L\right],\\ u\left(0,t\right)=0,u\left(L,t\right)=0,t\in \left[0,T\right].\end{array}$ (1)

其中， $f\left(x,t\right),u\left(x,t\right),\phi \left(x\right),\psi \left(x\right)$ 是足够光滑的函数， $L,T$ 是已知的正常数。

本篇文章的章节内容安排如下。在第二节中，我们对方程(1)进行紧致有限差分格式的构造。在第三节中，我们用能量不等式方法分析了离散格式的稳定性和误差估计。第四节是本文的数值实验部分，第五节是本文的总结。

2. 离散格式的构造

用有限差分方法求解波动方程的时候，必须把连续问题进行离散化。在对方程(1)进行离散格式构造之前，我们需要对空间区间 $\left[0,L\right]$ 和时间区间 $\left[0,T\right]$ 进行网格剖分。本篇文章中，对于给定的正整数M和N，记 $x_j=jh\left(j=0,1,2,\cdots ,M\right)$ ， $t_n=n\tau \left(n=0,1,2,\cdots ,N\right)$ 其中 $h=L/M$ ， $\tau =T/N$ 。分别称h和 $\tau$ 为空间和时间步长。

空间网格函数的定义如下：

$V_h=\left\{v|v=\left(v_0,v_1,\cdots ,v_M\right),v_0=v_M=0\right\}.$

紧致差分算子 $\mathcal{H}$ 的定义如下：

$\{\begin{array}{l}\frac{1}{12}\left(u_{j-1}^n+10u_j^n+u_{j+1}^n\right)=\left(1+\frac{h^2}{12}{\delta }_x^2\right)u_j^n,1\le j\le M-1,\\ u_j^n,j=0,M.\end{array}$ (2)

下面提出的符号和一些引理在后面的证明中将会用到

$\begin{array}{l}{\delta }_x^2{u}_j^n=\frac{1}{h^2}\left(u_{j+1}^n-2u_j^n+u_{j-1}^n\right),\langle u^n,v^n\rangle =h{\displaystyle \underset{j=1}{\overset{M-1}{\sum }}{u}_j^n}{v}_j^n,\\ \Vert u^n\Vert =\sqrt{\langle u^n,u^n\rangle },\Vert {\delta }_x^2{u}^n\Vert =\sqrt{\langle {\delta }_x^2{u}^n,{\delta }_x^2{u}^n\rangle },\\ \langle {\delta }_x^2{u}^k,u^n\rangle =-\langle {\delta }_x{u}^k,{\delta }_x{u}^n\rangle ,{\Vert u\Vert }_{\mathcal{H}}=\sqrt{{\left(u,u\right)}_{\mathcal{H}}},\\ {\left(u,v\right)}_{\mathcal{H}}=\langle {\delta }_{x}u,{\delta }_{x}v\rangle -\frac{h^2}{12}\langle {\delta }_x^{2}u,{\delta }_x^{2}v\rangle .\end{array}$

引理2.1. [12] 若 $p\left(x\right)\in C^6\left[0,L\right]$ ，并且 $x_i=ih$ ， $0\le i\le M$ 则

$\frac{1}{12}\left[p^{″}\left(x_{i+1}\right)+10p^{″}\left(x_i\right)+p^{″}\left(x_{i-1}\right)\right]-\frac{1}{h^2}\left[p\left(x_{i+1}\right)+10p\left(x_i\right)+p\left(x_{i-1}\right)\right]=\frac{h^4}{240}{p}^{\left(6\right)}(\; \xi \; i\; )$

其中， ${\xi }_i\in \left(x_{i-1},x_{i+1}\right)$ ， $1\le i\le M-1$ 。

对于离散格式的构造方面。时间方向采用中心有限差分法离散，空间方向在n-1层、n层和n+1层上用中心差商的权平均去离散，即

$\begin{array}{l}\frac{{\partial }^{2}u}{\partial t^2}=\frac{u_j^{n+1}-2u_j^n+u_j^{n-1}}{{\tau }^2}+O\left({\tau }^2\right),\\ \frac{{\partial }^{2}u}{\partial x^2}=\frac{1}{\text{12}}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_{n+1}\right)+\frac{1\text{0}}{\text{1}2}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_n\right)+\frac{1}{\text{12}}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_{n-1}\right).\end{array}$

因此方程(1)可以先离散为

$\frac{u_j^{n+1}-2u_j^n+u_j^{n-1}}{{\tau }^2}=\frac{a^2}{\text{12}}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_{n+1}\right)+\frac{\text{10}{a}^2}{\text{12}}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_n\right)+\frac{a^2}{\text{12}}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_{n-1}\right)+f_j^n+O\left({\tau }^2\right),$ (3)

接下来将 $\mathcal{H}$ 算子作用在(3)式的左右两端，我们可以得到

$\mathcal{H}\left(\frac{u_j^{n+1}-2u_j^n+u_j^{n-1}}{{\tau }^2}\right)=\frac{a^2}{\text{12}}\mathcal{H}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_{n+1}\right)+\frac{\text{10}{a}^2}{\text{1}2}\mathcal{H}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_n\right)+\frac{a^2}{\text{12}}\mathcal{H}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_{n-1}\right)+\mathcal{H}{f}_j^n+O\left({\tau }^2\right).$ (4)

由引理2.1，有

$\mathcal{H}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_{n+1}\right)={\delta }_x^2{u}_j^{n+1}+O\left(h^4\right),$

$\mathcal{H}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_n\right)={\delta }_x^2{u}_j^n+O\left(h^4\right),$

$\mathcal{H}\frac{{\partial }^{2}u}{\partial x^2}\left(x_j,t_{n-1}\right)={\delta }_x^2{u}_j^{n-1}+O\left(h^4\right).$

将上面三个离散的式子代入到(4)式中，可得到

$\{\begin{array}{l}\mathcal{H}\left(\frac{u_j^{n+1}-2u_j^n+u_j^{n-1}}{{\tau }^2}\right)=\frac{a^2}{\text{12}}{\delta }_x^2{u}_j^{n+1}+\frac{\text{10}{a}^2}{\text{1}2}{\delta }_x^2{u}_j^n+\frac{a^2}{\text{12}}{\delta }_x^2{u}_j^{n-1}+\mathcal{H}{f}_j^n+R_j^n,\\ 1\le j\le M-1,1\le n\le N-1,\\ u_j^0=\phi \left(x_j\right),1\le j\le M-1,\\ u_0^n=0,u_M^n=0,0\le n\le N,\end{array}$ (5)

在上式中，误差项 $R_j^n=O\left({\tau }^2+h^4\right)$ 。舍去截断误差项 $R_j^n$ 后，用 $U_j^n$ 代替精确解 $u_j^n$ ，可得到方程(1)的紧致有限差分格式如下

$\{\begin{array}{l}\mathcal{H}\left(\frac{U_j^{n+1}-2U_j^n+U_j^{n-1}}{{\tau }^2}\right)=\frac{a^2}{\text{12}}{\delta }_x^2{U}_j^{n+1}+\frac{\text{10}{a}^2}{\text{1}2}{\delta }_x^2{U}_j^n+\frac{a^2}{\text{12}}{\delta }_x^2{U}_j^{n-1}+\mathcal{H}{f}_j^n,\\ 1\le j\le M-1,1\le n\le N-1,\\ U_j^0=\phi \left(x_j\right),1\le j\le M-1,\\ U_0^n=0,U_M^n=0,0\le n\le N.\end{array}$ (6)

3. 稳定性分析和误差估计

接下来进行离散格式(6)的稳定性分析，然后对产生的误差进行估计。在接下来的分析中，常数C的取值随着使用位置的不同而变化。

第一步进行离散格式(6)的稳定性分析。在接下来的证明中，我们省略下标j。

引理3.5. 假设 $u,v\in V_h$ ，则

$-\langle {\delta }_x^{2}u,\mathcal{H}v\rangle ={\left(u,v\right)}_{\mathcal{H}}$ 。

证明. 根据(2)式中 $\mathcal{H}$ 的定义有

$-\langle {\delta }_x^{2}u,\mathcal{H}v\rangle =-\langle {\delta }_x^{2}u,\left(1+\frac{h^2}{12}{\delta }_x^2\right)v\rangle =-\langle {\delta }_x^{2}u,v\rangle -\frac{h^2}{12}\langle {\delta }_x^{2}u,{\delta }_x^{2}v\rangle$ ，

又根据 $\langle {\delta }_x^2{u}^k,u^n\rangle =-\langle {\delta }_x{u}^k,{\delta }_x{u}^n\rangle$ ，故上式可改写为

$-\langle {\delta }_x^{2}u,\mathcal{H}v\rangle =\langle {\delta }_{x}u,{\delta }_{x}v\rangle -\frac{h^2}{12}\langle {\delta }_x^{2}u,{\delta }_x^{2}v\rangle ={\left(u,v\right)}_{\mathcal{H}}$ ，

引理证明完毕。 □

引理3.6. [13] 令 $v\in V_h$ ，那么有如下不等式成立，

$\frac{1}{3}{\Vert v\Vert }^2\le {\Vert \mathcal{H}v\Vert }^2\le {\Vert v\Vert }^2.$

引理3.7. [14] 令 $u^{-1}=u^0-\tau \psi$ 和 ${\epsilon }^{-1}=u\left(x,t_{-1}\right)-u^{-1}$ ，那么有

$\left|{\epsilon }^{-1}\right|\le C{\tau }^2,$

在上式中，C是一个常数。

以上引理在接下来的证明中将会用到。

定理3.8. 令 $U^n$ 是差分格式(6)的数值解。如果有下面的不等式成立

$\Vert U^n\Vert \le C\left(\Vert U^0\Vert +\underset{0\le s\le n-1}{\max }\Vert f^s\Vert \right),$

则差分格式(6)是无条件稳定的，其中，C是一个正常数。

证明. 通过(6)式，我们有

$\mathcal{H}{\mathcal{U}}^{n+1}-\frac{a^2{\tau }^2}{12}{\delta }_x^2{U}^{n+1}=2\mathcal{H}{\mathcal{U}}_j^n-\mathcal{H}{\mathcal{U}}^{n-1}+\frac{10a^2{\tau }^2}{12}{\delta }_x^2{U}^n+\frac{a^2{\tau }^2}{12}{\delta }_x^2{U}^{n-1}+\mathcal{H}{f}^n,$ (7)

用 $\mathcal{H}{\mathcal{U}}_j^{n+1}$ 在(7)式的两边做内积，得到

$\begin{array}{l}\langle \mathcal{H}{\mathcal{U}}^{n+1},\mathcal{H}{\mathcal{U}}^{n+1}\rangle -\frac{a^2{\tau }^2}{12}\langle {\delta }_x^2{U}^{n+1},\mathcal{H}{\mathcal{U}}^{n+1}\rangle \\ =2\langle \mathcal{H}{\mathcal{U}}^n,\mathcal{H}{\mathcal{U}}^{n+1}\rangle -\langle \mathcal{H}{\mathcal{U}}^{n-1},\mathcal{H}{\mathcal{U}}^{n+1}\rangle +\frac{10a^2{\tau }^2}{12}\langle {\delta }_x^2{U}^n,\mathcal{H}{\mathcal{U}}^{n+1}\rangle \\ +\frac{a^2{\tau }^2}{12}\langle {\delta }_x^2{U}^{n-1},\mathcal{H}{\mathcal{U}}^{n+1}\rangle +{\tau }^2\langle \mathcal{H}{f}^n,\mathcal{H}{\mathcal{U}}^{n+1}\rangle .\end{array}$

我们以引理3.5为依据，可得

$\begin{array}{l}\langle \mathcal{H}{\mathcal{U}}^{n+1},\mathcal{H}{\mathcal{U}}^{n+1}\rangle \le 2\langle \mathcal{H}{\mathcal{U}}^n,\mathcal{H}{\mathcal{U}}^{n+1}\rangle -\langle \mathcal{H}{\mathcal{U}}^{n-1},\mathcal{H}{\mathcal{U}}^{n+1}\rangle -\frac{10a^2{\tau }^2}{12}{\left(U^n,U^{n+1}\right)}_{\mathcal{H}}\\ -\frac{a^2{\tau }^2}{12}{\left(U^{n-1},U^{n+1}\right)}_{\mathcal{H}}+{\tau }^2\langle \mathcal{H}{f}^n,\mathcal{H}{\mathcal{U}}^{n+1}\rangle .\end{array}$ (8)

对于(8)式，根据范数等价定理，Cauchy-Schwarz不等式以及引理3.6得到

$\begin{array}{l}{\Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert }^2\le 2\Vert \mathcal{H}{\mathcal{U}}^n\Vert \Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert -\Vert \mathcal{H}{\mathcal{U}}^{n-1}\Vert \Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert +\frac{10a^2{\tau }^2}{12}{\Vert {\mathcal{U}}^n\Vert }_{\mathcal{H}}{\Vert {\mathcal{U}}^{n+1}\Vert }_{\mathcal{H}}\\ +\frac{a^2{\tau }^2}{12}{\Vert {\mathcal{U}}^{n-1}\Vert }_{\mathcal{H}}{\Vert {\mathcal{U}}^{n+1}\Vert }_{\mathcal{H}}+{\tau }^2\Vert \mathcal{H}{f}^n\Vert \Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert \\ \le 2\Vert \mathcal{H}{\mathcal{U}}^n\Vert \Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert -\Vert \mathcal{H}{\mathcal{U}}^{n-1}\Vert \Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert +\frac{C_{1}10a^2{\tau }^2}{12}\Vert {\mathcal{U}}^n\Vert \Vert {\mathcal{U}}^{n+1}\Vert \\ +\frac{C_2{a}^2{\tau }^2}{12}\Vert {\mathcal{U}}^{n-1}\Vert \Vert {\mathcal{U}}^{n+1}\Vert +{\tau }^2\Vert \mathcal{H}{f}^n\Vert \Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert \\ \le 2\Vert \mathcal{H}{\mathcal{U}}^n\Vert \Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert -\Vert \mathcal{H}{\mathcal{U}}^{n-1}\Vert \Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert +\frac{30C_1{a}^2{\tau }^2}{12}\Vert \mathcal{H}{\mathcal{U}}^n\Vert \Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert \\ +\frac{3C_2{a}^2{\tau }^2}{12}\Vert \mathcal{H}{\mathcal{U}}^{n-1}\Vert \Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert +{\tau }^2\Vert \mathcal{H}{f}^n\Vert \Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert .\end{array}$

对于上式，可写为

$\begin{array}{l}\Vert \mathcal{H}{\mathcal{U}}^{n+1}\Vert \le \left(2+\frac{30C_1{a}^2{\tau }^2}{12}\right)\Vert \mathcal{H}{\mathcal{U}}^n\Vert +\left(1+\frac{3C_2{a}^2{\tau }^2}{12}\right)\Vert \mathcal{H}{\mathcal{U}}^{n-1}\Vert +{\tau }^2\Vert \mathcal{H}{f}^n\Vert \\ \le C\left(\Vert \mathcal{H}{\mathcal{U}}^n\Vert +\Vert \mathcal{H}{\mathcal{U}}^{n-1}\Vert +\Vert \mathcal{H}{f}^n\Vert \right).\end{array}$ (9)

接下来将通过归纳法证明。当 $n=0$ 时，根据(9)式和引理3.7，我们有

$\begin{array}{l}\Vert \mathcal{H}{\mathcal{U}}^1\Vert \le C\left(\Vert \mathcal{H}{U}^0\Vert +\Vert \mathcal{H}{U}^{-1}\Vert +\Vert \mathcal{H}{f}^0\Vert \right)\\ \le C\left(\Vert \mathcal{H}{U}^0\Vert +\frac{3}{2}\Vert \mathcal{H}{U}^0\Vert +\frac{1}{2}\Vert \mathcal{H}{U}^1\Vert +\Vert \mathcal{H}{f}^0\Vert \right)\\ \le C\left(\Vert \mathcal{H}{U}^0\Vert +\Vert \mathcal{H}{f}^0\Vert \right).\end{array}$ (10)

假设(10)式对于 $m=1,2,\cdots ,n-1$ 时上述结论也成立，即

$\Vert \mathcal{H}{\mathcal{U}}^m\Vert \le C\left(\Vert \mathcal{H}{\mathcal{U}}^0\Vert +\underset{1\le s\le n-1}{\max }\Vert \mathcal{H}{f}^s\Vert \right).$ (11)

现在，我们证明当 $m=n$ 时结论也成立。

$\begin{array}{l}\Vert \mathcal{H}{\mathcal{U}}^n\Vert \le C\left(\Vert \mathcal{H}{\mathcal{U}}^{n-1}\Vert +\Vert \mathcal{H}{\mathcal{U}}^{n-2}\Vert +\Vert \mathcal{H}{f}^{n-1}\Vert \right)\\ \le C\left(\Vert \mathcal{H}{\mathcal{U}}^0\Vert +\underset{1\le s\le n-1}{\max }\Vert \mathcal{H}{f}^s\Vert +\Vert \mathcal{H}{\mathcal{U}}^0\Vert +\underset{1\le s\le n-1}{\max }\Vert \mathcal{H}{f}^s\Vert +\Vert \mathcal{H}{f}^{n-1}\Vert \right)\\ \le C\left(\Vert \mathcal{H}{\mathcal{U}}^0\Vert +\underset{1\le s\le n-1}{\max }\Vert \mathcal{H}{f}^s\Vert \right).\end{array}$ (12)

使用引理3.6，可得

$\Vert U^n\Vert \le C\left(\Vert U^0\Vert +\underset{0\le s\le n-1}{\max }\Vert f^s\Vert \right),$

差分格式(6)是无条件稳定的，定理得证。 □

定理3.8说明，所构造离散格式的每一层值，都受初值的控制，因此证明了离散格式是无条件稳定的。

接着，进行离散格式(6)的误差估计。

定理3.9. 若 $u^n=u\left(x,t_n\right)$ 是方程(1)的精确解， $U^n$ 是离散格式(6)的数值解，并且 $e^n=u^n-U^n$ 。则有

$\Vert e^n\Vert \le C\left({\tau }^2+h^4\right)$ 。

其中，C是一个常数。

证明. 将(5)式和(6)式做差，可得

$\mathcal{H}{e}^{n+1}-\frac{a^2{\tau }^2}{12}{\delta }_x^2{e}^{n+1}=2\mathcal{H}{e}^n-\mathcal{H}{e}^{n-1}+\frac{10a^2{\tau }^2}{12}{\delta }_x^2{e}^n+\frac{a^2{\tau }^2}{12}{\delta }_x^2{e}^{n-1}+{\tau }^2{R}^n,$ (13)

其中 $\left|R^n\right|\le C\left({\tau }^2+h^4\right)$ 。对于上式左右两端同时乘以 $\mathcal{H}{e}^{n+1}$ 做内积，可得

$\begin{array}{l}\langle \mathcal{H}{e}^{n+1},\mathcal{H}{e}^{n+1}\rangle -\frac{a^2{\tau }^2}{12}\langle {\delta }_x^2{e}^{n+1},\mathcal{H}{e}^{n+1}\rangle \\ =2\langle \mathcal{H}{e}^n,\mathcal{H}{e}^{n+1}\rangle -\langle \mathcal{H}{e}^{n-1},\mathcal{H}{e}^{n+1}\rangle +\frac{10a^2{\tau }^2}{12}\langle {\delta }_x^2{e}^n,\mathcal{H}{e}^{n+1}\rangle +\frac{a^2{\tau }^2}{12}\langle {\delta }_x^2{e}^{n-1},\mathcal{H}{e}^{n+1}\rangle +{\tau }^2\langle R^n,\mathcal{H}{e}^{n+1}\rangle .\end{array}$

由于 $\langle {\delta }_x^2{u}^k,u^n\rangle =-\langle {\delta }_x{u}^k,{\delta }_x{u}^n\rangle$ 和引理3.5，我们可知

$\begin{array}{l}\langle \mathcal{H}{e}^{n+1},\mathcal{H}{e}^{n+1}\rangle \le 2\langle \mathcal{H}{e}^n,\mathcal{H}{e}^{n+1}\rangle -\langle \mathcal{H}{e}^{n-1},\mathcal{H}{e}^{n+1}\rangle -\frac{10a^2{\tau }^2}{12}{\left(e^n,e^{n+1}\right)}_{\mathcal{H}}\\ -\frac{a^2{\tau }^2}{12}{\left(e^{n-1},e^{n+1}\right)}_{\mathcal{H}}+{\tau }^2\langle R^n,\mathcal{H}{e}^{n+1}\rangle .\end{array}$ (14)

对于(14)式，根据范数等价定理，Cauchy-Schwarz不等式及引理3.6可得

$\begin{array}{l}{\Vert \mathcal{H}{e}^{n+1}\Vert }^2\le 2\Vert \mathcal{H}{e}^n\Vert \Vert \mathcal{H}{e}^{n+1}\Vert +\Vert \mathcal{H}{e}^{n-1}\Vert \Vert \mathcal{H}{e}^{n+1}\Vert +\frac{10a^2{\tau }^2}{12}{\Vert e^n\Vert }_{\mathcal{H}}{\Vert e^{n+1}\Vert }_{\mathcal{H}}\\ +\frac{a^2{\tau }^2}{12}{\Vert e^{n-1}\Vert }_{\mathcal{H}}{\Vert e^{n+1}\Vert }_{\mathcal{H}}+{\tau }^2\Vert R^n\Vert \Vert \mathcal{H}{e}^{n+1}\Vert \\ \le 2\Vert \mathcal{H}{e}^n\Vert \Vert \mathcal{H}{e}^{n+1}\Vert +\Vert \mathcal{H}{e}^{n-1}\Vert \Vert \mathcal{H}{e}^{n+1}\Vert +\frac{C_1{a}^2{\tau }^2}{12}\Vert e^n\Vert \Vert e^{n+1}\Vert \\ +\frac{C_2{a}^2{\tau }^2}{12}\Vert e^{n-1}\Vert \Vert e^{n+1}\Vert +{\tau }^2\Vert R^n\Vert \Vert \mathcal{H}{e}^{n+1}\Vert \\ \le 2\Vert \mathcal{H}{e}^n\Vert \Vert \mathcal{H}{e}^{n+1}\Vert +\Vert \mathcal{H}{e}^{n-1}\Vert \Vert \mathcal{H}{e}^{n+1}\Vert +\frac{3C_1{a}^2{\tau }^2}{12}\Vert \mathcal{H}{e}^n\Vert \Vert \mathcal{H}{e}^{n+1}\Vert \\ +\frac{3C_2{a}^2{\tau }^2}{12}\Vert \mathcal{H}{e}^{n-1}\Vert \Vert \mathcal{H}{e}^{n+1}\Vert +{\tau }^2\Vert R^n\Vert \Vert \mathcal{H}{e}^{n+1}\Vert .\end{array}$

上式化简得

$\begin{array}{l}\Vert \mathcal{H}{e}^{n+1}\Vert \le \left(2+\frac{3C_1{a}^2{\tau }^2}{12}\right)\Vert \mathcal{H}{e}^n\Vert +\left(1+\frac{3C_2{a}^2{\tau }^2}{12}\right)\Vert \mathcal{H}{e}^{n-1}\Vert +{\tau }^2\Vert R^n\Vert \\ \le C\left(\Vert \mathcal{H}{e}^n\Vert +\Vert \mathcal{H}{e}^{n-1}\Vert +\Vert R^n\Vert \right).\end{array}$ (15)

对于(15)式，当 $n=0$ 时，有

$\Vert \mathcal{H}{e}^1\Vert \le C\Vert \mathcal{H}{e}^0\Vert +C\Vert \mathcal{H}{e}^{-1}\Vert +\Vert R^0\Vert ,$ (16)

为了证明结论，首先需要对 $\Vert \mathcal{H}{e}^{-1}\Vert$ 进行估计。由引理3.7可得

$\begin{array}{c}\Vert e^{-1}\Vert =\Vert u\left(x,t_{-1}\right)-U^{-1}\Vert \\ =\Vert u\left(x,t_{-1}\right)-u^{-1}+u^{-1}-U^{-1}\Vert \\ =\Vert {\epsilon }^{-1}+u^0-\tau \psi -\left(U^0-\tau \psi \right)\Vert \\ \le C\left({\tau }^2\right).\end{array}$

故， $\Vert \mathcal{H}{e}^{-1}\Vert \le \Vert e^{-1}\Vert \le C\left({\tau }^2\right)$ 。

结合(16)式可知 $\Vert \mathcal{H}{e}^1\Vert \le C\left({\tau }^2+h^4\right)$ 。

现在我们假设当 $k=1,2,\cdots ,n-1$ 时，(15)描述的结论也成立，即

$\Vert \mathcal{H}{e}^k\Vert \le C\left({\tau }^2+h^4\right).$

接下来我们将考虑 $k=n$ ，根据(15)式及假设，可得

$\begin{array}{l}\Vert \mathcal{H}{e}^n\Vert \le C\left(\Vert \mathcal{H}{e}^{n-1}\Vert +\Vert \mathcal{H}{e}^{n-2}\Vert +\Vert R^{n-1}\Vert \right)\\ \le C\left({\tau }^2+h^4\right).\end{array}$

以引理3.6为依据，可得

$\Vert e^n\Vert \le \sqrt{3}\Vert \mathcal{H}{e}^n\Vert \le \sqrt{3}C\left({\tau }^2+h^4\right),$

因此

$\Vert e^n\Vert \le C\left({\tau }^2+h^4\right),$

以上定理证明完毕。 □

定理3.9说明了方程的数值解和精确解之间的误差随着时间和空间剖分的增大而减小，从侧面说明了数值解和精确解在理论分析上是很接近的。

4. 数值实验

在本节中，我们将进行数值实验去证明对紧致有限差分格式理论分析的正确性。数值实验是在Intel酷睿i5-8265U，内存4 GB和主频1.60 GHz的计算机环境下，使用MATLAB2017a进行数值实验。用 $L^{\infty }$ 范数来衡量数值误差。数值解的 $L^{\infty }$ 范数误差、时间和空间方向的收敛阶分别用下面的公式计算：

$L^{\infty }\text{-}误差=\underset{0\le n\le N}{\max }{\Vert u^n-U^n\Vert }_{\infty },$

$时间收敛阶=\log \left(e_1/e_2\right)/\log \left({\tau }_1/{\tau }_2\right),$

其中， $e_1$ 和 $e_2$ 分别为时间步长 ${\tau }_1$ 和 ${\tau }_2$ 所对应的误差。

$空间收敛阶=\log \left(e_1/e_2\right)/\log \left(h_1/h_2\right),$

其中， $e_1$ 和 $e_2$ 分别为空间步长 $h_1$ 和 $h_2$ 所对应的误差。

例1. 考虑下面的时间波动方程

$\{\begin{array}{l}\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}=a^2\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+f\left(x,t\right),x\in \left(0,1\right),t\in \left(0,1\right],\\ u\left(x,0\right)=\sin \left(2\pi x\right),u_t(x,0)=0,x\in \left[0,1\right],\\ u\left(0,t\right)=0,u\left(1,t\right)=0,t\in \left[0,1\right].\end{array}$

其中， $a=1$ ，精确解 $u\left(x,t\right)=\left(t^3+1\right)\sin \left(2\pi x\right)$ ，并且

$f\left(x,t\right)=\left(6t+4{\pi }^2\left(t^3+1\right)\right)\sin \left(2\pi x\right).$

例1的实验结果将会呈现在表1、表2和图1中。表1展示了紧致有限差分格式(6)的最大范数误差、时间收敛阶和CPU时间。我们选择合适的空间步长 $h=1/100$ ，时间步长取不同的剖分。我们很容易可以看到时间方向的收敛阶为 $O\left({\tau }^2\right)$ 。表2展示了紧致有限差分格式的最大范数误差、空间收敛阶和CPU时间。此时，空间步长和时间步长都取不同的剖分。我们很容易可以看到空间方向的收敛阶为 $O\left(h^4\right)$ ，这就验证了定理3.9。

Table 1. L ∞ -error, time convergence order and CPU time when M = 100

表1. 当M = 100时的 $L^{\infty }$ -误差，时间收敛阶和CPU时间

Table 2. L ∞ -error, spatial convergence order and CPU time when M and N are taken into different subdivisions

表2. 当M，N取不同剖分时的 $L^{\infty }$ -误差，空间收敛阶和CPU时间

图1展示了当 $N=100,M=5000$ 时的精确解、数值解、绝对误差和误差的等高线图。

(a) 精确解 (b) 数值解 (c) 绝对误差 (d) 绝对误差等高线

Figure 1. Contour plot of exact solution, numerical solution, absolute error and error when $N=100,M=5000$

图1. 当 $N=100,M=5000$ 时的精确解、数值解、绝对误差和误差的等高线图

通过上面四组图，我们很容易看出，离散格式的数值解是很接近原方程的精确解的。

5. 总结

本文针对满足一定的初边值条件的整数阶波动方程，结合紧致有限差分算子，构造了一种高精度的差分格式，并且用能量不等式的方法证明了差分格式的稳定性以及收敛性，通过理论分析得到差分格式的收敛精度在时间上达到二阶，在空间上达到四阶。最后通过数值实验证明了理论分析的正确性，即所构造差分格式的收敛阶为 $O\left({\tau }^2+h^4\right)$ 。

本文研究的模型虽然达到了高精度的收敛阶，但是研究的模型不具备足够的新颖性和可应用性。在接下来的研究中，作者将会致力于研究分数阶波动方程的紧致有限差分方法，在时间维度上将会考虑一维和二维。

基金项目

新疆应用职业技术学院科研项目(XYZY2023KQL016)。

NOTES

^*通讯作者。