一维扩展的Fisher-Kolmogorov方程的BDF2-直接间断Galerkin算法

The BDF2-Direct Discontinuous Galerkin Algorithm for One-Dimensional Extended Fisher-Kolmogorov Equation

摘要: 本文主要研究一维扩展的Fisher-Kolmogorov方程的有效数值算法。通过结合BDF2时间离散格式与直接间断有限元算法对一维扩展的Fisher-Kolmogorov方程进行求解。首先，引入辅助变量，将四阶的扩展的Fisher-Kolmogorov方程转化为低阶耦合方程，然后利用直接间断有限元求解耦合方程，最后使用BDF2方法，对时间格式进行离散。本文给出了详细的数值算法，并通过一个一维算例进行数值试验，验证了算法的有效性和收敛性。

Abstract: This paper mainly studies the effective numerical algorithm for the one-dimensional extended Fisher-Kolmogorov equation. By combining the BDF2 time discretization format with the direct discontinuous finite element algorithm, the one-dimensional extended Fisher-Kolmogorov equation is solved. Firstly, an auxiliary variable is introduced to transform the fourth-order extended Fisher-Kolmogorov equation into a low-order coupled equation. Then, the coupled equation is solved by using the direct discontinuous finite element method. Finally, the BDF2 method is used to discretize the time scheme. The detailed numerical algorithm is presented in this paper, and a one-dimensional example is used for numerical experiments to verify the effectiveness and convergence of the algorithm.

1. 引言

本文主要对以下具有周期边界条件的一维扩展的Fisher-Kolmogorov (EFK)方程

$u_t+\gamma u_{xxxx}-u_{xx}=-f\left(u\right),\left(x,t\right)\in \overline{\Omega }\times \left(0,T\right],$ (1)

进行了数值研究，其初值为

$u\left(x,0\right)=u_0\left(x\right),x\in \overline{\Omega },$ (2)

其中$u$ 和$u_{xx}$ 都满足周期边界，$f\left(u\right)=u^3-u$ ，$\overline{\Omega }=\left[a,b\right]$ 是空间区间，$\gamma$ ，$T$ 是正常数。

当$\gamma =0$ 时，得到标准的Fisher-Kolmogorov方程，它是由Fisher和Kolmogorov于1937年，在描述生物种群的扩散与适应间的相互作用时提出的。后来，由Coullet等人[1]、Dee和van Saroos [2]-[4]根据实际需要增加了四阶导数项得到扩展的Fisher-Kolmogorov (EFK)方程。目前该方程在相变动力学和模式形成问题中具有重要地位，并且被广泛地应用于双稳态系统的图式形成[2]、液晶中畴壁的传播问题[2]等方面。

对于EFK方程的求解，有以下数值方法。一维情形下，Danumjaya和Pani [5]开发并分析了二阶分裂法与正交三次样条配点法的组合方案，同时提出$C^1$ 连续有限元方法，提升了数值解的精度与连续性，为EFK方程的离散化提供了有效框架。Onyejekwe [6]利用边界元方法对其进行了数值分析。Danumjay [7]使用有限元解的半离散格式和Crank-Nicolson全离散格式，推导出最优误差估计。在二维情形下，Khiari和Omrani [8]、He [9]分别提出了Crank-Nicolson有限差分格式和一种二阶线性隐式有限差分格式，证明了数值解的存在唯一性，并得到了其在$L^{\infty }$ 范数下的收敛性结果。EFK的混合有限元方法也一直受到人们的关注，Danumjaya和Pani [10]利用协调线性三角形$P_1$ 元，通过引入中间变量$v=\Delta u$ 建立了一个低阶混合元格式，在半离散和Euler全离散格式下，证明了逼近解的存在唯一性和稳定性，并得到了在$H^1$ -模意义下和$v$ 在$L^2$ -模意义下$O\left(h^2\right)$ 阶的最优误差估计。Wang等人[11]提出一种新的混合有限元格式求解EFK方程，利用协调线性三角形$P_{{}^1}$ 元和分片常数$P_{{}^0}\times P_{{}^0}$ 元，借助投影算子，在半离散格式和Crank-Nicolson全离散格式下，导出$u$ 和$v$ 在$H^1$ -模意义下的最优误差估计。杨和张[12]通过引入中间变量$v=-\Delta u$ ，在内罚间断有限元框架下用Wilson元，得到原始变量和中间变量的$O\left(h^2\right)$ 和$O\left(h^2+\tau \right)$ 阶的最优误差估计结果。

传统有限差分法(FDM)在处理有间断的问题时存在固有缺陷，其基于网格点值的差分格式在间断处会引发数值振荡。而间断有限元(DG)方法允许单元间分片多项式的不连续性，突破了这一限制。该方法的核心优势在于：通过提升单元内插值多项式的阶数，实现任意阶精度，克服了传统有限元法(FEM)需全局连续性和有限体积法(FVM)受限于一阶精度的不足；在处理复杂边界条件时，DG方法在局部弱形式中处理边界通量，无需网格与边界严格匹配；通过引入数值通量函数，在保持高阶精度的同时有效抑制间断附近的数值振荡。直接间断有限元(DDG)方法相较其他DG方法，不是改变方程的弱格式，而是改变数值流通量。2009年，Liu和Yan [13]基于方程的弱形式，在单元边界处直接引入相容且守恒的数值通量，提出DDG方法求解非线性扩散方程，并且给出保证稳定性的一般形式的数值通量。Liu和Yan [13] [14]于2009年和2010年提出的两种DDG方法都缺乏对称性，因此很难得到$L^2$ 范数下的误差估计。2013年，Vidden和Yan [15]在弱形式中引入测试函数的导数的数值通量，且与$u_x$ 的数值通量相同，得到对称的DDG方法，数值实验表明对称的DDG方法对于数值通量中的系数选择并不敏感，存在大量可容许的数值通量且能够达到最优的$\left(k+1\right)$ 阶精度。

近年来，EFK方程的数值方法研究取得了显著进展。传统有限差分法虽然实现简单，但在处理四阶导数时需引入多步离散，导致计算复杂度增加，且易在间断处引发数值振荡[8] [9]。有限元法通过弱形式离散可部分缓解此问题[5]，但高阶连续性要求限制了其灵活性。局部间断有限元法[16]通过分片多项式提升空间离散的灵活性，但需额外处理数值通量稳定性问题。此外，混合有限元法通过变量降阶简化了四阶项的处理[10] [11]，但非线性项的稳定性仍是挑战。

近期，一些新方法被引入EFK方程，Sun等人[17]通过凸分裂技术与变步长BDF2格式的结合，有效提升了非线性项处理的数值稳定性，为长时间尺度模拟提供了理论保障；王[18]针对四阶算子的间断有限元离散，建立了最优误差估计框架，揭示了间断Galerkin方法在高阶导数问题中的收敛特性；张[19]提出非标准有限元构造方法，为非线性反应项的处理提供了新思路；而Ismail等人[20]则从全离散Galerkin格式的角度，系统分析了时空耦合误差的传播机制。这些工作共同推进了EFK方程数值解法的理论深度。

本文引入中间变量$w=u_{xx}$ 将四阶的EFK方程转化为低阶耦合方程，有效规避传统方法处理高阶导数的连续性难题；在空间上使用对称的DDG方法，克服早期DDG格式非对称性导致的稳定性限制，并实现对间断附近数值振荡的抑制，避免全局连续性约束，同时保证最优$\left(k+1\right)$ 阶空间收敛，时间上采用BDF2方法进行离散，以求解EFK方程。基于全离散DDG数值格式，本文给出了详细的数值算法，并通过一个一维算例进行了数值试验，验证了算法的有效性和收敛性。本文的布局如下：在第2节中，构造了全离散直接间断有限元数值格式；在第3节中，我们提供了一个详细的数值算法；在第4节中，通过数值计算结果，验证算法的有效性；在第5节中，我们对本文做一个简要的总结。

2. DDG全离散数值格式

在本节中，时间离散格式使用BDF2，构造全离散(1)的DDG数值格式。

首先引DDG方法中的常用符号。划分空间区间$\overline{\Omega }=\left[a,b\right]$ ，其中$a=x_{1/2}<x_{2/3}<\cdots <x_{N+1/2}=b$ ，每个单元区间表示为$I_j=\left(x_{j-1/2},x_{j+1/2}\right)$ ，其中$1\le j\le N$ 。单元区间中心和单元区间长度表示为 $x_j=\frac{1}{2}\left(x_{j-1/2}+x_{j+1/2}\right)$ ，$h_j=x_{j-1/2}-x_{j+1/2}$ ，$h={\max }_{1\le j\le N}{h}_j$ ，$u_{j+1/2}^{+}$ 和$u_{j+1/2}^{-}$ 分别表示为$I_{j+1}$ 左端点$x_{j+1/2}$ 处$u$ 的值和$I_j$ 右端点$x_{j+1/2}$ 处$u$ 的值。每个单元边界上$u$ 的跳跃值定义为$\left[u\right]=u^{+}-u^{-}$ ，均值定义为 $\overline{u}=\frac{1}{2}\left(u^{+}+u^{-}\right)$ 。

通过引入辅助变量$w=-u_{xx}$ 将问题(1)转化为以下耦合系统

$\{\begin{array}{l}{u}_t-\gamma w_{xx}+w=-f\left(u\right),\left(x,t\right)\in \overline{\Omega }\times \left(0,T\right],\hfill \\ w=-u_{xx},\left(x,t\right)\in \overline{\Omega }\times \left(0,T\right].\hfill \end{array}$ (3)

其中$f\left(u\right)=u^3-u$ 。然后，为了进一步得到全离散格式，将时间区间$\left[0,T\right]$ 划分为$0=t_0<t_1<\cdots <t_M=T$ 。时间步长为$\tau =T/M$ 且整数$M>0$ 。为方便起见，定义$u^{n+1}=u\left(x,t_{n+1}\right)$ 。

对于任意给定的单元$I_j,j=1,2,\cdots ,N$ ，以及相应的函数$v,r$ ，我们给出内积和$L^2$ 范数的定义如下：

${\left(u,v\right)}_{I_j}={\displaystyle {\int }_{I_j}u}v\text{d}x,{\Vert u\Vert }_{I_j}^2={\left(u,u\right)}_{I_j}.$

对方程(3)两侧同时乘以光滑函数$v$ 和$\phi \in H^1\left(\overline{\Omega }\right)$ ，并在每个单元上进行分部积分可以得到方程(3)的局部变分形式：求$u:\left[0,T\right]\to H^1\left(\overline{\Omega }\right)$ ，使得

$\{\begin{array}{l}{\left(u_t,v\right)}_{I_j}+\gamma {\left(w_x,v_x\right)}_{I_j}-\gamma {w_{x}v|}_{\partial I_j}+{\left(w,v\right)}_{I_j}+{\left(f\left(u\right),v\right)}_{I_j}=0,\hfill \\ {\left(w,\phi \right)}_{I_j}-{\left(u_x,{\phi }_x\right)}_{I_j}+{u_x\phi |}_{\partial I_j}=0.\hfill \end{array}$ (4)

对所有单元求和，可得在$t=t_{n+1}$ 处(3)的全局弱形式为

情形Ⅰ：$n=0$

$\{\begin{array}{l}\left({\delta }_t{u}^1,v\right)+\gamma \left(w_x^1,v_x\right)-\gamma {\displaystyle \sum _{j=1}^N{w_x^1\left[v\right]|}_{\partial I_j}}+\left(w^1,v\right)+\left(f\left(u^1\right),v\right)=\left(R_1^1,v\right),\hfill \\ \left(w^1,\phi \right)-\left(u_x^1,{\phi }_x\right)+{\displaystyle \sum _{j=1}^N{u_x^1\left[\phi \right]|}_{\partial I_j}}=0.\hfill \end{array}$ (5)

情形Ⅱ：$n\ge 1$

$\{\begin{array}{l}\left({\delta }_t{u}^{n+1},v\right)+\gamma \left(w_x^{n+1},v_x\right)-\gamma {\displaystyle \sum _{j=1}^N{w_x^{n+1}\left[v\right]|}_{\partial I_j}}+\left(w^{n+1},v\right)+\left(f\left(u^{n+1}\right),v\right)=\left(R_2^{n+1},v\right),\hfill \\ \left(w^{n+1},\phi \right)-\left(u_x^{n+1},{\phi }_x\right)+{\displaystyle \sum _{j=1}^N{u_x^{n+1}\left[\phi \right]|}_{\partial I_j}}=0.\hfill \end{array}$ (6)

其中时间方向采用BDF2离散形式，表示为：

${\delta }_t{u}^{n+1}=\{\begin{array}{ll}\frac{u^{n+1}-u^n}{\tau },\hfill & n=0,\hfill \\ \frac{3u^{n+1}-4u^n+u^{n-1}}{2\tau },\hfill & n\ge 1.\hfill \end{array}$

并且有

$R_1^1=\frac{u^1-u^0}{\tau }-u_t^1=O\left(\tau \right),R_2^{n+1}=\frac{3u^{n+1}-4u^n+u^{n-1}}{2\tau }-u_t^{n+1}=O\left({\tau }^2\right).$ (7)

为了进一步得到全离散数值格式，定义分片多项式空间

$V_h^k=\left\{v:v\in P^k\left(I_j\right),v\in \left(I_j\right),j=1,\cdots ,N\right\}$

其中$P^k\left(I_j\right)$ 表示每个单元$I_j$ 上次数不超过$k\left(k\ge 0\right)$ 的多项式全体。设$u_h^{n+1},w_h^{n+1}\in V_h^k$ 分别为$u^{n+1},w^{n+1}$ 的近似，为了保证稳定性，增加额外的边界项，对任意试探函数$v_h$ 和${\phi }_h\in V_h^k$ ，可以得到全离散的DDG格式。

情形Ⅰ：$n=0$

$\{\begin{array}{l}\left({\delta }_t{u}_h^1,v_h\right)+\gamma \left(w_{h,x}^1,v_{h,x}\right)+\gamma {\displaystyle \sum _{j=1}^N{\left(\widehat{w_{h,x}^1}\left[v_h\right]\right)}_{j+\frac{1}{2}}}+\gamma {\displaystyle \sum _{j=1}^N{\left(\left[w_h^1\right]\widehat{v_{h,x}}\right)}_{j+\frac{1}{2}}}+\left(w_h^1,v_h\right)+\left(f\left(u_h^1\right),v_h\right)=0,\\ \left(w_h^1,{\phi }_h\right)-\left(u_{h,x}^1,{\phi }_{h,x}\right)-{\displaystyle \sum _{j=1}^N{\left(\widehat{u_{h,x}^1}\left[{\phi }_h\right]\right)}_{j+\frac{1}{2}}}-{\displaystyle \sum _{j=1}^N{\left(\left[u_h^1\right]\widehat{{\phi }_{h,x}}\right)}_{j+\frac{1}{2}}}=0.\end{array}$ (8)

情形Ⅱ：$n\ge 1$

$\{\begin{array}{l}\left({\delta }_t{u}_h^{n+1},v_h\right)+\gamma \left(w_{h,x}^{n+1},v_{h,x}\right)+\gamma {\displaystyle \sum _{j=1}^N{\left(\widehat{w_{h,x}^{n+1}}\left[v_h\right]\right)}_{j+\frac{1}{2}}}+\gamma {\displaystyle \sum _{j=1}^N{\left(\left[w_h^{n+1}\right]\widehat{v_{h,x}}\right)}_{j+\frac{1}{2}}}+\left(w_h^{n+1},v_h\right)+\left(f\left(u_h^{n+1}\right),v_h\right)=0,\\ \left(w_h^{n+1},{\phi }_h\right)-\left(u_{h,x}^{n+1},{\phi }_{h,x}\right)-{\displaystyle \sum _{j=1}^N{\left(\widehat{u_{h,x}^{n+1}}\left[{\phi }_h\right]\right)}_{j+\frac{1}{2}}}-{\displaystyle \sum _{j=1}^N{\left(\left[u_h^{n+1}\right]\widehat{{\phi }_{h,x}}\right)}_{j+\frac{1}{2}}}=0.\end{array}$ (9)

其中，数值流通量选取为：

$\{\begin{array}{l}\widehat{u_{h,x}}=\frac{{\beta }_0}{h}\left[u_h\right]+{\overline{u}}_{h,x}+{\beta }_{1}h\left[u_{h,xx}\right],\hfill \\ \widehat{w_{h,x}}=\frac{{\beta }_0}{h}\left[w_h\right]+{\overline{w}}_{h,x}+{\beta }_{1}h\left[w_{h,xx}\right],\hfill \end{array}$ (10)

$\begin{array}{l}\{\begin{array}{l}{\left({\widehat{v}}_x\right)}_{j+\frac{1}{2}}=-\frac{{\beta }_0}{h}\left(v^{-}\right)+\frac{1}{2}{\left(v_x\right)}^{-}-{\beta }_{1}h{\left(v_{xx}\right)}^{-},\hfill \\ {\left({\widehat{v}}_x\right)}_{j-\frac{1}{2}}=\frac{{\beta }_0}{h}\left(v^{+}\right)+\frac{1}{2}{\left(v_x\right)}^{+}+{\beta }_{1}h{\left(v_{xx}\right)}^{+},\hfill \end{array}\\ \{\begin{array}{l}{\left({\widehat{\phi }}_x\right)}_{j+\frac{1}{2}}=-\frac{{\beta }_0}{h}\left({\phi }^{-}\right)+\frac{1}{2}{\left({\phi }_x\right)}^{-}-{\beta }_{1}h{\left({\phi }_{xx}\right)}^{-},\hfill \\ {\left({\widehat{\phi }}_x\right)}_{j-\frac{1}{2}}=\frac{{\beta }_0}{h}\left({\phi }^{+}\right)+\frac{1}{2}{\left({\phi }_x\right)}^{+}+{\beta }_{1}h{\left({\phi }_{xx}\right)}^{+}.\hfill \end{array}\end{array}$

3. 算法实现

引入如下的双线性格式$\mathbb{B}\left(u,v\right)$ ，

$\mathbb{B}\left(u,v\right)={\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{I_j}{u}_x}}{v}_x\text{d}x+{\displaystyle \sum _{j=1}^N{\left(\widehat{u_x}\left[v\right]\right)}_{j+1/2}}+{\displaystyle \sum _{j=1}^N{\left(\left[u\right]\widehat{v_x}\right)}_{j+1/2}},$

此时可得半离散格式为

$\{\begin{array}{l}\left(u_{h,t},v_h\right)+\gamma \mathbb{B}\left(w_h,v_h\right)+\left(w_h,v_h\right)+\left(f\left(u_h\right),v_h\right)=0,\hfill \\ \left(w_h,{\phi }_h\right)-\mathbb{B}(u_h,{\phi }_h)=0.\hfill \end{array}$ (11)

在第$j$ 个区间上，通过取$u_{h,j},w_{h,j}\in V_h^k$ ，我们有$u_{h,j}={\displaystyle \sum _{i=0}^k{u}_{i,j}}{\varphi }_{i,j},w_{h,j}={\displaystyle \sum _{p=0}^k{w}_{p,j}}{\varphi }_{p,j}$ ，将这些代入到(11)中，并且令$v_h={\varphi }_{m,j},{\phi }_h={\varphi }_{l,j}$ ，我们可以得到(11)的有限元格式的矩阵形式

$\{\begin{array}{l}blockM{U}_t+\gamma blockW+blockMW=F,\hfill \\ blockMW-blockU=0.\hfill \end{array}$ (12)

其中

$\begin{array}{l}blockM={\left[\begin{array}{cccc}{M}_1& & & \\ & M_2& & \\ & & \ddots & \\ & & & M_N\end{array}\right]}_{\left(k+1\right)N\times \left(k+1\right)N},\\ block=blockA+blockC,\\ F=blockMf\left(U\left(t\right)\right),\\ U\left(t\right)={\left(u_1,\cdots ,u_N\right)}^{\text{T}},W\left(t\right)={\left(w_1,\cdots ,w_N\right)}^{\text{T}},\\ u_j=\left(u_{0,j},\cdots ,u_{k,j}\right),w_j=\left(w_{0,j},\cdots ,w_{k,j}\right),j=1,\cdots ,N,\end{array}$

其中，结合边界条件

$blockA={\left[\begin{array}{cccc}{A}_1& & & \\ & A_2& & \\ & & \ddots & \\ & & & A_N\end{array}\right]}_{\left(k+1\right)N\times \left(k+1\right)N},$

$blockC={\left[\begin{array}{cccccc}C{2}_1& C{3}_1& & & & C{1}_1\\ C{1}_2& C{2}_2& C{3}_2& & & \\ & C{1}_3& C{2}_3& C{3}_3& & \\ & & & \ddots & & \\ & & & C{1}_{N-1}& C{2}_{N-1}& C{3}_{N-1}\\ C{3}_N& & & & C{1}_N& C{2}_N\end{array}\right]}_{\left(k+1\right)N\times \left(k+1\right)N}.$

$\begin{array}{l}C{1}_j=2C{B}_j+CB{1}_j+CB{2}_j-EB{1}_j^{\text{T}}-EB{2}_j^{\text{T}};\\ C{2}_j=C21_j+C22_j;\\ C21_j=2D{B}_j+DB{1}_j+DB{2}_j+DB{1}_j^{\text{T}}+DB{2}_j^{\text{T}};\\ C22_j=-2B_j-B{1}_j-B{2}_j-B{1}_j^{\text{T}}-B{2}_j^{\text{T}};\\ C{3}_j=-2E{B}_j-EB{1}_j-EB{2}_j+CB{1}_j^{\text{T}}+CB{2}_j^{\text{T}};\end{array}$

其中

$\begin{array}{l}{M}_j={\left(\left({\varphi }_{i,j},{\varphi }_{m,j}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\\ A_j={\left(\left(\nabla {\varphi }_{i,j},\nabla {\varphi }_{m,j}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\\ B_j=-\frac{{\beta }_0}{h}{\left({\Phi }_{i,j}\left(x_{j+\frac{1}{2}}\right){\Phi }_m\left(x_{j+\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\\ C{B}_j=-\frac{{\beta }_0}{h}{\left({\Phi }_{i,j-1}\left(x_{j-\frac{1}{2}}\right){\Phi }_m\left(x_{j-\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\\ D{B}_j=\frac{{\beta }_0}{h}{\left({\Phi }_{i,j}\left(x_{j-\frac{1}{2}}\right){\Phi }_m\left(x_{j-\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\\ E{B}_j=\frac{{\beta }_0}{h}{\left({\Phi }_{i,j+1}\left(x_{j+\frac{1}{2}}\right){\Phi }_m\left(x_{j+\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\end{array}$

$\begin{array}{l}B{1}_j=\frac{1}{2}{\left({\Phi }_{i,j,x}\left(x_{j+\frac{1}{2}}\right){\Phi }_m\left(x_{j+\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\\ CB{1}_j=\frac{1}{2}{\left({\Phi }_{i,j,x}\left(x_{j-\frac{1}{2}}\right){\Phi }_m\left(x_{j-\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\\ DB{1}_j=\frac{1}{2}{\left({\Phi }_{i,j,x}\left(x_{j-\frac{1}{2}}\right){\Phi }_m\left(x_{j-\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\\ EB{1}_j=\frac{1}{2}{\left({\Phi }_{i,j,x}\left(x_{j+\frac{1}{2}}\right){\Phi }_m\left(x_{j+\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\end{array}$

$\begin{array}{l}B{2}_j=-{\beta }_{1}h{\left({\Phi }_{i,j,xx}\left(x_{j+\frac{1}{2}}\right){\Phi }_m\left(x_{j+\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\\ CB{2}_j=-{\beta }_{1}h{\left({\Phi }_{i,j,xx}\left(x_{j-\frac{1}{2}}\right){\Phi }_m\left(x_{j-\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\\ DB{2}_j={\beta }_{1}h{\left({\Phi }_{i,j,xx}\left(x_{j-\frac{1}{2}}\right){\Phi }_m\left(x_{j-\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)},\\ EB{2}_j=-{\beta }_{1}h{\left({\Phi }_{i,j,xx}\left(x_{j+\frac{1}{2}}\right){\Phi }_m\left(x_{j+\frac{1}{2}}\right)\right)}_{\left(k+1\right)\times \left(k+1\right)}\end{array}$

接下来，我们给出(12)在$n\ge 1$ 的全离散格式，

$\{\begin{array}{l}blockM\frac{3u_h^{n+1}-4u_h^n+u_h^{n-1}}{2\tau }+\gamma block{w}_h^{n+1}+blockM{w}_h^{n+1}=F\left(u_h^n\right),\hfill \\ blockM{w}_h^{n+1}-block{u}_h^{n+1}=0.\hfill \end{array}$ (13)

我们可以将(13)简化为

$\{\begin{array}{l}\left(\gamma block+blockM\right)w_h^{n+1}+\frac{3}{2\tau }blockM{u}_h^{n+1}=F\left(u_h^n\right)+\frac{2}{\tau }blockM{u}_h^n-\frac{1}{2\tau }blockM{u}_h^{n-1},\hfill \\ blockM{w}_h^{n+1}-block{u}_h^{n+1}=0.\hfill \end{array}$ (14)

矩阵形式为

$\left[\begin{array}{cc}\gamma block+blockM& \frac{3}{2\tau }blockM\\ blockM& -block\end{array}\right]\left[\begin{array}{c}{w}_h^{n+1}\\ u_h^{n+1}\end{array}\right]=\left[\begin{array}{c}F\left(u_h^n\right)+\frac{2}{\tau }blockM{u}_h^n-\frac{1}{2\tau }blockM{u}_h^{n-1}\\ 0\end{array}\right]$

4. 数值实验

在本节中，我们将提供一个一维算例来验证算法的有效性和收敛性。

$u_t+\gamma u_{xxxx}-u_{xx}+f\left(u\right)=g\left(x,t\right),\left(x,t\right)\in \overline{\Omega }\times \left(0,T\right],$

初值为

$u\left(x,0\right)=u_0\left(x\right),x\in \overline{\Omega },$

其中$u$ 和$u_{xx}$ 都满足周期边界，非线性项$f\left(u\right)=u^3-u$ ，$\overline{\Omega }=\left[a,b\right]$ 是空间区间，$\gamma$ 是正常数。

在本算例中，我们取$\overline{\Omega }=\left[-4,4\right]$ ，$T=1$ ，初值为$u_0=\sin \left(\pi x\right)$ ，精确解$u\left(x,t\right)={\text{e}}^{-5t}\sin \left(\pi x\right)$ ，源项$g\left(x,t\right)={\text{e}}^{-5t}\sin \left(\pi x\right)\left(-6+\gamma {\pi }^4+{\pi }^2+{\text{e}}^{-10t}{\sin }^2\left(\pi x\right)\right)$ ，表1显示了$\gamma =0.001$ 的计算结果，固定时间步长 $\tau =1/3000$ ，空间步长分别为h = 8/10，8/20，8/40，8/80，8/160，8/320给出算法的误差以及空间收敛阶。表1的数值结果表明全离散格式的空间收敛阶为$O\left(h^{k+1}\right)$ 。

Table 1. Spatial convergence results with a fixed time step $\tau =1/3000$ , $\gamma =0.001$

表1. 固定时间步长$\tau =1/3000$ ，$\gamma =0.001$ 空间收敛结果

$k$	$N$	${\Vert u-u_h\Vert }_{L^2}$	收敛阶
1	10	6.9474E−03	–
	20	3.1494E−03	1.1414
	40	1.0133E−03	1.6360
	80	2.7253E−04	1.8946
	160	6.9451E−05	1.9724
	320	1.7456E−05	1.9922
2	10	1.9750E−03	–
	20	1.8752E−04	3.3967
	40	2.0135E−05	3.2193
	80	2.4098E−06	3.0627
	160	3.1039E−07	2.9568
	320	4.0216E−08	2.9482

Figure 1. $u$ and its numerical solution when $\gamma =0.001,k=2,h=10$

图1. $\gamma =0.001,k=2,h=10$ 时$u$ 和其数值解

图1展示了$\gamma =0.001,k=2,h=8/10$ 时精确解与数值解的对比，可见数值解在$\overline{\Omega }$ 内与精确解基本一致，但在每个小区间上仍存在偏差。图2进一步减小空间步长至$h=8/80$ ，数值解与精确解几乎完全重合，验证了算法的空间收敛性。

Figure 2. $u$ and its numerical solution when $\gamma =0.001,k=2,h=80$

图2. $\gamma =0.001,k=2,h=80$ 时$u$ 和其数值解

5. 总结

本文中，将BDF2时间离散格式与空间上的DDG方法相结合求解EFK方程。通过引入一个辅助变量，将EFK方程转换成一个低阶耦合系统，使用BDF2时间离散格式与直接间断有限元方法得到全离散格式。我们给出了数值算法以及如何进行数值计算。最后通过数值算例验证算法的有效性。

NOTES

^*通讯作者。