粘弹性方程的时间间断Galerkin时空有限元方法

A Time Discontinuous Galerkin Space-Time Finite Element Method for the Viscoelastic Equation

摘要: 粘弹性方程是研究地震波、断裂力学问题的主要数学物理方程之一。为时空统一高精度求解位移$u$ 和速度$u_t$ 的近似解，针对方程特性，引入中间变量$\sigma =u_t$ 将原问题转换为方程组系统，并建立时间间断时空有限元格式去获得位移和速度的近似解。本文采用有限元分析，结合基于Radau积分节点的Lagrange插值，推导出位移关于$L^{\infty }\left(\left[0,T\right];H^1\left(\Omega \right)\right)$ –模和速度关于$L^{\infty }\left(\left[0,T\right];L^2\left(\Omega \right)\right)$ –模最优误差估计。最后，通过一个数值例子证明了方法的可行性和误差分析结果的合理性。

Abstract: The viscoelastic equation is one of the main mathematical and physical equations for studying seismic waves and fracture mechanics problems. To uniformly and precisely solve for the approximate solutions of displacement $u$ and velocity $u_t$ in space-time, by analyzing the equation’s properties, intermediate variables $\sigma =u_t$ are introduced to reformulate the original problem as a system of equations. Then, a temporally discontinuous space-time finite element framework is established to obtain the approximate solutions of both displacement and velocity fields. By using finite element analysis and combining with Lagrange interpolation based on Radau integration nodes, the optimal error estimation of displacement with respect to the $L^{\infty }\left(\left[0,T\right];H^1\left(\Omega \right)\right)$ -norm and velocity with respect to the $L^{\infty }\left(\left[0,T\right];L^2\left(\Omega \right)\right)$ -norm is derived. Finally, a numerical example was provided to demonstrate the feasibility of the method and the rationality of the theoretical results of the error analysis.

1. 引言

研究下述的粘弹性方程初边值问题

$\{\begin{array}{l}{u}_{tt}-\epsilon \Delta u_t-\gamma \Delta u=f\left(x,t\right)\left(x,t\right)\in \Omega \times J\\ u\left(x,t\right)=0\left(x,t\right)\in \partial \Omega \times J\\ u\left(x,0\right)={\Phi }_1\left(x\right),u_t\left(x,0\right)={\Phi }_2\left(x\right)x\in \Omega \end{array}$ (1)

其中$\Omega \subset R^2$ 是凸多边形有界区域，边界为$\partial \Omega$ ，$J=\left[0,T\right]$ ，$T\left(0<T<\infty \right)$ 是总时间。$u\left(x,t\right)$ 是未知函数，$u_{tt}=\frac{{\partial }^{2}u}{\partial t^2}$ ，$u_t=\frac{\partial u}{\partial t}$ ，$\gamma$ 和$\epsilon$ 分别是弹性系数和粘性系数，$f\left(x,t\right)$ 是源函数，${\Phi }_1\left(x\right)$ 和${\Phi }_2\left(x\right)$ 是足够光滑的已给定初值函数。

粘弹性方程是研究地震波、断裂力学问题的主要数学物理方程之一[1]。在实际问题中，该方程其本身较为复杂而难以获得解析解，因此求其数值解来研究该方程成为一种可行的途径[2]。文献[3]针对一类半线性粘弹性方程采用非标准Hermite元推导出$H^1\left(\Omega \right)$ –模意义下的超收敛性质。文献[4]建立了一种非协调变网格有限元方法求解粘弹性方程问题。文献[5] [6]针对粘弹性方程分别运用半离散$H^1$ -Galerkin混合有限元方法和标准有限元方法进行数值模拟研究，推导出了最优收敛阶的误差估计结果。文献[7]应用分裂正定混合有限元方法推导出相应的误差估计，并利用数值实验证明了理论结果。Wang等[8]用不连续分段多项式的离散弱梯度算子构造了粘弹性方程的弱Galerkin有限元格式，并证明了格式的稳定性和$H^1$ –模最优误差估计。为了降低计算成本、提高计算速度，Luo等[9]提出了粘弹性方程的POD基降维有限元方法，并给出了降维格式误差估计。然而，在粘弹性方程的数值求解中，上述方法均采用有限差分技术对时间导数进行离散化处理，且该离散化过程与空间导数的离散化是独立进行的。这种处理方法难于获得时空统一高精度数值解。

时空有限元(TSFE)方法作为一种能够有效统一时空离散的高精度数值方法，在Nickell Seckman等人的文章中最早见其雏形，并不断延伸应用于流体力学与热传导等领域[10] [11]。文献[12]针对非线性薛定谔方程采用时间间断时空有限元法(TDG-TSFEM)进行求解，第一次提出将有限元与Lagrange插值相结合，证明了格式的最优收敛性。文献[13]为线性时间间断时空有限法的推导和实施过程提出一种通用的计算框架。文献[14]为了提高求解非线性Volterra积分微分方程的间断Galerkin (DG)时间步进法的全局精度，提出了一种简洁高效的后处理技术。这些研究表明，时间间断时空有限元法具有较小的数值耗散和非物理振荡，可以实现时空统一高精度等优势。然而，据我们所知，目前还没有该方法对粘弹性方程初边值问题的研究报道。为此，本文针对方程特性引入中间变量$\sigma =u_t$ 将原问题转化为方程组系统，并构造TDG-TSFE格式同时求解位移$u$ 和速度$u_t$ 的时空高精度数值解。进一步，通过引入相关算子和投影，采用有限元分析并结合基于Radau积分节点的Lagrange插值，推导出位移关于$L^{\infty }\left(\left[0,T\right];H^1\left(\Omega \right)\right)$ –模和速度关于$L^{\infty }\left(\left[0,T\right];L^2\left(\Omega \right)\right)$ –模最优误差估计。最后，通过一个数值例子证明了该方法的可行性和误差分析结果的合理性。

本文涉及到的$c_i$ ，$C_i$ 均是正常数，不依赖于时间和空间步长，且在不同位置可能具有不同取值。

2. 预备知识及数值格式

本文中使用标准Sobolev空间$H^m\left(\Omega \right)$ ，其范数表示为${\Vert v\Vert }_m$ ；$L^2\left(\Omega \right)=H^0\left(\Omega \right)$ 是Lebesgue空间，$\left(v,w\right)$ 表示其内积，$\Vert v\Vert$ 表示其范数，并且${\Vert \cdot \Vert }_{0,\infty }$ 表示$L^{\infty }\left(\Omega \right)$ 空间的范数。设$l$ ，$m$ 都是非负整数，分别给出时空域上Sobolev空间$H^l\left(J;H^m\left(\Omega \right)\right)$ 和$L^{\infty }\left(J;H^m\left(\Omega \right)\right)$ 的定义为

$H^l\left(J;H^m\left(\Omega \right)\right)=\left\{u\in H^m\left(\Omega \right):{\displaystyle {\int }_J{\displaystyle \sum _{i=0}^l{\Vert \frac{{\text{d}}^i}{\text{d}{t}^i}u\left(x,t\right)\Vert }_m^2}\text{d}x<\infty }\right\}$ ,

$L^{\infty }\left(J;H^m\left(\Omega \right)\right)=\left\{u\in H^m\left(\Omega \right):\text{ess}\underset{J}{\sup }{\Vert u\left(x,t\right)\Vert }_m<\infty \right\}$ ,

并分别赋予如下范数

${\Vert u\Vert }_{H^l\left(J;H^m\left(\Omega \right)\right)}={\left({\displaystyle {\int }_J{\displaystyle \sum _{i=0}^l{\Vert \frac{{\text{d}}^i}{\text{d}{t}^i}u\left(x,t\right)\Vert }_m^2}\text{d}x}\right)}^{\frac{1}{2}}$ , ${\Vert u\Vert }_{L^{\infty }\left(J;H^m\left(\Omega \right)\right)}=\underset{J}{\max }{\Vert u\left(x,t\right)\Vert }_m$ .

特别地，当$J=I_n$ 时，简记${\Vert u\Vert }_{L^2\left(I_n;H^m\left(\Omega \right)\right)}={\left|\Vert u\Vert \right|}_{m,n}$ ；当$m=0$ 时，简记${\left|\Vert u\Vert \right|}_{m,n}={\left|\Vert u\Vert \right|}_n$ ，即${\left|\Vert u\Vert \right|}_n:={\left({\displaystyle {\int }_J{\Vert u\left(x,t\right)\Vert }^2\text{d}t}\right)}^{1/2}$ 。在有限维Hilbert向量空间${\left(X\right)}^q$ 上定义内积

$\left(\left(\Phi ,\Psi \right)\right)={\displaystyle \sum _{i=1}^q\left({\varphi }_i,{\psi }_i\right)}$ ，其中$\Phi ={\left({\varphi }_1,{\varphi }_2,\cdots ,{\varphi }_q\right)}^{\text{T}}$ ，$\Psi ={\left({\psi }_1,{\psi }_2,\cdots ,{\psi }_q\right)}^{\text{T}}\in {\left(X\right)}^q$ ，并定义相应的范数${\left({\displaystyle \sum _{i=1}^q{\Vert {\psi }_i\Vert }^2}\right)}^{1/2}$ ，记为$\left|\Vert \cdot \Vert \right|$ 。

为得到方程(1)的TDG-TSFE数值格式，将区间$J=\left[0,T\right]$ 划分。设$0=t^0<t^1<\cdots <t^N=T$ ，$I_n=\left(t^n,t^{n+1}\right]$ ，$k_n$ 是时间步长，$k_n=t^{n+1}-t^n$ ，$k=\underset{n}{\max }{k}_n$ ，$\left(n=0,1,\cdots ,N-1\right)$ 。定义时空片$Q_n:=\Omega \times I_n$ ，设$T_h^n=\left\{K\right\}$ 是空间$\Omega$ 的一种拟一致剖分，$h_K$ 表示单元$K$ 的长度，$h_n=\underset{K\in T_h^n}{\max }{h}_k$ ，$h=\underset{n}{\max }{h}_n$ ，且在不同$Q_n$ 中允许空间$\Omega$ 采用不同的剖分。建立有限元空间

$S_h^n=\left\{\chi \in H_0^1\left(\Omega \right):{\chi |}_K\in P_r\left(K\right),K\in T_h^n\right\}$ ,

其中$P_p\left(K\right)$ 是由在单元$K$ 上变量$x$ 的次数最高不超过$p$ 次的多项式函数构成的。对应$t^0$ 的空间记为$S_h^{-1}$ ，为了方便，置$S_h^{-1}=S_h^0$ 。在下文中，对于$s,m=0,1,\cdots$ 和$u\in H^m\left(\Omega \right)$ ，设${\Vert h_n^{s}u\Vert }_{m,h}:={\left({\displaystyle \sum _{K\in T_h^n}{h}_K^{2s}{\Vert u\Vert }_{m,K}^2}\right)}^{1/2}$ 。

设$q>0$ 是整数。令$V_{hk}^n=V_{hk}^n\left(q\right)$ 表示由分块多项式$\phi :\Omega \times \left(0,T\right]\to R$ 组成的空间，即

$V_{hk}=\left\{\phi :{\phi |}_{\Omega \times I_n}={\displaystyle \sum _{j=0}^{q-1}{t}^j{\chi }_j\left(x\right),{\chi }_j\left(x\right)\in S_h^n}\right\}$ ,

由$V_{hk}$ 定义可知，对于$\forall x\in \Omega$ ，函数$\phi$ 是关于$t$ 的$q-1$ 次分片多项式，在$t^n\left(n=1,2,\cdots ,N-1\right)$ 处允许间断。而对$\forall t\in I_n$ ，函数$\phi$ 是关于$x$ 的次数最高不超过$r$ 的连续多项式。进而令

$V_{hk}^n=\left\{{\phi |}_{\Omega \times I_n}:\phi \in V_{hk}\right\}$ ,

则由$V_{hk}^n$ 定义可知，${\phi }^n={\phi }^{n-}$ ，其中${\phi }^{n-}=\underset{t\to t^{n-}}{\lim }\phi \left(t\right)$ ，${\phi }^{n+}=\underset{t\to t^{n+}}{\lim }\phi \left(t\right)$ 。

假设问题(1)在$\left[0,T\right]$ 上存在唯一的光滑解$u$ 。引入辅助变量$\sigma =u_t$ 得以下方程组

$\{\begin{array}{l}\left(\text{a}\right){\sigma }_t-\epsilon \Delta \sigma -\gamma \Delta u=f,\left(x,t\right)\in \Omega \times \left[0,T\right]\\ \left(\text{b}\right)u_t=\sigma ,\left(x,t\right)\in \Omega \times \left[0,T\right]\\ \left(\text{c}\right)u(x,t)=0,\left(x,t\right)\in \partial \Omega \times \left[0,T\right]\\ \left(\text{d}\right)u\left(x,0\right)={\Phi }_1\left(x\right),\sigma \left(x,0\right)={\Phi }_2\left(x\right),x\in \Omega \end{array}$ (2)

在(2)(a)两端乘以$v\in H^1\left(I_n;H_0^1\left(\Omega \right)\right)$ ，在(2)(b)两端乘以$w\in H^1\left(I_n;H_0^1\left(\Omega \right)\right)$ ，并在$Q_n$ 上积分。利用Green公式可得方程组(2)的弱形式为：求$\left\{u,\sigma \right\}\in L^2\left(I_n;H_0^1\left(\Omega \right)\right)\times L^2\left(I_n;H_0^1\left(\Omega \right)\right)$ 使得

$\{\begin{array}{l}\left(\text{a}\right)\left({\sigma }^{n+1},v^{n+1}\right)-{\displaystyle {\int }_{I_n}\left(\sigma ,v_t\right)\text{d}t}+\epsilon {\displaystyle {\int }_{I_n}\left(\nabla \sigma ,\nabla v\right)\text{d}t}+\gamma {\displaystyle {\int }_{I_n}\left(\nabla u,\nabla v\right)\text{d}t}=\left({\sigma }^n,v^{n+}\right)+{\displaystyle {\int }_{I_n}\left(f,v\right)\text{d}t}\\ \left(\text{b}\right)\left(u^{n+1},w^{n+1}\right)-{\displaystyle {\int }_{I_n}\left(u,w_t\right)\text{d}t}-{\displaystyle {\int }_{I_n}\left(\sigma ,w\right)\text{d}t}=\left(u^n,w^{n+}\right)\end{array}$ (3)

通过将其有限元近似函数代入弱形式(3)，构造相应的TDG-TSFE数值格式为：求$\left\{U,Z\right\}\in V_{hk}\times V_{hk}$ ，使得对$\forall v\in V_{hk}^n$ ，$\forall w\in V_{hk}^n\left(n=0,1,\cdots ,N-1\right)$ ，满足

$\{\begin{array}{l}\left(\text{a}\right)\left(Z^{n+1},v^{n+1}\right)-{\displaystyle {\int }_{I_n}\left(Z,v_t\right)\text{d}t}+\epsilon {\displaystyle {\int }_{I_n}\left(\nabla Z,\nabla v\right)\text{d}t}+\gamma {\displaystyle {\int }_{I_n}\left(\nabla U,\nabla v\right)\text{d}t}=\left(Z^n,v^{n+}\right)+{\displaystyle {\int }_{I_n}\left(f,v\right)\text{d}t}\\ \left(\text{b}\right)\left(U^{n+1},w^{n+1}\right)-{\displaystyle {\int }_{I_n}\left(U,w_t\right)\text{d}t}-{\displaystyle {\int }_{I_n}\left(Z,w\right)\text{d}t}=\left(U^n,w^{n+}\right)\end{array}$ (4)

3. 辅助引理

本文采用基于Radau积分公式节点的Lagrange插值与有限元分析结合的技巧，分析TDG-TSFE法的最优误差估计，而Radau求积方法适合近似计算定义在半开半闭区间上的积分，与TDG-TSFEM基函数的时间剖分节点处允许间断的特点相符合。为此，现引入本文所涉及的符号和定义，简要介绍Radau求积公式及其相关的基本引理。

Radau求积公式[12]：对于任意的$q\ge 1$ ，如果求积公式

${\displaystyle {\int }_0^{1}g\left(\tau \right)\text{d}\tau \cong {\displaystyle \sum _{j=1}^q{\omega }_{j}g}}\left({\tau }_j\right),0<{\tau }_1<{\tau }_2<\cdots <{\tau }_q=1$ .

对于所有次数$\le 2q-2$ 次的多项式都是准确成立的，称其为Radau求积公式，其中权函数${\omega }_j={\displaystyle {\int }_0^1{l}_j^2\left(t\right)\text{d}t}$ 。

对任意的$q\ge 1$ ，设${\left\{l_i\left(\tau \right)\right\}}_{i=1}^q$ 是以Radau点${\tau }_1,{\tau }_2,\cdots ,{\tau }_q$ 为节点的$q-1$ 次Lagrange插值多项式，即

$l_i\left(\tau \right)={\displaystyle \prod _{j=1,i\ne j}^q\frac{\tau -{\tau }_j}{{\tau }_i-{\tau }_j}}$ ，利用$t=t^n+\tau k_n$ 将[0, 1]映射到$\overline{I_n}$ ，

$\begin{array}{l}{t}^{n,j}=t^n+{\tau }_j{k}_n,j=1,2,\cdots ,q,\\ l_{n,i}\left(t\right)=l_i\left(\tau \right),t=t^n+\tau k_n,t^{n,q}=t^{n+1},\\ {\omega }_{n,i}={\displaystyle {\int }_{t^n}^{t^{n+1}}{l}_{n,i}\left(t\right)\text{d}t}=k_n{\displaystyle {\int }_0^1{l}_i\left(\tau \right)\text{d}\tau }=k_n{\omega }_i\end{array}$ (5)

设$M,N$ 分别是$q\times q$ 阶的矩阵，$N_{ij}={\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right)={\omega }_j{l^{\prime }}_i\left({\tau }_j\right)$ ，$M=e_q{e}_q^{\text{T}}-N$ ，$e_q^{\text{T}}=\left(0,\cdots ,1\right)\in R^q$ ，其中$M,N$ 与$k_n$ 无关。若$y={\left(y^{n,1},\cdots ,y^{n,q}\right)}^{\text{T}}\in R^q$ ，有

$y^{\text{T}}My={\displaystyle \sum _{i=1}^q{\delta }_{qi}}{y}^{n,q}{y}^{n,i}-{\displaystyle \sum _{i,j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right)y^{n,j}{y}^{n,i}}$ .

引理3.1 [12] 设对角阵$D=diag\left\{{\tau }_1,{\tau }_2,\cdots ,{\tau }_q\right\}$ ，则矩阵$\tilde{M}=D^{-1/2}M{D}^{1/2}$ 正定，即存在$\lambda :=\frac{1}{2}\min \left\{\frac{{\omega }_1}{{\tau }_1},\cdots ,\frac{{\omega }_{q-1}}{{\tau }_{q-1}},1+{\omega }_q\right\}>0$ ，使得

$x^{\text{T}}\tilde{M}x\ge \lambda {\left|x\right|}^2=\lambda \left({\displaystyle \sum _{i=1}^q{x}_i^2}\right),\forall x={\left(x_1,x_2,\cdots ,x_q\right)}^{\text{T}}\in R^q$. (6)

由于${\left\{l_{n,i}\left(t\right)\right\}}_{i=1}^q$ 是多项式空间$P_{q-1}\left(I_n\right)$ 的一组基函数，故在$I_n$ 内有

$U\left(x,t\right)={\displaystyle \sum _{j=1}^q{l}_{n,j}{U}^{n,j}},Z\left(x,t\right)={\displaystyle \sum _{j=1}^q{l}_{n,j}{Z}^{n,j}}$ , (7)

其中$U^{n,j}=U^{n,j}\left(x\right)\in S_h^n$ ，$Z^{n,j}=Z^{n,j}\left(x\right)\in S_h^n$ 。

将(7)代入格式(4)，并在格式(4)中分别取$v=l_{n,i}\left(t\right)\psi$ 和$w=l_{n,i}\left(t\right)\varphi$ $\left(\psi ,\varphi \in S_h^n\right)$ ，再采用Radau积分公式离散格式(4)中时间积分项，得到格式(4)的Radau积分公式离散的等价形式(Ⅰ)

$\{\begin{array}{l}\left(\text{a}\right){\delta }_{qi}\left(Z^{n,q},\psi \right)-{\displaystyle \sum _{j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right)}\left(Z^{n,j},\psi \right)+\epsilon k_n{\omega }_i\left(\nabla Z^{n,i},\nabla \psi \right)+\gamma k_n{\omega }_i\left(\nabla U^{n,i},\nabla \psi \right)\\ =l_{n,i}\left(t^n\right)\left(Z^n,\psi \right)+{\displaystyle {\int }_{I_n}\left(l_{n,i}\left(t\right)f,\psi \right)\text{d}t}\\ \left(\text{b}\right){\delta }_{qi}\left(U^{n,q},\varphi \right)-{\displaystyle \sum _{j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right)}\left(U^{n,j},\varphi \right)-k_n{\omega }_i\left(Z^{n,i},\varphi \right)=l_{n,i}\left(t^n\right)\left(U^n,\varphi \right)\end{array}$ (8)

进一步，为了充分利用矩阵$\tilde{M}$ 的性质，记${\tilde{U}}^{n,j}={\tau }_j^{-1/2}{U}^{n,j}\left(x\right)\in S_h^n$ ，${\tilde{Z}}^{n,j}={\tau }_j^{-1/2}{Z}^{n,j}\left(x\right)\in S_h^n$ 。则${Z|}_{I_n}$ 和${U|}_{I_n}$ 可表示为$U\left(x,t\right)={\displaystyle \sum _{j=1}^q{\tau }_j^{1/2}{l}_{n,j}{\tilde{U}}^{n,j}}$ ，$Z\left(x,t\right)={\displaystyle \sum _{j=1}^q{\tau }_j^{1/2}{l}_{n,j}{\tilde{Z}}^{n,j}}$ ，将该表达式代入格式(4)，同时在格式(4)分别取$v={\tau }_i^{-1/2}{l}_{n,i}\left(t\right)\psi$ 和$w={\tau }_i^{-1/2}{l}_{n,i}\left(t\right)\varphi$ $\left(\psi ,\varphi \in S_h^n\right)$ ，得Radau积分公式离散的等价形式(Ⅱ)

$\{\begin{array}{l}\left(\text{a}\right){\delta }_{qi}\left({\tilde{Z}}^{n,q},\psi \right)-{\displaystyle \sum _{j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right){\tau }_j^{1/2}{\tau }_i^{-1/2}\left({\tilde{Z}}^{n,j},\psi \right)}+\epsilon k_n{\omega }_i\left(\nabla {\tilde{Z}}^{n,i},\nabla \psi \right)+\gamma k_n{\omega }_i\left(\nabla {\tilde{U}}^{n,i},\nabla \psi \right)\\ ={\tau }_i^{-1/2}\left(l_{n,i}\left(t^n\right)\left(Z^n,\psi \right)+{\displaystyle {\int }_{I_n}\left(l_{n,i}f,\psi \right)\text{d}t}\right)\\ \left(\text{b}\right){\delta }_{qi}\left({\tilde{U}}^{n,q},\varphi \right)-{\displaystyle \sum _{j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right){\tau }_j^{1/2}{\tau }_i^{-1/2}\left({\tilde{U}}^{n,j},\varphi \right)}-k_n{\omega }_i\left({\tilde{Z}}^{n,i},\varphi \right)={\tau }_i^{-1/2}{l}_{n,i}\left(t^n\right)\left(U^n,\varphi \right)\end{array}$ (9)

定义3.1 [15] 定义椭圆投影算子$P_h^n:H_0^1\to S_h^n$ 使得

$\left(\nabla \left(u-P_h^{n}u\right),\nabla \phi \right)=0,\forall \phi \in S_h^n$ , (10)

则当$u\in H^s\left(\Omega \right)\cap H_0^1\left(\Omega \right)\left(2\le s\le r+1\right)$ 时，$P_h^n$ 满足如下估计式

$\Vert u-P_h^{n}u\Vert \le C{\Vert h_n^{s}u\Vert }_{s,h},{\Vert u-P_h^{n}u\Vert }_1\le C{\Vert h_n^{s-1}u\Vert }_{s,h}$ . (11)

定义3.2 [12] 定义Lagrange插值算子${{\rm I}}_{q-1}^n:C\left(I_n\right)\to P_{q-1}\left(I_n\right)$ ，满足

$\left({{\rm I}}_{q-1}^{n}y\right)\left(t^{n,j}\right)=y\left(t^{n,j}\right),j=1,\cdots ,q$ (12)

其中$t^{n,j}$ 定义在(5)。注意到${{\rm I}}_{q-1}^{n}u\left(x,\cdot \right)\in P_{q-1}\left({{\rm I}}_n\right)$ ，${{\rm I}}_{q-1}^{n}u\left(x,t^{n+1}\right)=u\left(x,t^{n+1}\right),x\in \Omega$ 。

定义函数$W_0:\left(0,T\right]\to H_0^1\left(\Omega \right)$ 和$W_1:\left(0,T\right]\to H_0^1\left(\Omega \right)$ ，满足

$W_0\left(x,t\right)={{\rm I}}_{q-1}^n{P}_h^{n}u\left(x,t\right),W_1\left(x,t\right)={{\rm I}}_{q-1}^n{P}_h^n\sigma \left(x,t\right),\left(x,t\right)\in \Omega \times {{\rm I}}_n$ , (13)

其中$W_0\left(x,t\right)$ ，$W_1\left(x,t\right)$ 均是$V_{hk}$ 的元素。在$I_n$ 上，仍置${W_0|}_{I_n}=W_0,{W_1|}_{I_n}=W_1$ 。为了进行误差分析，将$U-u$ 和$Z-\sigma$ 分解

$U-u=E_0-{\rho }_0,E_0=U-W_0,{\rho }_0=u-W_0$ ,

$Z-\sigma =E_1-{\rho }_1,E_1=Z-W_1,{\rho }_1=\sigma -W_1$ , (14)

引理3.2 [15] 对$2\le s\le r+1$ 和$\kappa =0,1$ ，误差${\rho }_0$ 和${\rho }_1$ 有估计式

$\begin{array}{l}\left(\text{a}\right){\left|\Vert {\rho }_0\Vert \right|}_{\kappa ,n}\le c{k}_n^q{\left|\Vert u^{\left(q\right)}\Vert \right|}_{\kappa ,n}+c{k}_n^{1/2}\underset{I_n}{\max }{\Vert h_n^{s-\kappa }u\Vert }_{s,h},\\ \left(\text{b}\right){\left|\Vert {\rho }_1\Vert \right|}_{\kappa ,n}\le c{k}_n^q{\left|\Vert {\sigma }^{\left(q\right)}\Vert \right|}_{\kappa ,n}+c{k}_n^{1/2}\underset{I_n}{\max }{\Vert h_n^{s-\kappa }\sigma \Vert }_{s,h},\\ \left(\text{c}\right)\underset{I_n}{\max }{\Vert {\rho }_0\Vert }_{\kappa }\le c{k}_n^q\underset{I_n}{\max }{\Vert u^{\left(q\right)}\Vert }_{\kappa }+c\underset{I_n}{\max }{\Vert h_n^{s-\kappa }u\Vert }_{s,h},\\ \left(\text{d}\right)\underset{I_n}{\max }{\Vert {\rho }_1\Vert }_{\kappa }\le c{k}_n^q\underset{I_n}{\max }{\Vert {\sigma }^{\left(q\right)}\Vert }_{\kappa }+c\underset{I_n}{\max }{\Vert h_n^{s-\kappa }\sigma \Vert }_{s,h}.\end{array}$ (15)

将(14)代到(4)得到${E_0|}_{I_n}$ 和${E_1|}_{I_n}$ 满足的方程组

$\{\begin{array}{l}\left(\text{a}\right)\left(E_1^{n+1},v^{n+1}\right)-{\displaystyle {\int }_{I_n}\left(E_1,v_t\right)\text{d}t}+\epsilon {\displaystyle {\int }_{I_n}\left(\nabla E_1,\nabla v\right)\text{d}t}+\gamma {\displaystyle {\int }_{I_n}\left(\nabla E_0,\nabla v\right)\text{d}t}\\ =\left(E_1^n,v^{n+}\right)+{\displaystyle {\int }_{I_n}\left(f,v\right)\text{d}t}-\{\left(W_1^{n+1},v^{n+1}\right)-{\displaystyle {\int }_{I_n}\left(W_1,v_t\right)\text{d}t}\\ +\epsilon {\displaystyle {\int }_{I_n}\left(\nabla W_1,\nabla v\right)\text{d}t}+\gamma {\displaystyle {\int }_{I_n}\left(\nabla W_0,\nabla v\right)\text{d}t}-\left(W_1^n,v^{n+}\right)\}\\ \left(\text{b}\right)\left(E_0^{n+1},w^{n+1}\right)-{\displaystyle {\int }_{I_n}\left(E_0,w_t\right)\text{d}t}-{\displaystyle {\int }_{I_n}\left(E_1,w\right)\text{d}t}\\ =\left(E_0^n,w^{n+}\right)-\{\left(W_0^{n+1},w^{n+1}\right)-{\displaystyle {\int }_{I_n}\left(W_0,w_t\right)\text{d}t}-{\displaystyle {\int }_{I_n}\left(W_1,w\right)\text{d}t}-\left(W_0^n,w^{n+}\right)\}\end{array}$ (16)

这里取初始值$W_0^0=P_h^0{\Phi }_1\left(x\right),W_1^0=P_h^0{\Phi }_2\left(x\right),U^0=P_h^0{\Phi }_1\left(x\right),Z^0=P_h^0{\Phi }_2\left(x\right)$ ，故$E_0^0=U^0-W_0^0=0$ ，$E_1^0=0$ 。另外，引入函数$w_0,w_1,{\eta }_0,{\eta }_1$ ，使得

$w_0=P_h^{n}u\left(x,t\right),{\eta }_0=u-w_0,\left(x,t\right)\in \Omega \times I_n$ ,

$w_1=P_h^n\sigma \left(x,t\right),{\eta }_1=\sigma -w_1,\left(x,t\right)\in \Omega \times I_n\left(n=0,\cdots ,N-1\right)$ .

显然，这些函数在区间$I_n$ 内关于时间t连续，并且只有当$S_h^{n-1}\ne S_h^n$ 时，函数在处$t^n$ 点间断。因此${W_{\kappa }|}_{I_n}$ $\left(\kappa =0,1\right)$ 是$w_{\kappa }$ 在点$t^{n,j}$ 处关于t的插值，即

$W_{\kappa }\left(x,t\right)={\displaystyle \sum _{j=1}^q{l}_{n,j}\left(t\right)w_{\kappa }^{n,j}}\left(x\right),E_{\kappa }\left(x,t\right)={\displaystyle \sum _{j=1}^q{l}_{n,j}\left(t\right)E_{\kappa }^{n,j}}\left(x\right),\left(x,t\right)\in \Omega \times I_n$ . (17)

在(16)中分别取$v=l_{n,i}\left(t\right)\psi$ 和$w=l_{n,i}\left(t\right)\Delta \varphi$ $\left(\psi ,\varphi \in S_h^n\right)$ ，并将(17)代入(16)得

$\begin{array}{l}\left(\text{a}\right){\delta }_{qi}\left(E_1^{n,q},\psi \right)-{\displaystyle \sum _{j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right)}\left(E_1^{n,j},\psi \right)+\epsilon k_n{\omega }_i\left(\nabla E_1^{n,i},\nabla \psi \right)+\gamma k_n{\omega }_i\left(\nabla E_0^{n,i},\nabla \psi \right)\\ =l_{n,i}\left(t^n\right)\left(E_1^n,\psi \right)-l_{n,i}\left(t^n\right)\left(J\left[{\eta }_1^n\right],\psi \right)+\left({\theta }_1^{n,i}+A_1^{n,i}+\gamma B_1^{n,i}-\epsilon B_0^{n,i},\psi \right)\end{array}$

$\begin{array}{l}\left(\text{b}\right){\delta }_{qi}\left(\nabla E_0^{n,q},\nabla \varphi \right)-{\displaystyle \sum _{j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right)}\left(\nabla E_0^{n,j},\nabla \varphi \right)-k_n{\omega }_i\left(\nabla E_1^{n,i},\varphi \right)\\ =l_{n,i}\left(t^n\right)\left(\nabla E_0^n,\nabla \varphi \right)+\left(A_0^{n,i},\nabla \varphi \right)+\left(B_0^{n,i},\varphi \right)\end{array}$ (18)

其中${\eta }_{\kappa }^{n,i}={\eta }_{\kappa }\left(\cdot ,t^{n,i}\right),\left(\kappa =0,1\right)$ ，

$\begin{array}{l}{\theta }_1^{n,i}:={\delta }_{qi}{\eta }_1^{n,q}-{\displaystyle \sum _{j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right){\eta }_1^{n,j}}-l_{n,i}\left(t^n\right){\eta }_1^{n+},\\ A_1^{n,i}:={\displaystyle \sum _{j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right){\sigma }^{n,j}}-{\displaystyle {\int }_{I_n}{l^{\prime }}_{n,i}\sigma \text{d}t,}{B}_1^{n,i}:=k_n{\omega }_i\Delta u^{n,i}-{\displaystyle {\int }_{I_n}{l}_{n,i}\Delta u\text{d}t,}\\ A_0^{n,i}:={\displaystyle \sum _{j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right)\nabla u^{n,j}}-{\displaystyle {\int }_{I_n}{l^{\prime }}_{n,i}\nabla u\text{d}t,}{B}_0^{n,i}:={\displaystyle {\int }_{I_n}{l}_{n,i}\Delta \sigma \text{d}t}-k_n{\omega }_i\Delta {\sigma }^{n,i},\\ J\left[{\eta }_1^n\right]={\eta }_1^n-{\eta }_1^{n+}=w_1^{n+}-w_1^n=\left(P_h^n-P_h^{n-1}\right)\sigma \left(t^n\right)\end{array}$

显然，当$S_h^{n-1}=S_h^n$ 时，$J\left[{\eta }_1^n\right]=0$ ，否则，$J\left[{\eta }_1^n\right]\ne 0$ 。

为获得时间最大模、空间$L^2\left(\Omega \right)$ 模最优误差估计，下面给出误差分析中起重要作用的引理。

引理3.3 [15] 对任意的$n\left(0\le n\le N-1\right)$ 和$i=1,2,\cdots ,q$ ，有如下误差估计式

$\begin{array}{l}\left(\text{a}\right)\Vert {\theta }_1^{n,i}\Vert \le c{k}_n^{\frac{1}{2}}{\left({\displaystyle {\int }_{I_n}{\Vert h_n^{r+1}{\sigma }_t\Vert }_{r+1,h}^2\text{d}t}\right)}^{\frac{1}{2}},\\ \left(\text{b}\right)\Vert A_1^{n,i}\Vert \le c{k}_n^{q+\frac{1}{2}}{\left|\Vert {\sigma }^{\left(q+1\right)}\Vert \right|}_n,\left(\text{c}\right)\Vert B_1^{n,i}\Vert \le c{k}_n^{q+\frac{1}{2}}{\left|\Vert \Delta u^{\left(q\right)}\Vert \right|}_n,\\ \left(\text{d}\right)\Vert A_0^{n,i}\Vert \le c{k}_n^{q+\frac{1}{2}}{\left|\Vert \nabla u^{\left(q+1\right)}\Vert \right|}_n,\left(\text{e}\right)\Vert B_0^{n,i}\Vert \le c{k}_n^{q+\frac{1}{2}}{\left|\Vert \Delta {\sigma }^{\left(q\right)}\Vert \right|}_n.\end{array}$ (19)

该引理的证明过程与文献[15]中证明类似，故这里省略。

引理3.4对于每一个$n\left(0\le n\le N-1\right)$ ，假设$S_h^n\subset S_h^{n-1}$ ，则有估计式

$\begin{array}{l}{\Vert E_1^{n+1}\Vert }^2+{\Vert \nabla E_0^{n+1}\Vert }^2\le \left(1+{\beta }_n\right)\left\{{\Vert E_1^n\Vert }^2+{\Vert \nabla E_0^n\Vert }^2\right\}+c\left\{{\left|\Vert E_1\Vert \right|}_n^2+{\left|\Vert \nabla E_0\Vert \right|}_n^2\right\}+c{k}_n^{2q}\{{\left|\Vert {\sigma }^{\left(q+1\right)}\Vert \right|}_n^2+{\left|\Vert \Delta u^{\left(q\right)}\Vert \right|}_n^2\\ +{\left|\Vert \nabla u^{\left(q+1\right)}\Vert \right|}_n^2+{\left|\Vert \Delta {\sigma }^{\left(q\right)}\Vert \right|}_n^2\}+c{k}_n\underset{I_n}{\max }{\Vert h_n^{r+1}{\sigma }_t\Vert }_{r+1,h}^2+M_n{\Vert J\left[{\eta }_1^n\right]\Vert }^2,\end{array}$

其中$M_n$ 表示与n有关的常数(其具体取值参见后面步骤)。

证 在(18)(a)中取$\psi =E_1^{n,i}$ ，且i从1到q求和，所得方程的左侧表达式和(16)(a)中取$v=E_1$ 时的左侧表达式相等。类似地，在(18)(b)中取$\varphi =E_0^{n,i}$ 得到

$\begin{array}{l}\left(\text{a}\right)\frac{1}{2}{\Vert E_1^{n+1}\Vert }^2+\frac{1}{2}{\Vert E_1^{n+}\Vert }^2+\epsilon {\displaystyle {\int }_{I_n}\left(\nabla E_1,\nabla E_1\right)\text{d}t}+\gamma {\displaystyle {\int }_{I_n}\left(\nabla E_0,\nabla E_1\right)\text{d}t}\\ =\left(E_1^n,E_1^{n+}\right)+{\displaystyle \sum _{i=1}^q\left({\theta }_1^{n,i}+A_1^{n,i}+\gamma B_1^{n,i}-\epsilon B_0^{n,i},E_1^{n,i}\right)}-\left(J\left[{\eta }_1^n\right],E_1^{n+}\right)\\ \left(\text{b}\right)\frac{1}{2}{\Vert \nabla E_0^{n+1}\Vert }^2+\frac{1}{2}{\Vert \nabla E_0^{n+}\Vert }^2-{\displaystyle {\int }_{I_n}\left(\nabla E_1,\nabla E_0\right)\text{d}t}\\ =\left(\nabla E_0^n,\nabla E_0^{n+}\right)+{\displaystyle \sum _{i=1}^q\left(A_0^{n,i},\nabla E_0^{n,i}\right)}+{\displaystyle \sum _{i=1}^q\left(B_0^{n,i},E_0^{n,i}\right)}\end{array}$ (20)

将(a)和(b)式相加，并应用Poincaré不等式和Cauchy-Schwarz不等式，结合引理3.2和引理3.3的结论进行分析，同时注意到${\left|\Vert E_{\kappa }\Vert \right|}_n^2=k_n{\displaystyle \sum _{j=1}^q{\omega }_j{\Vert E_{\kappa }^{n,j}\Vert }^2}$ $\left(\kappa =0,1\right)$ ，可得

$\begin{array}{l}{\Vert E_1^{n+1}\Vert }^2+{\Vert \nabla E_0^{n+1}\Vert }^2\le \left(1+{\beta }_n\right)\left\{{\Vert E_1^n\Vert }^2+{\Vert \nabla E_0^n\Vert }^2\right\}+c\left\{{\left|\Vert E_1\Vert \right|}_n^2+{\left|\Vert \nabla E_0\Vert \right|}_n^2\right\}+c{k}_n^{2q}\{{\left|\Vert {\sigma }^{\left(q+1\right)}\Vert \right|}_n^2+{\left|\Vert \Delta u^{\left(q\right)}\Vert \right|}_n^2\\ +{\left|\Vert \nabla u^{\left(q+1\right)}\Vert \right|}_n^2+{\left|\Vert \Delta {\sigma }^{\left(q\right)}\Vert \right|}_n^2\}+c{k}_n\underset{I_n}{\max }{\Vert h_n^{r+1}{\sigma }_t\Vert }_{r+1,h}^2+M_n{\Vert J\left[{\eta }_1^n\right]\Vert }^2,\end{array}$

其中，$M_n$ 表示与n有关的常数，对于固定的n，取$M_m=M=N_c\left(n\right)\left(m=1,\cdots ,n\right)$ 。后续说明$N_c$ 的定义。

此外，当$S_h^n=S_h^{n-1}$ 时，${\beta }_n=0$ ；当$S_h^n\ne S_h^{n-1}$ 时，${\beta }_n=\frac{1}{M_n-1}$ 。

引理3.5 对任意的$n\left(0\le n\le N-1\right)$ ，假设$S_h^n\subset S_h^{n-1}$ 且$k_n$ 充分小，则有下面的估计式

$\begin{array}{l}{\left|\Vert E_1\Vert \right|}_n^2+{\left|\Vert \nabla E_0\Vert \right|}_n^2\le c{k}_n\left\{{\Vert E_1^n\Vert }^2+{\Vert \nabla E_0^n\Vert }^2+{\Vert J\left[{\eta }_1^n\right]\Vert }^2\right\}+c{k}_n\{k_n\underset{I_n}{\max }{\Vert h_n^{r+1}{\sigma }_t\Vert }_{r+1,h}^2\\ +k_n^{2q+1}\left[{\left|\Vert {\sigma }^{\left(q+1\right)}\Vert \right|}_n^2+{\left|\Vert \Delta u^{\left(q\right)}\Vert \right|}_n^2+{\left|\Vert \nabla u^{\left(q+1\right)}\Vert \right|}_n^2+{\left|\Vert \Delta {\sigma }^{\left(q\right)}\Vert \right|}_n^2\right]\}.\end{array}$

证 式(16)取$E_{\kappa }={\displaystyle \sum _{j=1}^q{\tau }_j^{1/2}{l}_{n,j}\left(t\right){\tilde{E}}_{\kappa }^{n,j}\left(x\right)}$ ，其中${\tilde{E}}_{\kappa }^{n,j}={\tau }_j^{-1/2}{E}_{\kappa }^{n,j}\left(\kappa =0,1\right)$ ，取$v={\tau }_i^{-1/2}{l}_{n,i}\left(t\right)\psi$ 和$w=-{\tau }_i^{-1/2}{l}_{n,i}\left(t\right)\Delta \varphi$ $\left(\psi ,\varphi \in S_h^n\right)$ ，假设$S_h^n\subset S_h^{n-1}$ ，则有

$\begin{array}{l}\left(\text{a}\right){\delta }_{qi}\left({\tilde{E}}_1^{n,q},\psi \right)-{\displaystyle \sum _{j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right){\tau }_j^{1/2}{\tau }_i^{-1/2}}\left({\tilde{E}}_1^{n,j},\psi \right)+\epsilon k_n{\omega }_i\left(\nabla {\tilde{E}}_1^{n,i},\nabla \psi \right)+\gamma k_n{\omega }_i\left(\nabla {\tilde{E}}_0^{n,i},\nabla \psi \right)\\ ={\tau }_i^{-1/2}{l}_{n,i}\left(t^n\right)\left({\tilde{E}}_1^n,\psi \right)-l_{n,i}\left(t^n\right)\left(J\left[{\eta }_1^n\right],\psi \right)+\left({\theta }_1^{n,i}+A_1^{n,i}+\gamma B_1^{n,i}-\epsilon B_0^{n,i},\psi \right)\\ \left(\text{b}\right){\delta }_{qi}\left(\nabla {\tilde{E}}_0^{n,q},\nabla \varphi \right)-{\displaystyle \sum _{j=1}^q{\omega }_{n,j}{l^{\prime }}_{n,i}\left(t^{n,j}\right){\tau }_j^{1/2}{\tau }_i^{-1/2}}\left(\nabla {\tilde{E}}_0^{n,j},\nabla \varphi \right)-k_n{\omega }_i\left(\nabla {\tilde{E}}_1^{n,i},\nabla \varphi \right)\\ ={\tau }_i^{-1/2}{l}_{n,i}\left(t^n\right)\left(\nabla E_0^n,\nabla \varphi \right)+\left(A_0^{n,i},\nabla \varphi \right)+\left(B_0^{n,i},\varphi \right)\end{array}$ (21)

在上式的(a)中取$\psi ={\tilde{E}}_1^{n,i}$ ，(b)中取$\varphi ={\tilde{E}}_0^{n,i}$ ，且i从1到q求和，借助引理3.2和引理3.3的结果，Cauchy-Schwarz不等式和Young不等式，以及${\left|\Vert E_{\kappa }\Vert \right|}_n^2=k_n{\displaystyle \sum _{j=1}^q{\omega }_j{\Vert E_{\kappa }^{n,j}\Vert }^2}\left(\kappa =0,1\right)$ 得到

$\begin{array}{l}{\left|\Vert {\tilde{E}}_1\Vert \right|}^2+{\left|\Vert \nabla {\tilde{E}}_0\Vert \right|}^2\le c\left\{{\Vert {\tilde{E}}_1^n\Vert }^2+{\Vert \nabla {\tilde{E}}_0^n\Vert }^2+{\Vert J\left[{\eta }_1^n\right]\Vert }^2\right\}+c\{k_n\underset{I_n}{\max }{\Vert h_n^{r+1}{\sigma }_t\Vert }_{r+1,h}^2\\ +k_n^{2q+1}\left[{\left|\Vert {\sigma }^{\left(q+1\right)}\Vert \right|}_n^2+{\left|\Vert \Delta u^{\left(q\right)}\Vert \right|}_n^2+{\left|\Vert \nabla u^{\left(q+1\right)}\Vert \right|}_n^2+{\left|\Vert \Delta {\sigma }^{\left(q\right)}\Vert \right|}_n^2\right]\}.\end{array}$

进一步利用${\left|\Vert E_{\kappa }\Vert \right|}_n\le c{k}_n\left|\Vert E_{\kappa }\Vert \right|$ 可证明

$\begin{array}{l}{\left|\Vert E_1\Vert \right|}_n^2+{\left|\Vert \nabla E_0\Vert \right|}_n^2\le c{k}_n\left\{{\Vert E_1^n\Vert }^2+{\Vert \nabla E_0^n\Vert }^2+{\Vert J\left[{\eta }_1^n\right]\Vert }^2\right\}+c{k}_n\{k_n\underset{I_n}{\max }{\Vert h_n^{r+1}{\sigma }_t\Vert }_{r+1,h}^2\\ +k_n^{2q+1}\left[{\left|\Vert {\sigma }^{\left(q+1\right)}\Vert \right|}_n^2+{\left|\Vert \Delta u^{\left(q\right)}\Vert \right|}_n^2+{\left|\Vert \nabla u^{\left(q+1\right)}\Vert \right|}_n^2+{\left|\Vert \Delta {\sigma }^{\left(q\right)}\Vert \right|}_n^2\right]\}.\end{array}$

4. 误差估计

定理4.1 设$\left\{u,\sigma \right\}$ 是方程组(3)的解，$\left\{U,Z\right\}$ 是时间间断TSFE格式(4)的解，并且初始值$U^0=P_h^0{\Phi }_1\left(x\right)$ ，$Z^0=P_h^0{\Phi }_2\left(x\right)$ ，$S_h^n\subset S_h^{n-1}$ ，则有$\sigma =u_t$ 的$L^{\infty }\left(\left[0,T\right];L^2\left(\Omega \right)\right)$ ——模误差估计

$\begin{array}{l}\underset{t\in \left(0,T\right]}{\max }\Vert \sigma -Z\Vert \le c{\text{e}}^{cT}\{\underset{0\le m\le n}{\max }\underset{I_m}{\max }\left({\Vert h_m^{r+1}{\sigma }_t\Vert }_{r+1,h}+{\Vert h_m^{r+1}\sigma \Vert }_{r+1,h}\right)+\underset{0\le m\le n}{\max }{k}_m^q\underset{I_m}{\max }(\Vert {\sigma }^{\left(q+1\right)}\Vert \\ +\Vert \Delta u^{\left(q\right)}\Vert +\Vert \nabla u^{\left(q+1\right)}\Vert +\Vert \Delta {\sigma }^{\left(q\right)}\Vert +\Vert {\sigma }^{\left(q\right)}\Vert )+N_c\underset{1\le m\le n}{\max }\Vert J\left[{\eta }_1^m\right]\Vert \}\end{array}$ (22)

和$u$ 的$L^{\infty }\left(\left[0,T\right];H^1\left(\Omega \right)\right)$ ——模误差估计

$\begin{array}{l}\underset{t\in \left(0,T\right]}{\max }{\Vert u-U\Vert }_1\le c{\text{e}}^{cT}\{\underset{0\le m\le n}{\max }\underset{I_m}{\max }\left({\Vert h_m^{r+1}{\sigma }_t\Vert }_{r+1,h}+{\Vert h_m^{r}u\Vert }_{r+1,h}\right)+\underset{0\le m\le n}{\max }{k}_m^q\underset{I_m}{\max }(\Vert {\sigma }^{\left(q+1\right)}\Vert \\ +\Vert \Delta u^{\left(q\right)}\Vert +\Vert \nabla u^{\left(q+1\right)}\Vert +\Vert \Delta {\sigma }^{\left(q\right)}\Vert +{\Vert u^{\left(q\right)}\Vert }_1)+N_c\underset{1\le m\le n}{\max }\Vert J\left[{\eta }_1^m\right]\Vert \}\end{array}$ (23)

其中$N_c$ 表示$S_h^j\ne S_h^{j-1}\left(j=1,2,\cdots ,N-1\right)$ 的总数，$J\left[w_1^n\right]=w_1^{n+}-w_1^n$ 表示投影$w_1=P_h^n\sigma \left(x,t\right)$ 在点$t^n$ 处的跳跃项$\left(P_h^n-P_h^{n-1}\right)\sigma \left(x,t^n\right)$ 。

证 将引理3.5的结果应用到引理3.4结论的右端，并对n进行迭代，且$E_{\kappa }^0=0\left(\kappa =0,1\right)$ 。整理得

$\begin{array}{l}{\Vert E_1^{n+1}\Vert }^2+{\Vert \nabla E_0^{n+1}\Vert }^2\le c{\displaystyle \sum _{m=0}^n\left({\displaystyle \prod _{j=m+1}^n\left(1+{\beta }_j+c{k}_j\right)}\right)}(\left(c{k}_m+M_m\right){\Vert J\left[{\eta }_1^m\right]\Vert }^2+c{k}_m\underset{I_m}{\max }{\Vert h_m^{r+1}{\sigma }_t\Vert }_{r+1,h}^2\\ +c{k}_m^{2q}\left\{{\left|\Vert {\sigma }^{\left(q+1\right)}\Vert \right|}_m^2+{\left|\Vert \Delta u^{\left(q\right)}\Vert \right|}_m^2+{\left|\Vert \nabla u^{\left(q+1\right)}\Vert \right|}_m^2+{\left|\Vert \Delta {\sigma }^{\left(q\right)}\Vert \right|}_m^2\right\}\end{array}$

其中$M_m$ 是与n有关的常数，$N_c$ 表示$S_h^j\ne S_h^{j-1}\left(j=1,2,\cdots ,N-1\right)$ 的总数，对于固定的n，取$M_m=M=N_c\left(n\right)$ $\left(m=1,\cdots ,n\right)$ ，当$N_c=0,1$ 时，取$M=2$ ，且${\beta }_j$ 满足

${\displaystyle \prod _{j=0}^n\left(1+{\beta }_j+c{k}_j\right)}\le {\displaystyle \prod _{j=0,\beta <c{k}_j}^n\left(1+2c{k}_j\right)}{\displaystyle \prod _{j=0,\beta \ge c{k}_j}^n\left(1+2\beta \right)}\le {\displaystyle \prod _{j=0}^n\left(1+2c{k}_j\right)}{\left(1+2\beta \right)}^M\le 3{\text{e}}^{2c{t}_{n+1}+2}$

现在，置$C_n:={\left(3{\text{e}}^{2c{t}_n+2}\right)}^{1/2}$ 。由于$\Vert J\left[{\eta }_1^0\right]\Vert =0$ 和$k_m\le M$ ，得到

$\begin{array}{l}\Vert E_1^{n+1}\Vert +\Vert \nabla E_0^{n+1}\Vert \le c{C}_{n+1}\{{\displaystyle \sum _{m=0}^n(c{k}_m\underset{I_m}{\max }{\Vert h_m^{r+1}{\sigma }_t\Vert }_{r+1,h}^2+c{k}_m^{2q}({\left|\Vert {\sigma }^{\left(q+1\right)}\Vert \right|}_m^2+{\left|\Vert \Delta u^{\left(q\right)}\Vert \right|}_m^2}\\ +{\left|\Vert \nabla u^{\left(q+1\right)}\Vert \right|}_m^2+{{\left|\Vert \Delta {\sigma }^{\left(q\right)}\Vert \right|}_m^2))\}}^{1/2}+c{C}_{n+1}\sqrt{M}{\left\{{\displaystyle \sum _{m=1}^n{\Vert J\left[{\eta }_1^m\right]\Vert }^2}\right\}}^{1/2}\end{array}$

利用引理3.5和(11)式得

$\begin{array}{l}{\left|\Vert E_1\Vert \right|}_n+{\left|\Vert \nabla E_0\Vert \right|}_n\le c{k}_n^{1/2}{C}_n\{{\displaystyle \sum _{m=0}^n(c{k}_m\underset{I_m}{\max }{\Vert h_m^{r+1}{\sigma }_t\Vert }_{r+1,h}^2+c{k}_m^{2q}({\left|\Vert {\sigma }^{\left(q+1\right)}\Vert \right|}_m^2+{\left|\Vert \Delta u^{\left(q\right)}\Vert \right|}_m^2}\\ +{\left|\Vert \nabla u^{\left(q+1\right)}\Vert \right|}_m^2+{{\left|\Vert \Delta {\sigma }^{\left(q\right)}\Vert \right|}_m^2))\}}^{1/2}+c{k}_n^{1/2}{C}_n\sqrt{M}{\left\{{\displaystyle \sum _{m=1}^n{\Vert J\left[{\eta }_1^m\right]\Vert }^2}\right\}}^{1/2}\end{array}$

因${E_0|}_{I_n}\in V_{hk}^n$ ，${E_1|}_{I_n}\in V_{hk}^n$ ，利用逆不等式

$\underset{I_n}{\max }\Vert y\left(t\right)\Vert \le c_I{k}_n^{-1/2}{\left({\displaystyle {\int }_{I_n}{\Vert y\left(t\right)\Vert }^2}\text{d}t\right)}^{1/2},\forall y\in P_{q-1}\left(I_n\right),c_I>0$

得

$\begin{array}{l}\underset{I_n}{\max }\left[\Vert E_1\Vert +\Vert \nabla E_0\Vert \right]\le c{C}_n\{{\displaystyle \sum _{m=0}^n(c{k}_m\underset{I_m}{\max }{\Vert h_m^{r+1}{\sigma }_t\Vert }_{r+1,h}^2+c{k}_m^{2q}({\left|\Vert {\sigma }^{\left(q+1\right)}\Vert \right|}_m^2+{\left|\Vert \Delta u^{\left(q\right)}\Vert \right|}_m^2}\\ +{\left|\Vert \nabla u^{\left(q+1\right)}\Vert \right|}_m^2+{{\left|\Vert \Delta {\sigma }^{\left(q\right)}\Vert \right|}_m^2))\}}^{1/2}+c{C}_n\sqrt{M}{\left\{{\displaystyle \sum _{m=1}^n{\Vert J\left[{\eta }_1^m\right]\Vert }^2}\right\}}^{1/2}\end{array}$ .

因为$\Vert \sigma -\text{Z}\Vert \le \Vert E_1\Vert +\Vert {\rho }_1\Vert$ ，利用引理3.2得

$\begin{array}{l}\underset{I_n}{\max }\Vert \sigma -Z\Vert \le c{C}_n\{\underset{0\le m\le n}{\max }\underset{I_m}{\max }\left({\Vert h_m^{r+1}{\sigma }_t\Vert }_{r+1,h}+{\Vert h_m^{r+1}\sigma \Vert }_{r+1,h}\right)+\underset{0\le m\le n}{\max }{k}_m^q\underset{I_m}{\max }(\Vert {\sigma }^{\left(q+1\right)}\Vert \\ +\Vert \Delta u^{\left(q\right)}\Vert +\Vert \nabla u^{\left(q+1\right)}\Vert +\Vert \Delta {\sigma }^{\left(q\right)}\Vert +\Vert {\sigma }^{\left(q\right)}\Vert )+N_c\underset{1\le m\le n}{\max }\Vert J\left[{\eta }_1^m\right]\Vert \}\end{array}$ .

故式(22)得证。类理可证(23)，在$H_0^1\left(\Omega \right)$ 中，$\nabla v$ 与${\Vert v\Vert }_1$ 等价，$\Vert \nabla \left(u-U\right)\Vert \le \Vert \nabla E_0\Vert +\Vert \nabla {\rho }_0\Vert$ ，利用引理3.2得

$\begin{array}{l}\underset{t\in \left(0,T\right]}{\max }{\Vert u-U\Vert }_1\le c{C}_n\{\underset{0\le m\le n}{\max }\underset{I_m}{\max }\left({\Vert h_m^{r+1}{\sigma }_t\Vert }_{r+1,h}+{\Vert h_m^{r}u\Vert }_{r+1,h}\right)+\underset{0\le m\le n}{\max }{k}_m^q\underset{I_m}{\max }(\Vert {\sigma }^{\left(q+1\right)}\Vert \\ +\Vert \Delta u^{\left(q\right)}\Vert +\Vert \nabla u^{\left(q+1\right)}\Vert +\Vert \Delta {\sigma }^{\left(q\right)}\Vert +{\Vert u^{\left(q\right)}\Vert }_1)+N_c\underset{1\le m\le n}{\max }\Vert J\left[{\eta }_1^m\right]\Vert \}\end{array}$ .

5. 数值算例

给定如下一维线性粘弹性初边值问题作为数值算例，从而验证该方法误差分析结果的合理性

$\{\begin{array}{l}{u}_{tt}-u_{xxt}-u_{xx}=f\left(x,t\right),x\in \left[0,0.5\right],t\in \left(0,1\right]\\ u\left(0,t\right)=0,u\left(1,t\right)=0,t\in \left[0,1\right]\\ u\left(x,0\right)=\sin \left(2\pi x\right),x\in \left[0,0.5\right]\end{array}$ (24)

其中精确解为$u\left(x,t\right)={\text{e}}^{-t}\sin \left(2\pi x\right)$ ，$f\left(x,t\right)={\text{e}}^{-t}\sin \left(2\pi x\right)$ 。在数值计算过程中，空间上采用一次多项式、时间上采用线性多项式作为基函数。

表1中粘弹性问题(24)中，当时间上取固定步长$k=0.001$ ，而空间上取步长分别为$h=1/8,1/16,1/32$

时，位移项${\Vert u-U\Vert }_{L^{\infty }\left(\left[0,T\right],H^1\left(\Omega \right)\right)}$ 和速度项${\Vert \sigma -Z\Vert }_{L^{\infty }\left(\left[0,T\right],L^2\left(\Omega \right)\right)}$ 的收敛阶和误差估计。由表可知，随着空间步长h不断对半缩减时，${\Vert u-U\Vert }_{L^{\infty }\left(\left[0,T\right],H^1\left(\Omega \right)\right)}$ 的收敛阶接近一阶，${\Vert \sigma -Z\Vert }_{L^{\infty }\left(\left[0,T\right],L^2\left(\Omega \right)\right)}$ 的收敛阶接近二阶，均趋于最

优收敛，与理论推导结果相吻合。

Table 1. The time step $k=0.001$ , the error and convergence order in the spatial direction

表1. 时间步长$k=0.001$ 时，空间方向上误差和收敛阶

h	${\Vert u-U\Vert }_{L^{\infty }\left(\left[0,T\right],H^1\left(\Omega \right)\right)}$	收敛阶	${\Vert \sigma -Z\Vert }_{L^{\infty }\left(\left[0,T\right],L^2\left(\Omega \right)\right)}$	收敛阶
1/8	7.0555e−01		2.7778e−02
1/16	3.5529e−01	0.9897	7.0153e−03	1.9854
1/32	1.7796e−01	0.9974	1.7583e−03	1.9963

表2中取固定空间步长$h=0.0005$ 时，而时间步长分别取$k=1/2,1/4,1/8$ 时，分别给出了位移项${\Vert u-U\Vert }_{L^{\infty }\left(\left[0,T\right],H^1\left(\Omega \right)\right)}$ 和速度项${\Vert \sigma -Z\Vert }_{L^{\infty }\left(\left[0,T\right],L^2\left(\Omega \right)\right)}$ 的收敛阶以及误差估计。由表可知，随着时间步长k不断对半缩减时，${\Vert u-U\Vert }_{L^{\infty }\left(\left[0,T\right],H^1\left(\Omega \right)\right)}$ 和${\Vert \sigma -Z\Vert }_{L^{\infty }\left(\left[0,T\right],L^2\left(\Omega \right)\right)}$ 的收敛阶均趋近二阶最优收敛，与理论分析结果相吻合。数值实验数据表明，本文所提出的TDG-TSFE格式对粘弹性问题求解效果显著，验证了该理论推导的正确性。此外，图1和图2分别是取$h=1/40,k=1/20$ 时，解析解和数值解的数值图形及对比图。

Table 2. The spatial step $h=0.0005$ , the error and convergence order in the time direction

表2. 空间步长$h=0.0005$ 时，时间方向上误差和收敛阶

k	${\Vert u-U\Vert }_{L^{\infty }\left(\left[0,T\right],H^1\left(\Omega \right)\right)}$	收敛阶	${\Vert \sigma -Z\Vert }_{L^{\infty }\left(\left[0,T\right],L^2\left(\Omega \right)\right)}$	收敛阶
1/2	9.6107e−02		1.5454e−02
1/4	2.8251e−02	1.7663	4.5121e−03	1.7761
1/8	8.1389e−03	1.7954	1.2198e−03	1.8871

Figure 1. Comparison chart of the exact solution u and the numerical solution U

图1. 精确解u和数值解U对比图

Figure 2. Comparison chart of the exact solution $\sigma$ and the numerical solution Z

图2. 精确解$\sigma$ 和数值解Z对比图

6. 结论

本文针对粘弹性方程初边值问题构造能同时高精度数值求解位移$u$ 和速度$u_t$ 的一种TDG-TSFE格式。在时间方向引入Radau求积公式和基于Radau积分节点的Lagrange插值基函数，将TDG-TSFE格式转换为不包含时间积分项的有限元格式。基于此，引入时空投影算子和正定性引理，推导出位移$u$ 的$L^{\infty }\left(\left[0,T\right];H^1\left(\Omega \right)\right)$ –模误差估计和速度$\sigma =u_t$ 的$L^{\infty }\left(\left[0,T\right];L^2\left(\Omega \right)\right)$ –模最优误差估计，并通过一个数值例子证明了误差分析结果的合理性。目前本文只考虑了粘弹性方程的线性情形。后期我们将该方法推广应用于非线性粘弹性问题，并利用Brouwer不动点定理、Gronwall引理和Lipschitz条件等进行稳定性分析和误差估计。进一步通过对偶问题的稳定性估计、Galerkin正交性原理，研究时空格式在时间节点上的超收敛估计和后验误差估计，基于此讨论其时空自适应算法。

基金项目

国家自然科学基金(12161034, 12161063)；包头师范学院青年科研项目(BSYKJ2022-ZQ06)。

NOTES

^*通讯作者。