非线性延迟波动方程的紧致有限差分格式

A Compact Finite Difference Scheme for Nonlinear Wave Equation with Delay

摘要: 本文主要研究了一类带有Dirichlet边界条件的非线性延迟波动方程，并建立了一个紧致有限差分格式。运用离散能量法证明了该差分格式在L_∞范数下满足时间二阶、空间四阶的收敛率。最后通过数值算例验证了算法的精度和有效性。

Abstract: A compact finite difference scheme is established for a class of nonlinear wave equations with delay with Dirichlet boundary value conditions. By using the discrete energy method, the proposed compact finite difference scheme is proved temporally second-order convergence rate and spatially fourth-order convergence rate in L_∞ norm. Finally, numerical results have confirmed the accuracy and effectiveness of the algorithm.

1. 引言

在描述自然科学和社会科学中的各种现象时，常常需要利用系统过去时刻的状态，描述即将发生的状态，延迟微分方程模型由此诞生。延迟微分方程广泛应用于经济学、物理学、生态学、医学、交通调度、工程控制、计算机辅助设计、核工程等许多科学领域，因此研究此类问题具有非常重要的意义和实用价值。随着对延迟微分方程的研究，其中的一个分支，延迟偏微分方程逐渐受到大家的重视，而其数值方法的研究可以极大弥补理论研究方面的不足。

近些年来，在延迟抛物方程的数值研究方面，已有不少成果。孙志忠等 [1] 研究了非线性延迟抛物方程的Crank-Nicolson差分格式，池永日 [2] 和范乐乐等 [3] 分别研究了非线性延迟抛物方程的紧差分格式。然而，人们对延迟波动方程的研究关注不多，陈景良等 [4] 研究了非线性延迟波动方程的显式有限差分格式，该格式具有时间二阶、空间二阶的精度。为了更有效地解决问题，人们对于计算结果的精度要求越来越高，传统的显式差分格式必须在增加网格点数量的前提下，才能提高计算精度，但是这也增加了计算时间，而紧致有限差分格式仅仅利用三个网格点，就可以达到四阶的收敛精度，所以具有很高的研究意义，因此本文研究非线性延迟波动方程的高阶紧致有限差分方法。

文章安排如下，在第二部分给出了非线性延迟波动方程的紧致有限差分格式，并对其收敛性进行理论性分析。在第三部分，通过数值算例验证了理论结果。最后一部分，给出了文章的结论。

2. 紧致有限差分格式

考虑如下非线性延迟波动方程：

$\frac{{\partial }^{2}u}{\partial t^2}=a^2\frac{{\partial }^{2}u}{\partial x^2}+f\left(x,t,u\left(x,t-s\right)\right),b_1\le x\le b_2,0<t\le T,$ (1a)

其初边值条件为：

$u\left(x,t\right)=\varphi \left(x,t\right),\frac{\partial u}{\partial t}\left(x,0\right)=\phi \left(x\right),b_1\le x\le b_2,-s\le t\le 0,$ (1b)

$u\left(b_1,t\right)=\alpha \left(t\right),u\left(b_2,t\right)=\beta \left(t\right),0<t\le T.$ (1c)

其中 $a>0$ 为扩散系数，s为延迟系数。

2.1. 记号

为建立非线性延迟波动方程(1a)~(1c)的差分格式，首先对 $\Omega =\left\{\left(x,t\right)|b_1\le x\le b_2,-s\le t\le T\right\}$ 做均匀网格剖分。取空间步长 $h_x=\left(b_2-b_1\right)/m$，取时间步长 $h_t=s/n_1=T/n$。这里 $m,n_1,n$ 为整数。记 $x_i=b_1+i{h}_x$， $t_k=k{h}_t$， $0\le i\le m$， $-n_1\le k\le n$，这里 $i,k$ 均为整数。分别记结点 $\left(x_i,t_k\right)$ 处精确解和数值解为 $u_i^k,U_i^k$。记网格剖分区域

${\Omega }_h=\left\{\left(x_i,t_k\right)|0\le i\le m,-n_1\le k\le n\right\},$

定义网格函数空间

$U_h=\left\{U|U=U_i,0\le i\le m,U_0=U_m=0\right\},$

对任意 $U_i^k\in U_h$，定义算子、内积和范数，

${\delta }_x{U}_{i+\frac{1}{2}}^k=\frac{1}{2h_x}\left(U_{i+1}^k-U_i^k\right),{\delta }_t{U}_i^{k+\frac{1}{2}}=\frac{1}{h_t}\left(U_i^{k+1}-U_i^k\right),{\delta }_t^2{U}_i^k=\frac{1}{h_t^2}\left(U_i^{k+1}-2U_i^k+U_i^{k-1}\right),$

${\delta }_x^2{U}_i^k=\frac{1}{h_x^2}\left(U_{i+1}^k-2U_i^k+U_{i-1}^k\right),\left(U,V\right)=h_x\underset{i=1}{\overset{m-1}{{\displaystyle \sum }}}{U}_i{V}_i,{\left(U,V\right)}_1=h_x\underset{i=1}{\overset{m-1}{{\displaystyle \sum }}}\left({\delta }_x{U}_{i-\frac{1}{2}}\right)\left({\delta }_x{V}_{i-\frac{1}{2}}\right),$

$\Vert U\Vert =\sqrt{\left(U,U\right)},{\left|U\right|}_1=\sqrt{{\left(U,U\right)}_1}.$

定义紧致差分算子，

${\left(\mathcal{A}U\right)}_i=\{\begin{array}{l}\frac{1}{12}\left(U_{i+1}+10U_i+U_{i-1}\right),1\le i\le m-1,\\ U_i,i=0,m.\end{array}$

下面构造紧致有限差分格式。

2.2. 差分格式的建立

在点 $\left(x_i,t_k\right)$ 处考虑方程(1a)有

$\frac{{\partial }^{2}u}{\partial t^2}\left(x_i,t_k\right)-a^2\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_k\right)=f\left(x_i,t_k,u_i^{k-n_1}\right),$ (2)

由泰勒公式可知，

$\frac{{\partial }^{2}u}{\partial t^2}\left(x_i,t_k\right)={\delta }_t^2{u}_i^k-\frac{h_t^2}{12}\frac{{\partial }^{4}u}{\partial t^4}\left(x_i,{\theta }_{ik}\right),{\theta }_{ik}\in \left(t_{k-1},t_{k+1}\right),$

$\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_k\right)=\frac{1}{2}\left[\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_{k+1}\right)+\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_{k-1}\right)\right]-\frac{h_t^2}{2}\frac{{\partial }^{4}u}{\partial x^2\partial t^2}\left(x_i,\widetilde{{\theta }_{ik}}\right),\widetilde{{\theta }_{ik}}\in \left(t_{k-1},t_{k+1}\right).$

将上式代入方程(2)得，

$\begin{array}{l}{\delta }_t^2{u}_i^k-\frac{1}{2}{a}^2\left[\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_{k+1}\right)+\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_{k-1}\right)\right]\\ =f\left(x_i,t_k,u_i^{k-n_1}\right)+\left(\frac{1}{12}\frac{{\partial }^{4}u}{\partial t^4}\left(x_i,{\theta }_{ik}\right)-\frac{a^2}{2}\frac{{\partial }^{4}u}{\partial x^2\partial t^2}\left(x_i,\widetilde{{\theta }_{ik}}\right)\right)h_t^2.\end{array}$

将紧算子 $\mathcal{A}$ 作用于上式两端得，

$\begin{array}{l}\mathcal{A}{\delta }_t^2{u}_i^k-\frac{1}{2}{a}^2\left[\mathcal{A}\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_{k+1}\right)+\mathcal{A}\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_{k-1}\right)\right]\\ =\mathcal{A}f\left(x_i,t_k,u_i^{k-n_1}\right)+\mathcal{A}\left(\frac{1}{12}\frac{{\partial }^{4}u}{\partial t^4}\left(x_i,{\theta }_{ik}\right)-\frac{a^2}{2}\frac{{\partial }^{4}u}{\partial x^2\partial t^2}\left(x_i,\widetilde{{\theta }_{ik}}\right)\right)h_t^2.\end{array}$ (3)

根据泰勒公式，可得，

$\mathcal{A}\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_k\right)={\delta }_x^2{u}_i^k+\frac{h_x^4}{240}\frac{{\partial }^{6}u}{\partial x^6}\left({\epsilon }_{ik},t_k\right),{\epsilon }_{ik}\in \left(x_{i-1},x_{i+1}\right),$

则

$\frac{1}{2}\mathcal{A}\left[\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_{k+1}\right)+\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_{k-1}\right)\right]={\delta }_x^2{u}_i^k+\frac{h_x^4}{240}\frac{{\partial }^{6}u}{\partial x^6}\left({\epsilon }_{ik},t_{ik}\right),{\epsilon }_{ik}\in \left(x_{i-1},x_{i+1}\right),t_{ik}\in \left(t_{k-1},t_{k+1}\right).$ (4)

将方程(4)代入方程(3)中，得

$\mathcal{A}{\delta }_t^2{u}_i^k-a^2{\delta }_x^2{u}_i^k=\mathcal{A}f\left(x_i,t_k,u_i^{k-n_1}\right)+R_{ik},1\le i\le m-1,0\le k\le n,$ (5)

其中

$R_{ik}=\mathcal{A}\left[\frac{1}{12}\frac{{\partial }^{4}u}{\partial t^4}\left(x_i,{\theta }_{ik}\right)-\frac{a^2}{2}\frac{{\partial }^{4}u}{\partial t^4}\left(x_i,\widetilde{{\theta }_{ik}}\right)\right]h_t^2+\frac{a^2}{240}\frac{{\partial }^{6}u}{\partial x^6}\left(\widetilde{{\epsilon }_{ik}},\widetilde{t_{ik}}\right)h_x^4.$

假设解在充分光滑的情况下，则 $R_{ik}=C\left(h_t^2+h_x^4\right)$，本文中C为不同的任意常数。

用数值解 $U_i^k$ 代替精确解 $u_i^k$，并省略 $R_{ik}$，得紧致差分格式：

$\mathcal{A}{\delta }_t^2{U}_i^k-a^2{\delta }_x^2{U}_i^k=\mathcal{A}f\left(x_i,t_k,U_i^{k-n_1}\right),1\le i\le m-1,0\le k\le n.$ (6)

由于紧致差分格式(6)为三层格式，而我们仅知道初始条件中 $U_i^0$ 的值，下面讨论如何求解 $U_i^1$。

由方程(1b)可知，

$u_i^0=\phi \left(x_i\right),1\le i\le m-1,$

在点 $\left(x_i,t_0\right)$ 处考虑方程(1a)可知，

$\frac{{\partial }^{2}u}{\partial t^2}\left(x_i,t_0\right)=a^2\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_0\right)+f\left(x_i,t_0,u_i^{-n_1}\right)=a^2{\varphi }^{″}\left(x\right)+f\left(x_i,t_0,u_i^{-n_1}\right),$

由泰勒公式及方程(1b)可知，

$\begin{array}{l}{u}_i^1=u\left(x_i,t_0\right)+h_t\frac{\partial u}{\partial t}\left(x_i,t_0\right)+\frac{1}{2}{h}_t^2\frac{{\partial }^{2}u}{\partial t^2}\left(x_i,t_0\right)+\frac{1}{2}{\displaystyle {\int }_0^{h_t}{\left(h_t-t\right)}^2\frac{{\partial }^{3}u}{\partial t^3}\left(x_i,t_0\right)\text{d}t}\\ =\varphi \left(x_i\right)+h_t\phi \left(x_i\right)+\frac{1}{2}{h}_t{}^2\left[a^2{\varphi }^{″}\left(x_i\right)+f\left(x_i,t_0,u_i^{-n_1}\right)\right]+R_{i0},1\le i\le m-1.\end{array}$

其中

$R_{i0}=\mathcal{A}\frac{1}{2}{\displaystyle {\int }_0^{h_t}{\left(h_t-t\right)}^2\frac{{\partial }^{3}u}{\partial t^3}\left(x_i,t_0\right)\text{d}t}=C\left(h_t^2+h_x^4\right).$

再次用数值解 $U_i^k$ 代替精确解 $u_i^k$，并省略 $R_{i0}$，得

$U_i^1=\varphi \left(x_i\right)+h_t\phi \left(x_i\right)+\frac{1}{2}{h}_t{}^2\left[a^2{\varphi }^{″}\left(x_i\right)+f\left(x_i,t_0,U_i^{-n_1}\right)\right],1\le i\le m-1.$ (7)

根据上面的推导过程，得到如下紧致有限差分格式，

$\mathcal{A}{\delta }_t^2{U}_i^k-a^2{\delta }_x^2{U}_i^k=\mathcal{A}f\left(U_i^{k-n_1},x_i,t_k\right),1\le i\le m-1,0\le k\le n,$ (8a)

$U_i^1=\varphi \left(x_i\right)+h_t\phi \left(x_i\right)+\frac{1}{2}{h}_t^2\left[a^2{\varphi }^{″}\left(x_i\right)+f\left(U_i^{-n_1},x_i,t_0\right)\right],1\le i\le m-1,$ (8b)

$U_i^k=\varphi \left(x_i,t_k\right),1\le i\le m-1,$ (8c)

$U_0^k=\alpha \left(t_k\right),U_m^k=\beta \left(t_k\right),1\le t\le n.$ (8d)

注：紧致差分算子 $\mathcal{A}$ 作用于 $U^k$，可写为矩阵形式，即

$\mathcal{A}{U}^k=\left[\begin{array}{ccccccc}10& 1& 0& & 0& 0& 0\\ 1& 10& 1& \cdots & 0& 0& 0\\ 0& 1& 10& & 0& 0& 0\\ & ⋮& & \ddots & & ⋮& \\ 0& 0& 0& & 10& 1& 0\\ 0& 0& 0& \cdots & 1& 10& 1\\ 0& 0& 0& & 0& 1& 10\end{array}\right]\left[\begin{array}{c}{U}_i^k\\ U_i^k\\ U_i^k\\ ⋮\\ U_i^k\\ U_i^k\\ U_i^k\end{array}\right],$

可见，紧致差分算子 $\mathcal{A}$ 是对角占优的，因此 $\mathcal{A}$ 是可逆算子，即满足 $\frac{{\partial }^{2}u}{\partial x^2}\left(x_i,t_k\right)={\mathcal{A}}^{-1}{\delta }_x^2{u}_i^k+O\left(h^4\right)$。

2.3. 差分格式的收敛性分析

下面介绍三个重要引理来研究紧致差分格式(8a)~(8d)的收敛性。

引理1 [5] 设 $v\in U_h$，则有下列不等式成立

$\left(-{\delta }_x^{2}v,v\right)={\Vert {\delta }_{x}v\Vert }^2,{\Vert v\Vert }_{\infty }\le \frac{\sqrt{b_2-b_1}}{2}{\left|v\right|}_1,$

$\Vert v\Vert \le \frac{b_2-b_1}{\sqrt{6}}{\left|v\right|}_1,h_x{\Vert {\delta }_{x}v\Vert }^2\le 4{\Vert v\Vert }^2.$

引理2 [6] 设 $\left\{F^k|k\ge 0\right\}$ 是非负序列且满足 $F^{k+1}\le A+B{h}_t{\displaystyle {\sum }_{k=1}^K{F}^k}$，则

$\underset{0\le k\le K+1}{\max }{F}^k\le A\exp \left(2B\left(K+1\right)h_t\right),$

其中 $A,B$ 为非负常数。

引理3 [7] 对于对称正定矩阵 ${\mathcal{A}}^{-1}$，我们有

${\mathcal{A}}^{-1}=R_1^{\text{T}}{R}_1$ 和 ${\mathcal{C}}_0\Vert U^n\Vert \le \Vert R_1{U}^n\Vert \le {\mathcal{C}}_1\Vert U^n\Vert ,$

其中 $R_1=\mathcal{C}hol\left({\mathcal{A}}^{-1}\right)$， ${\mathcal{C}}_i\left(i=0,1\right)$ 为正常数。

定理1假设问题(1a)~(1c)的精确解u足够光滑，令紧致差分格式(8a)~(8d)的数值解为 $U_i^k$。记 $e_i^k=u_i^k-U_i^k$，当 $\frac{a^2{h}_t^2}{h_x^2}<1$ 且 $f\left(u,x,t\right)$ 满足局部Lipschitz条件时，则

${\left|e^k\right|}_1^2\le C{\left(h_t^2+h_x^4\right)}^2,-n_1\le k\le n,$ (9)

${\Vert e^k\Vert }_{\infty }\le C\left(h_t^2+h_x^4\right),-n_1\le k\le n,$ (10)

成立。

证明：将方程(5)与(6)相减，得到

$\mathcal{A}{\delta }_t^2{e}_i^k-a^2{\delta }_x^2{e}_i^k=\mathcal{A}f\left(x_i,t_k,u_i^{k-n_1}\right)-\mathcal{A}f\left(x_i,t_k,U_i^{k-n_1}\right)+R_{ik},1\le i\le m-1,0\le k\le n,$ (11)

方程(11)的两端同时乘以 ${\mathcal{A}}^{-1}$，运用引理3，可得误差方程

${\delta }_t^2{e}_i^k-R_1{R}_1^{\text{T}}{a}^2{\delta }_x^2{e}_i^k=f\left(x_i,t_k,u_i^{k-n_1}\right)-f\left(x_i,t_k,U_i^{k-n_1}\right)+R_1^{\text{T}}{R}_1{R}_i^k,1\le i\le m,1\le m<n,$ (11a)

$e_i^k=0,0\le i\le m,-n_1\le k\le 0,$ (11b)

$e_i^k=0,i=0,m,0\le k\le n,$ (11c)

当 $-n_1\le k\le 0$ 时，由方程(11b)可知， $e_i^k=0$。显然误差满足方程(9)~(10).

下面我们利用数学归纳法证明误差。假设当 $k=l$ 时，方程(9)~(10)均成立，应用引理1，则

${\Vert e^k\Vert }_{\infty }\le \frac{\sqrt{b_2-b_1}}{2}{\left|e^k\right|}_1\le \epsilon ,1\le k\le l.$

方程(11a)的两端同时与 $2{\delta }_{\hat{t}}{e}_k$ 做内积，则

$\left({\delta }_t^2{e}_i^k,2{\delta }_{\hat{t}}{e}_k\right)-\left(R_1{R}_1^{\text{T}}{a}^2{\delta }_x^2{e}_i^k,2{\delta }_{\hat{t}}{e}_k\right)=\left(f\left(x_i,t_k,u_i^{k-n_1}\right)-f\left(x_i,t_k,U_i^{k-n_1}\right),2{\delta }_{\hat{t}}{e}_k\right)+\left(R_1^{\text{T}}{R}_1{R}_i^k,2{\delta }_{\hat{t}}{e}_k\right).$

运用引理1，得

$\begin{array}{c}\left({\delta }_t^2{e}_i^k,2{\delta }_{\hat{t}}{e}_k\right)=\left(\frac{e^{k+1}-2e^k+e^{k-1}}{h_t^2},\frac{e^{k+1}-e^{k-1}}{h_t}\right)\\ =\frac{1}{h_t^3}\left(e^{k+1}-e^k-\left(e^k-e^{k-1}\right),e^{k+1}-e^k+\left(e^k-e^{k-1}\right)\right)\\ =\frac{1}{h_t}\left({\delta }_t{\Vert e^{k+\frac{1}{2}}\Vert }^2-{\delta }_t{\Vert e^{k-\frac{1}{2}}\Vert }^2\right)\end{array}$ (12)

和

$\begin{array}{c}\left(R_1{R}_1^{\text{T}}{a}^2{\delta }_x^2{e}_i^k,2{\delta }_{\hat{t}}{e}_k\right)=-\left(R_1{\delta }_x{e}^k,2R_1{\delta }_x{\delta }_{\hat{t}}{e}_k\right)=-\left(R_1{\delta }_x{e}^k,\frac{R_1{\delta }_x{e}^{k+1}-R_1{\delta }_x{e}^{k-1}}{h_t}\right)\\ =-\frac{1}{h_t}\left[\left(R_1{\delta }_x{e}^k,R_1{\delta }_x{e}^{k+1}\right)-\left(R_1{\delta }_x{e}^{k-1},R_1{\delta }_x{e}^k\right)\right].\end{array}$ (13)

另外，

$\left(f\left(x_i,t_k,u_i^{k-n_1}\right)-f\left(x_i,t_k,U_i^{k-n_1}\right),2{\delta }_{\hat{t}}{e}_k\right)\le \frac{1}{2}{\Vert f\left(x_i,t_k,u_i^{k-n_1}\right)-f\left(x_i,t_k,U_i^{k-n_1}\right)\Vert }^2+\frac{1}{2}{\Vert 2{\delta }_{\hat{t}}{e}_k\Vert }^2.$

由于 $f\left(u,x,t\right)$ 满足局部Lipschitz条件，则

${\Vert f\left(x_i,t_k,u_i^{k-n_1}\right)-f\left(x_i,t_k,U_i^{k-n_1}\right)\Vert }^2\le C{\Vert e^{k-n_1}\Vert }^2.$ (14)

由柯西不等式，计算得

$\frac{1}{2}{\Vert 2{\delta }_{\hat{t}}{e}_k\Vert }^2\le {\Vert {\delta }_t{e}^{k+\frac{1}{2}}\Vert }^2+{\Vert {\delta }_t{e}^{k-\frac{1}{2}}\Vert }^2,$

$\left(R_1^{\text{T}}{R}_1{R}_i^k,2{\delta }_{\hat{t}}{e}_k\right)\le C{\left(h_t^2+h_x^4\right)}^2+{\Vert {\delta }_t{e}^{k+\frac{1}{2}}\Vert }^2+{\Vert {\delta }_t{e}^{k-\frac{1}{2}}\Vert }^2.$ (15)

记 $L^k={\Vert {\delta }_t{e}^{k+\frac{1}{2}}\Vert }^2+a^2\left(R_1{\delta }_x{e}^k,R_1{\delta }_x{e}^{k+1}\right)$，结合方程(12)~(15)，可得

$\frac{1}{h_t}\left(L^k-L^{k-1}\right)\le C{\Vert e^{k-n_1}\Vert }^2+2{\Vert {\delta }_t{e}^{k+\frac{1}{2}}\Vert }^2+2{\Vert {\delta }_t{e}^{k-\frac{1}{2}}\Vert }^2+C{\left(h_t^2+h_x^4\right)}^2.$ (16)

下面我们分析 $L^k$。当 $q=\frac{a^2{h}_t^2}{h_x^2}<1$ 时，有

$\begin{array}{c}{L}^k={\Vert {\delta }_t{e}^{k+\frac{1}{2}}\Vert }^2+a^2{\left|R_1{e}^{k+\frac{1}{2}}\right|}_1^2-\left[a^2{\left|R_1{e}^{k+\frac{1}{2}}\right|}_1^2-a^2\left(R_1{\delta }_x{e}^k,R_1{\delta }_x{e}^{k+1}\right)\right]\\ =\left(1-q\right){\Vert {\delta }_t{e}^{k+\frac{1}{2}}\Vert }^2+a^2{\left|R_1{e}^{k+\frac{1}{2}}\right|}_1^2+q{\Vert {\delta }_t{e}^{k+\frac{1}{2}}\Vert }^2-\left[a^2{\left|R_1{e}^{k+\frac{1}{2}}\right|}_1^2-a^2\left(R_1{\delta }_x{e}^k,R_1{\delta }_x{e}^{k+1}\right)\right]\\ \ge \left(1-q\right){\Vert {\delta }_t{e}^{k+\frac{1}{2}}\Vert }^2+a^2{\left|R_1{e}^{k+\frac{1}{2}}\right|}_1^2.\end{array}$ (17)

由上式易知

${\Vert {\delta }_t{e}^{k+\frac{1}{2}}\Vert }^2\le \frac{1}{1-q}{L}^k,{\left|e^{k+\frac{1}{2}}\right|}_1^2\le \frac{1}{a^2}{L}^k.$ (18)

因为 $e_i^{k+1}=e_i^{k+\frac{1}{2}}+h_t{\delta }_t{e}_i^{k+\frac{1}{2}}/2$，则

$\begin{array}{c}{\left({\delta }_x{e}_{i+\frac{1}{2}}^{k+1}\right)}^2={\left({\delta }_x{e}_{i+\frac{1}{2}}^{k+\frac{1}{2}}+\frac{h_t}{2}{\delta }_x{\delta }_t{e}_{i+\frac{1}{2}}^{k+1}\right)}^2={\left[{\delta }_x{e}_{i+\frac{1}{2}}^{k+\frac{1}{2}}+\frac{h_t}{2h_x}\left({\delta }_t{e}_{i+1}^{k+\frac{1}{2}}-{\delta }_t{e}_i^{k+\frac{1}{2}}\right)\right]}^2\\ \le 2{\left({\delta }_x{e}_{i+\frac{1}{2}}^{k+\frac{1}{2}}\right)}^2+\frac{q}{a^2}\left[{\left({\delta }_t{e}_{i+1}^{k+\frac{1}{2}}\right)}^2+{\left({\delta }_t{e}_i^{k+\frac{1}{2}}\right)}^2\right],0\le i\le m-1,\end{array}$ (19)

在方程(19)两端同时乘以 $a^2{h}_x$ 并对i从1到 $m-1$ 求和，应用方程(18)可得

$a^2{\left|e^{k+1}\right|}_1^2\le 2L^k+\frac{2}{1-q}{L}^k\le \frac{2}{1-q}{L}^k,1\le k\le l.$ (20)

将方程(20)直接代入方程(16)，得

$\frac{1}{h_t}\left(L^k-L^{k-1}\right)\le \frac{2}{1-q}\left(L^k+L^{k-1}\right)+C{\Vert e^{k-n_1}\Vert }^2+C{\left(h_t^2+h_x^4\right)}^2,$

运用引理1，

$\begin{array}{c}\frac{1}{h_t}\left(L^k-L^{k-1}\right)\le \frac{2}{1-q}\left(L^k+L^{k-1}\right)+C{\left|e^{k-n_1}\right|}_1^2+C{\left(h_t^2+h_x^4\right)}^2\\ \le \frac{2}{1-q}\left(L^k+L^{k-1}\right)+C{L}^{k-n_1-1}+C{\left(h_t^2+h_x^4\right)}^2.\end{array}$

在上式两端同时乘以 $h_t$，将上式中k用p替换，并对p从−1到k求和得

$L^k\le \frac{2h_t}{1-q}\underset{p=0}{\overset{k}{{\displaystyle \sum }}}\left(L^p+L^{p-1}\right)+C{h}_t\underset{p=0}{\overset{k-1}{{\displaystyle \sum }}}{L}^p+CT{\left(h_t^2+h_x^4\right)}^2,$

整理可得，

$L^k\le C{h}_t\underset{p=0}{\overset{k-1}{{\displaystyle \sum }}}{L}^p+CT{\left(h_t^2+h_x^4\right)}^2,$

运用引理2，得

$L^k\le C\left(h_t^2+h_x^4\right),0\le k\le l$ (21)

在方程(21)中，运用方程(20)和引理1，得

${\left|e^{l+1}\right|}_1^2\le \frac{2}{a^2\left(1-q^2\right)}{L}^l\le C{\left(h_t^2+h_x^4\right)}^2,$

${\Vert e^{l+1}\Vert }_{\infty }\le \frac{\sqrt{b_2-b_1}}{2}{\left|e^{l+1}\right|}_1\le C{\left(h_t^2+h_x^4\right)}^2.$

所以当 $k=l+1$ 时，方程(9)~(10)成立，由数学归纳法得证。

3. 数值实验

应用紧致差分格式(8a)~(8d)计算如下非线性延迟波动方程的初边值问题

$u_{tt}-2u_{xx}=f\left(x,t,u\left(x,t-0.1\right)\right),\left(x,t\right)\in \left[1,2\right]\times \left(0,1\right],$

$u\left(x,t\right)=\sin \left(xt\right),\left(x,t\right)\in \left[1,2\right]\times \left[-0.1,0\right),$

$u\left(1,t\right)=\sin t,u\left(2,t\right)=\sin \left(2t\right),t\in \left[0,1\right].$

其中 $f\left(x,t,u\left(x,t-0.1\right)\right)=-u^2\left(x,t-0.1\right)+{\sin }^2\left(x\left(t-0.1\right)\right)+\left(2t^2-x^2\right)\sin xt$。该问题的精确解为 $u\left(x,t\right)=\sin \left(xt\right)$。

令

$E_{\infty }\left(h_x,h_t\right)=\max \left|u_i^k-U_i^k\right|,$

$order={\log }_2{E}_{\infty }\left(2h_x,2h_t\right)/E_{\infty }\left(h_x,h_t\right).$

Table 1. When h t = h x 2 , space step error and convergence order

表1. 当 $h_t=h_x^2$ 时，空间步长的误差和收敛阶

Figure 1. Exact solution surface and numerical solution surface

图1. 精确解曲面和数值解曲面图

表1为紧致差分格式(8a)~(8d)在 $h_t=h_x^2$ 时取不同步长计算数值解得到的最大模误差，CPU为程序运行时间。从表1可以看到，我们所提的紧致差分格式(8a)~(8d)满足空间四阶的收敛率，这与理论结果一致。

图1给出了空间步长 $h_x=\frac{1}{40}$，时间步长 $h_t=\frac{1}{1600}$ 时的精确解曲面和数值解曲面。可以看到，精确解和数值解的三维图像是完全一致的。

4. 结论

非线性延迟波动方程属于时滞微分方程，所求的未知函数 $u\left(x,t\right)$ 在某确定时刻的导数由前面时刻函数所确定，所以它也是一类特殊的泛函偏微分方程。本文根据文献 [4] 讨论的传统显式差分格式，提出了一种高阶紧致有限差分格式，其中时间方向利用二阶中心差分方法、空间方向利用四阶紧算子进行离散，该格式的误差阶为 $O\left(h_t^2+h_x^4\right)$。最后通过数值算例，说明该方法具有一定的可行性和实用性。

基金项目

大学生创新创业计划项目(项目编号IEYCAUC2021020)。

NOTES

^*通讯作者。