Chebyshev谱法求解电报方程

Chebyshev Spectral Method for Solving Telegraph Equation

摘要: 双曲偏微分方程是重要的偏微分方程之一。提出求解电报方程的Chebyshev谱法，采用Chebyshev-Gauss-Lobatto配点，利用Chebyshev多项式构造导数矩阵，将电报方程近似为常微分方程，证明了电报方程的离散Chebyshev谱法的误差估计，采用Runge-Kutta进行求解。将该法得到的数值结果与精确解进行比较，验证了方法的有效性，数据结果的误差与其他方法相比有较高的精确度。

Abstract: Hyperbolic partial differential equation is one of the important partial differential equations. The Chebyshev spectral method is proposed to solve the telegraph equation. Chebyshev-gauss-lobatto is used to assign points, the derivative matrix is constructed by Chebyshev polynomial, and the telegraph equation is approximated as an ordinary differential equation. The error estimation of the discrete Chebyshev spectral method for the telegraph equation was proved. Runge-Kutta was used to solve the problem. The numerical results obtained by the method are compared with the exact solution, and the effectiveness of the method is verified. The error of the data results is more accurate than that of other methods.

1. 引言

谱方法[1] [2]是一种既古老又新兴的求解偏微分方程的数值方法。早在1820年，Navier就运用双重三角级数求解弹性薄板问题。谱方法有很多优点，最重要的是他的“无穷阶”收敛速度。Chebyshev谱法主要用于求解非周期性问题和具有光滑系数和简单域的数值微分方程[3]-[5]。迄今为止，该方法已经成功应用于流体力学[6] [7]，量子力学[8] [9]等领域的微分方程求解。

双曲偏微分方程已被用于描述原子物理不同领域的现象，与普通扩散方程相比，电报方程在模拟反应扩散方程方面更具优势[10]。电报方程又名传输线方程，是描述传输线上任意点电压、电流与传输线一次参数之间关系的微分方程组，用于描述电信号在传输线电缆中的传播和波现象。特别是，它们可以适用于模拟物理学和生物学中反应扩散之间的相互作用。

线性电报方程的标准形式是

$\{\begin{array}{l}\frac{{\partial }^{2}u}{\partial t^2}+a\frac{\partial u}{\partial t}+bu=\frac{{\partial }^{2}u}{\partial x^2}+f\left(x,t\right),\left(x,t\right)\in \Omega \times \left[0,T\right]\\ u\left(x,0\right)=u_0\left(x\right),\frac{\partial u}{\partial t}\left(x,0\right)=u_1\left(x\right),x\in \Omega \\ u\left(x,t\right)=g\left(x,t\right),\left(x,t\right)\in \partial \Omega \times \left[0,T\right]\end{array}$ (1)

其中u是未知函数，T是最终时间。

目前已有许多方法被用来求解电报方程。Zhou和Luo [11]基于Chebyshev多项式建立二维电报方程的Crank-Nicolson配点模型，并讨论了其数值解的存在性、唯一性、稳定性和收敛性。Dehghan和Shokri [12]使用配点法求解二阶双曲电报方程，并使用薄板样条径向基函数直接近似数值解。Wang和Ahmad等[13]是通过将时间变量视为正常空间变量，给出了径向基函数求解二阶一维双曲电报方程的直接无网格法。Ahmad和Seadawy等[14]利用径向基函数的局部无网格法研究三维二阶双曲电报方程，利用显式时间积分技术在时间方向上对模型进行半离散化。Mardani和Hooshmandasl等[15]研究了一类变系数非线性双曲电报方程，在空间域和时域中分别采用移动最小二乘法(MLS)近似和有限差分法(FDM)，给出了该问题的半离散解。Reutskiy和Zhang [16]等将三层的Crank-Nicolson格式应用于时间方向的导数，将电报方程转化为具有混合导数和变系数的二阶椭圆型偏微分方程。朱妍红和宋灵宇等[17]分别使用时空MQ径向基无网格法和时空MQ径向基耦合时空多项式基的无网格法求解电报方程。

本文在此基础上，推导出一维电报方程的Chebyshev谱格式，并对其进行收敛性和稳定性证明。之后再采用Chebyshev-Gauss-Lobatto节点及Chebyshev导数矩阵将电报方程转化为常微分方程，使用Runge-Kutta法对其进行求解，结果表明Chebyshev谱法更有效。

2. 基础知识

设$\Omega \in {{\rm P}}^2$ 是一个有界开域，其边界为$\partial \Omega$ ，$L^2$ 表示在其上定义的所有平方可积函数的集合，其上的内积和范数定义为

$\left(u,v\right)={\displaystyle {\int }_{\Omega }uv\text{d}x\text{d}y},{\Vert u\Vert }_0={\left({\displaystyle {\int }_{\Omega }{\left|u\right|}^2\text{d}x\text{d}y}\right)}^{1/2},\forall u,v\in L^2\left(\Omega \right)$ (2)

定义

$H^m\left(\Omega \right)=\left\{u\in L^2\left(\Omega \right):D^{\alpha }u\in L^2\left(\Omega \right),0\le \left|\alpha \right|\le m\right\}$ (3)

其中$m$ 是非负整数，$\alpha =\left({\alpha }_1,{\alpha }_2\right)$ ，并且${\alpha }_i\ge 0$ 是整数$\left|\alpha \right|={\alpha }_1+{\alpha }_2$ 。

$H^m\left(\Omega \right)$ 上的范数和半范数分别定义为

${\Vert u\Vert }_m={\left({\displaystyle \sum _{0\le \left|\alpha \right|\le m}{\Vert D^{\alpha }u\Vert }_0^2}\right)}^{1/2},{\left|u\right|}_m={\left({\displaystyle \sum _{\left|\alpha \right|=m}{\Vert D^{\alpha }u\Vert }_0^2}\right)}^{1/2}$ (4)

其中$D^{\alpha }u=\frac{{\partial }^{\left|\alpha \right|}u}{\partial x^{{\alpha }_1}\partial y^{{\alpha }_2}}$ 。

设$H_0^m\left(\Omega \right)=\left\{u\in H^m\left(\Omega \right):{D^{\alpha }u\left(x\right)|}_{\partial \Omega }=0,\left|\alpha \right|\le m\right\}$ ，$H^{-m}$ 是$H_0^m\left(\Omega \right)$ 的对偶空间。

此外，令$\omega =:\omega \left(x,y\right)=\omega \left(x\right)\omega \left(y\right)=1/\sqrt{\left(1-x^2\right)\left(1-y^2\right)}$ ，$\Omega ={\left(-1,1\right)}^2$ ，设$L_{\omega }^2\left(\Omega \right)$ 表示其上定义的所有平方可积函数的集合，其上的范数为

${\Vert u\Vert }_{0,\omega }={\left({\displaystyle {\int }_{\Omega }{\left|u\right|}^2\omega dxdy}\right)}^{1/2}$ (5)

并且设$H_{\omega }^m\left(\Omega \right):=\left\{u\in L_{\omega }^2\left(\Omega \right):D^{\alpha }u\in L_{\omega }^2\left(\Omega \right),0\le \left|\alpha \right|\le m\right\}$ 是带有Chebyshev-Gauss-Lobatto正交加权函数的加权Sobolev空间，其上的范数为

${\Vert u\Vert }_{m,\omega }={\left({\displaystyle \sum _{0\le \left|\alpha \right|\le m}{\Vert D^{\alpha }u\Vert }_{0,\omega }^2}\right)}^{1/2}$ (6)

令$H_{0,\omega }^m\left(\Omega \right)=\left\{u\in H_{\omega }^1\left(\Omega \right):{u|}_{\partial \Omega }=0\right\}$ ，${\left(\cdot ,\cdot \right)}_{\omega }$ 表示$L_{\omega }^2\left(\Omega \right)=H_{\omega }^0\left(\Omega \right)$ 的加权互积，设${\Vert \cdot \Vert }_{H^l\left(H_{\omega }^m\right)}$ 是以下空间中的范数

$H^l\left(0,T;H_{\omega }^m\left(\Omega \right)\right)=\left\{v\left(t\right)\in H_{\omega }^m\left(\Omega \right):{\Vert v\Vert }_{H^l\left(H_{\omega }^m\right)}^2={\displaystyle {\int }_0^T{\displaystyle \sum _{i=0}^l{\Vert \frac{d^i}{d{t}^i}v\left(t\right)\Vert }_{m,\omega }^{2}dt<\infty }}\right\}$ (7)

引理2.1 [1] 设$S_N$ 为一维或二维空间中的插值子空间。引入正交映射$P_N:L_{\omega }^2\left(I\right)\to S_N$ ，则对任意的$u\in H_{0,\omega }^1\left(\Omega \right)$ ，都有

$\left(\nabla \left(P_{N}u-u\right),\nabla v\right)=0,\forall v\in S_N$ (8)

引理2.2 [1] 假设$u\in H^m\left(\Omega \right)$ ，对任意的$0\le \mu \le m$ ，存在一个与$u$ 和$N$ 无关的常数$C$ ，使得

${\Vert u-P_{N}u\Vert }_{\mu }\le C{N}^{\mu -m}{\Vert u\Vert }_m,{\Vert P_{N}u\Vert }_{\mu }\le C{\Vert u\Vert }_{\mu }$ (9)

引理2.3 [1] 假设X为Hilbert空间，$B\left(\cdot ,\cdot \right):X\times X\to {\rm P}$ 是一个双线性泛函，且存在常数$\alpha >0,\beta \ge 0$ ，使得

$\left|B\left(u,v\right)\right|\le \alpha \Vert u\Vert \Vert v\Vert ,\forall u,v\in X$ (10)

$\beta {\Vert u\Vert }^2\le B\left(u,u\right),\forall u\in X$ (11)

那么对于X上的任意连续线性泛函F，存在一个唯一的$v\in X$ ，使得

$B\left(u,v\right)=F\left(v\right),\forall v\in X$ (12)

引理2.4 [2] 令${\eta }_N(t)=P_{N}u-u$ ，且$u,u_t,u_{tt}\in L^{\infty }\left(J,H_{\omega }^m\left(\Omega \right)\cap H_{\omega }^m\left(\Omega \right)\right)$ ，对于$t\in \overline{J}$ ，有

$\begin{array}{l}{\rm A}\left({\eta }_N,v\right)=0\\ N\Vert {\eta }_N\left(t\right)\Vert +\left|{\eta }_N\left(t\right)\right|\le C{N}^{1-m}{\Vert u\left(t\right)\Vert }_{L^{\infty }\left(J;H_w^m\left(\Omega \right)\right)}\\ N\Vert {\eta }_{Nt}\left(t\right)\Vert +\left|{\eta }_{Nt}\left(t\right)\right|\le C{N}^{1-m}{\Vert u_t\left(t\right)\Vert }_{L^{\infty }\left(J;H_w^m\left(\Omega \right)\right)}\\ N\Vert {\eta }_{Ntt}\left(t\right)\Vert +\left|{\eta }_{Ntt}\left(t\right)\right|\le C{N}^{1-m}{\Vert u_{tt}\left(t\right)\Vert }_{L^{\infty }\left(J;H_w^m\left(\Omega \right)\right)}\end{array}$ (13)

3. 电报方程的Chebyshev谱法

Chebyshev多项式是在区间[−1, 1]上关于切比雪夫权函数$w\left(x\right)={\left(1-x^2\right)}^{-\frac{1}{2}}$ 的正交多项式系。设$T_k\left(x\right)=\cos \left(k\arccos \left(x\right)\right)$ 为k次的第一类切比雪夫多项式，选取N + 1个Chebyshev-Gauss-Lobatto节点

$x_j=\cos \left(\frac{j\pi }{N}\right),j=0,1,\cdots ,N$ (14)

如果$u\left(x\right)$ 是一个连续可微函数，那么它可以在Chebyshev-Gauss-Lobatto插值点处近似为：

$u_N\left(x\right)={\displaystyle \sum _{k=0}^N{a}_k{L}_k\left(x\right)}$ (15)

$L_k\left(x\right)=\frac{2}{N{\mu }_k}{\displaystyle \sum _{j=0}^N\frac{1}{{\mu }_j}{T}_k}\left(x_j\right)T_k\left(x\right),k=0,1,\cdots ,N,a_k=u\left(x_k\right)$ (16)

其中

${\mu }_j=\{\begin{array}{cc}2,& j=0,N\\ 1,& 1\le j\le N-1\end{array}$

值得注意的是，对于$k,j=0,1,\cdots ,N$ Lagrange插值具有Kronecker内积

$L_k\left(x_j\right)=\{\begin{array}{cc}0,& j\ne k,\\ 1,& j=k.\end{array}$ (17)

为了近似$u\left(x\right)$ 的一阶导数，通过微分插值函数得到

$u_{Nx}\left(x\right)={\displaystyle \sum _{k=0}^{N}u\left(x_k\right)\frac{\text{d}{L}_k\left(x_j\right)}{\text{d}x}}={\displaystyle \sum _{k=0}^N{D}_{kj}u\left(x_k\right)}$ (18)

得到$\left(N+1\right)\times \left(N+1\right)$ 维的导数矩阵$D_N$

$D_N=\{\begin{array}{ll}\frac{2N^2+1}{6}\hfill & k=j=0,\hfill \\ \frac{{\mu }_k}{{\mu }_j}\frac{{\left(-1\right)}^{i+j}}{t_k-t_j}\hfill & k\ne j,\hfill \\ \frac{-x_j}{2\left(1-t_j^2\right)}\hfill & 1\le k=j\le N-1,\hfill \\ -\frac{2N^2+1}{6}\hfill & k=j=N,\hfill \end{array}$ (19)

3.1. 电报方程的谱格式及收敛性分析

在这里，通过使用变换$x^{\prime }=-1+2\left(x-a\right)/\left(x-b\right)$ 确保区间$\left[a,b\right]\leftrightarrow \left[-1,1\right]$ ，电报方程的Chebyshev谱法的半离散格式为

$\left(\frac{{\partial }^2}{\partial t^2}{u}_N,v\right)+a\left(\frac{\partial }{\partial t}{u}_N,v\right)+b\left(u_N,v\right)=\left(\Delta u_N,v\right)+\left(f\left(x,t_n\right),v\right),v\in V_N$ (20)

即

$A\left(u_N,v_N\right)=F\left(v_N\right)$ (21)

其中$u_N\in H_{0,\omega }^1\left(\Omega \right)\cap S_N$ ，使得对任何$v\in S_N$ 成立。$A\left(\cdot ,\cdot \right)$ 是上的有界正定双线性泛函，对于给定的$f\left(t_n\right)$ ，$F(\cdot )$ 是$H_{0,\omega }^1\left(\Omega \right)\cap S_N$ 上的有边界线性泛函。因此，根据Lax-Milgram定理，上式有一个唯一解$u_N\in H_{0,\omega }^1\left(\Omega \right)\cap S_N$ 。

定理1：设$u_N\in H_{0,\omega }^1\left(\Omega \right)\cap S_N$ ，存在一个常数c，使得方程(20)有唯一的解$u_N\left(t\right)$ 满足下列不等式：

${\Vert u_{Nt}\left(t\right)\Vert }^2+{\Vert u_N\left(t\right)\Vert }^2+{\Vert \nabla u_N\left(t\right)\Vert }^2\le {\text{e}}^t\left[\left({\Vert u_{Nt}\left(0\right)\Vert }^2+{\Vert u_N\left(0\right)\Vert }^2+{\Vert \nabla u_N\left(0\right)\Vert }^2\right)+{\displaystyle {\int }_0^t{\Vert f\left(x,s_n\right)\Vert }^2}\text{d}s\right]$ (22)

证明：令(20)式中$v=u_{Nt}$ ，得

$\left(u_{Ntt},u_{Nt}\right)+a\left(u_{Nt},u_{Nt}\right)+b\left(u_N,u_{Nt}\right)+\left(\nabla u_N,\nabla u_{Nt}\right)=\left(f\left(x,t_n\right),u_{Nt}\right)$ (23)

即

$\frac{1}{2}\frac{d}{dt}{\Vert u_{Nt}\Vert }^2+a{\Vert u_{Nt}\Vert }^2+\frac{1}{2}\frac{d}{dt}{\Vert u_N\Vert }^2+\frac{1}{2}\frac{d}{dt}{\Vert \nabla u_{Nt}\Vert }^2=\left(f\left(x,t_n\right),u_{Nt}\right)$ (24)

使用Cauchy-Schwarz不等式并对其在$\left(0,t\right)$ 上进行积分

$\begin{array}{l}\left({\Vert u_{Nt}\left(t\right)\Vert }^2+{\Vert u_N\left(t\right)\Vert }^2+{\Vert \nabla u_N\left(t\right)\Vert }^2\right)-\left({\Vert u_{Nt}\left(0\right)\Vert }^2+{\Vert u_N\left(0\right)\Vert }^2+{\Vert \nabla u_N\left(0\right)\Vert }^2\right)\\ \le {\displaystyle {\int }_0^t{\Vert f\left(x,s_n\right)\Vert }^2+\left(1-2a\right){\Vert u_{Ns}\Vert }^2}\text{d}s\end{array}$ (25)

由Gronwall不等式得

${\Vert u_{Nt}\left(t\right)\Vert }^2+{\Vert u_N\left(t\right)\Vert }^2+{\Vert \nabla u_N\left(t\right)\Vert }^2\le {\text{e}}^t\left[\left({\Vert u_{Nt}\left(0\right)\Vert }^2+{\Vert u_N\left(0\right)\Vert }^2+{\Vert \nabla u_N\left(0\right)\Vert }^2\right)+{\displaystyle {\int }_0^t{\Vert f\left(x,s_n\right)\Vert }^2}\text{d}s\right]$ (26)

定理2：假设$u\left(x,t\right)$ 是问题(1)的解，$u_N\left(x,t\right)$ 是近似方程(20)的解，那么存在一个常数C，使得

$\Vert u_N\left(t\right)-u\left(t\right)\Vert \le C{N}^{-m}\Vert u\Vert +C{N}^{1-m}\Vert u_t\Vert$ (27)

证明：用(20)式减去(1)式，令${\eta }_N\left(t\right)=u\left(t\right)-P_{N}u\left(t\right),e_N\left(t\right)=P_{N}u\left(t\right)-u_N\left(t\right)$ ，得

$\left(e_{Ntt}+{\eta }_{Ntt},v_N\right)+a\left(e_{Nt}+{\eta }_{Nt},v_N\right)+b\left(e_N+{\eta }_N,v_N\right)+\nabla \left(e_N+{\eta }_N,v_N\right)=0$ (28)

则

$\begin{array}{l}\frac{d}{dt}\left(e_{Nt}+{\eta }_{Nt},v_N\right)-\left(e_{Nt}+{\eta }_{Nt},v_{Nt}\right)+a\frac{d}{dt}\left(e_N+{\eta }_N,v_N\right)\\ -a\left(e_N+{\eta }_N,v_{Nt}\right)+b\left(e_N+{\eta }_N,v_N\right)+\nabla \left(e_N+{\eta }_N,v_N\right)=0\end{array}$ (29)

即

$\begin{array}{l}-\left(e_{Nt},v_{Nt}\right)-a\left(e_N,v_{Nt}\right)+b\left(e_N+{\eta }_N,v_N\right)+\nabla \left(e_N+{\eta }_N,v_N\right)\\ =\left({\eta }_{Nt},v_{Nt}\right)-\frac{d}{dt}\left(e_{Nt}+{\eta }_{Nt},v_N\right)-a\frac{d}{dt}\left(e_N+{\eta }_N,v_N\right)+a\left({\eta }_N,v_{Nt}\right)\end{array}$ (30)

令上式中$v_{Nt}=e_N$ ，有

$\begin{array}{l}\left(e_{Nt},e_N\right)+a\left(e_N,e_N\right)-b\left(e_N+{\eta }_N,v_N\right)-\nabla \left(e_N+{\eta }_N,v_N\right)\\ =-\left({\eta }_{Nt},e_N\right)+\frac{d}{dt}\left(e_{Nt}+{\eta }_{Nt},v_N\right)+a\frac{d}{dt}\left(e_N+{\eta }_N,v_N\right)-a\left({\eta }_N,e_N\right)\end{array}$ (31)

即

$\begin{array}{l}\frac{1}{2}\frac{d}{dt}{\Vert e_N\Vert }^2+a{\Vert e_N\Vert }^2-b\left(e_N+{\eta }_N,v_N\right)-\nabla \left(e_N+{\eta }_N,v_N\right)\\ =-\left({\eta }_{Nt},e_N\right)+\frac{d}{dt}\left(e_{Nt}+{\eta }_{Nt},v_N\right)+a\frac{d}{dt}\left(e_N+{\eta }_N,v_N\right)-a\left({\eta }_N,e_N\right)\end{array}$ (32)

对(32)式使用Cauchy-Schwarz不等式得

$\begin{array}{l}\frac{1}{2}\frac{d}{dt}{\Vert e_N\Vert }^2+a{\Vert e_N\Vert }^2-b\left(e_N+{\eta }_N,v_N\right)-\nabla \left(e_N+{\eta }_N,v_N\right)\\ \le \frac{1}{2}{\Vert {\eta }_{Nt}\Vert }^2+\frac{1}{2}{\Vert e_N\Vert }^2+\frac{d}{dt}\left(e_{Nt}+{\eta }_{Nt},v_N\right)+a\frac{d}{dt}\left(e_N+{\eta }_N,v_N\right)+\frac{a}{2}{\Vert {\eta }_N\Vert }^2+\frac{a}{2}{\Vert e_N\Vert }^2\end{array}$ (33)

对(33)式在$\left(0,t\right)$ 上进行积分

$\begin{array}{c}{\Vert e_N\left(t\right)\Vert }^2-{\Vert e_N\left(0\right)\Vert }^2\le {\displaystyle {\int }_0^t{\Vert {\eta }_{Ns}\Vert }^2}+\left(a+b\right){\Vert {\eta }_N\Vert }^2+\left(1-a+b\right){\Vert e_N\Vert }^2+b{\Vert v_N\Vert }^2+{\left|e_N\right|}^2+{\left|{\eta }_N\right|}^2+{\left|v_N\right|}^2\text{d}s\\ +\left(e_{Nt}\left(t\right)+{\eta }_{Nt}\left(t\right),v_N\left(t\right)\right)-\left(e_{Nt}\left(0\right)+{\eta }_{Nt}\left(0\right),v_N\left(0\right)\right)\\ +\left(e_N\left(t\right)+{\eta }_N\left(t\right),v_N\left(t\right)\right)-\left(e_N\left(0\right)+{\eta }_N\left(0\right),v_N\left(0\right)\right)\end{array}$ (34)

其中

$\begin{array}{c}{\displaystyle {\int }_0^t{\Vert e_N\Vert }^2\text{d}s}={\displaystyle {\int }_0^t{\Vert P_{N}u\left(t\right)-u_N\left(t\right)\Vert }^2\text{d}s}\\ ={\displaystyle {\int }_0^t{\Vert P_{N}u\left(t\right)-u\left(t\right)+u\left(t\right)-u_N\left(t\right)\Vert }^2\text{d}s}\\ \le {\displaystyle {\int }_0^t{\Vert P_{N}u\left(t\right)-u\left(t\right)\Vert }^2\text{d}s}+{\displaystyle {\int }_0^t{\Vert u\left(t\right)-u_N\left(t\right)\Vert }^2\text{d}s}\\ \le Ct{N}^{-2m}{\Vert u\Vert }^2\end{array}$ (35)

因为${\Vert v_N\Vert }^2=\left({\displaystyle {\int }_0^{\tau }{e}_N\left(s\right)\text{d}s},{\displaystyle {\int }_0^{\tau }{e}_N\left(s\right)\text{d}s}\right)={\displaystyle {\int }_{\Omega }{\left({\displaystyle {\int }_0^{\tau }{e}_N\left(s\right)\text{d}s}\right)}^2\text{d}x}$

由${\displaystyle {\int }_0^t{\Vert e_N\Vert }^2\text{d}s}$ 的估计可以得到

${\displaystyle {\int }_0^t{\Vert v_N\Vert }^2\text{d}s}\le C{N}^{-2m}{\Vert u\Vert }^2$ (36)

由定义${\left(u-P_{N}u,\varphi \right)}_{\omega }=0,\forall \varphi \in P_N$ ，且$e_N+{\eta }_N=u_N-u$ 知

$\left(e_{Nt}\left(t\right)+{\eta }_{Nt}\left(t\right),v_N\left(t\right)\right)=0,\left(e_N\left(t\right)+{\eta }_N\left(t\right),v_N\left(t\right)\right)=0$ (37)

所以(34)式变为

${\Vert e_N\left(t\right)\Vert }^2-{\Vert e_N\left(0\right)\Vert }^2\le C{N}^{-2m}{\Vert u\Vert }^2+{\displaystyle {\int }_0^t{\Vert {\eta }_{Ns}\Vert }^2}+\left(a+b\right){\Vert {\eta }_N\Vert }^2+\left(1-a+b\right){\Vert e_N\Vert }^2+{\left|e_N\right|}^2+{\left|{\eta }_N\right|}^2\text{d}s$ (38)

由引理2.4得

${\Vert e_N\left(0\right)\Vert }^2={\Vert P_{N}u\left(0\right)-u_N\left(0\right)\Vert }^2\le {\Vert P_{N}u\left(0\right)-u\left(0\right)\Vert }^2+{\Vert u\left(0\right)-u_N\left(0\right)\Vert }^2\le C{N}^{-2m}{\Vert u\Vert }^2$ (39)

综合上式，我们得到

${\Vert e_N\Vert }^2\le C{N}^{-2m}{\Vert u\Vert }^2+C{N}^{2-2m}{\Vert u_t\Vert }^2$ (40)

因此

$\Vert u_N\left(t\right)-u\left(t\right)\Vert \le \Vert e_N\left(t\right)\Vert +\Vert {\eta }_N\left(t\right)\Vert \le C{N}^{-m}\Vert u\Vert +C{N}^{1-m}\Vert u_t\Vert$ (41)

3.2. 电报方程时间方向上的离散

此时，为了保证方程在要求的区域中求解，令$x^{\prime }=\left(a+b\right)/2+x\left(\left(b-a\right)/2\right)$ ，$D^{\prime }=D/\left(\left(b-a\right)/2\right)$ 。对于线性电报方程，通过Chebyshev导数矩阵可以将方程写成矩阵的形式

$\{\begin{array}{l}\frac{\partial u}{\partial t}=v\\ \frac{\partial v}{\partial t}+av+bu=D^{2}u+f\left(x,t\right)\end{array}$ (42)

其中$u={\left(u\left(x_0\right),u\left(x_0\right),\cdots ,u\left(x_0\right)\right)}^{\text{T}}$ ，$x={\left(x_0,x_1,\cdots ,x_N\right)}^{\text{T}}$ ，使用四阶Runge-Kutta法对其进行求解。

4. 数值结果

在这里，通过以下三种误差来衡量计算结果的准确性，

$L_2=\sqrt{{\displaystyle \sum _{i=0}^N{\left|u_i-{\tilde{u}}_i\right|}^2}}$ ，$L_{\infty }=\underset{1\le i\le N}{\max }\left|u_i-{\tilde{u}}_i\right|$ ，$RMS=\sqrt{\frac{1}{N}\left({\displaystyle \sum _{i=1}^N\left|u_i-{\tilde{u}}_i\right|}\right)}$

其中$u_i$ 和${\tilde{u}}_i$ 分别表示为精确解和近似解。

例1 考虑具有初始条件的双曲电报方程

$\frac{{\partial }^{2}u}{\partial t^2}+20\frac{\partial u}{\partial t}+25u=\frac{{\partial }^{2}u}{\partial x^2}+f\left(x,t\right),\left(x,t\right)\in \Omega \times \left[0,T\right]$

$u\left(x,t\right)={\text{e}}^{2t}\sin \left(\pi x\right)$

我们取$\Omega =\left[0,1\right]$ ，T = 1，a = 20，b = 25，边界条件为$u\left(x,t\right)=g\left(x,t\right)$ ，$\left(x,t\right)\in \partial \Omega \times \left[0,1\right]$ ，由精确解可得

$f\left(x,t\right)=\left(69+{\pi }^2\right){\text{e}}^{2t}\sin \left(\pi x\right)$ 。

Figure 1. Numerical solution

图1. 数值解

图1，图2是N = 400，dt = 0.001，T = 1时的数值解和精确解，表1是在dt = 0.001时，不同时间的三种误差，可以发现比文献[15]、[17]中的误差小，表2是误差随Chebyshev节点数量的变化情况，可以发现随着Chebyshev节点的增大，误差在减小。

Figure 2. Exact solution

图2. 精确解

Table 1. dt = 0.001, error of t taking different values

表1. dt = 0.001，t取不同值的误差

t	$L_2$	$L_{\infty }$	RMS
0.5	8.8521e−10	6.1579e−10	4.4316e−11
1	1.8671e−09	1.6175e−09	9.3472e−11
1.5	5.6464e−09	4.3637e−09	2.8268e−10
2	2.0547e−08	1.9698e−08	1.0286e−09
3	1.4408e−07	1.4118e−07	7.2128e−09

Table 2. The error of t = 1 for different values of N

表2. N取不同值时t = 1的误差

N	$L_2$	$L_{\infty }$	RMS
50	8.4442e−07	8.2159e−07	1.2063e−07
100	7.7198e−08	7.0772e−08	7.7587e−09
200	1.1377e−08	8.3661e−09	8.0652e−10
400	1.8671e−09	1.6175e−09	9.3472e−11

例2 考虑非线性电报方程

$\frac{{\partial }^{2}u}{\partial t^2}+2\frac{\partial u}{\partial t}+u^2=\frac{{\partial }^{2}u}{\partial x^2}+f\left(x,t\right),\left(x,t\right)\in \Omega \times \left[0,T\right]$

初始值为

$u\left(x,0\right)=\cosh x,x\in \Omega$

$\frac{\partial u\left(x,0\right)}{\partial t}=-\cosh x,x\in \Omega$

边界条件为

$u\left(-4,t\right)={\text{e}}^{-t}\cosh \left(-4\right)$ ，$u\left(4,t\right)={\text{e}}^{-t}\cosh 4$

精确解为

$u={\text{e}}^{-t}\cosh x$

由精确解可得

$f\left(x,t\right)={\text{e}}^{-2t}{\cosh }^{2}x-2{\text{e}}^{-t}\cosh x$

这里$\Omega \in \left[0,1\right]$ 。

图3和图4是N = 300，dt = 0.001，T = 3时的数值解和精确解，表3是在dt = 0.001时，不同时间的三种误差，可以发现比文献[15]中的误差小，表4是误差随着Chebyshev节点数量的变化情况，可以发现随着Chebyshev节点的增大，误差在减小。

Figure 3. Numerical solution

图3. 数值解

Figure 4. Exact solution

图4. 精确解

Table 3. dt = 0.001, error of t taking different values

表3. dt = 0.001，t取不同值的误差

t	$L_2$	$L_{\infty }$	RMS
0.5	3.6898e−05	3.0630e−05	5.1668e−06
1	3.0670e−05	2.2985e−05	4.2947e−06
1.5	1.1252e−05	1.0758e−05	1.5756e−06
3	9.7158e−06	7.2337e−06	1.3605e−06
5	1.2927e−06	1.0662e−06	1.8102e−07

Table 4. The error of t = 1 for different values of N

表4. N取不同值时t = 1的误差

N	$L_2$	$L_{\infty }$	RMS
50	3.0670e−05	2.2985e−05	4.2947e−06
100	6.3692e−06	5.3988e−06	6.3378e−07
200	2.6677e−06	1.8957e−06	1.5376e−07

5. 本章小结

本章提出了一种基于Chebyshev-Gauss-Lobatto伪谱格式的求解非线性电报方程的有效数值解，并分析了Chebyshev谱法的收敛性和稳定性，将Chebyshev谱法得到的误差与其他方法进行比较，发现该方法误差更小，相较于其他的方法具有的精度更高。

基金项目

国家自然科学基金青年项目(编号：12101482)，中国博士后科学基金面上项目(2022M722604)，陕西数理基础科学研究项目(23JSQ042)，甘肃省科技计划项目(23CXGL0018)，西安市科技局高校院所科技人员服务企业项目(24GXFW0038)。

NOTES

^*通讯作者。