The Reduced-Dimension Method for Crank-Nicolson Mixed Finite Element Solution Coefficient Vectors of the Extended Fisher-Kolmogorov Equation

1. Introduction

The primary goal of this research is to apply the POD-based reduced-dimension method to solve the following 2D nonlinear extended Fisher-Kolmogorov (EFK) equation

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{u}_t+\gamma {\mathsf{\Delta }}^{2}u-\mathsf{\Delta }u+f(u)=0,\hfill & (\mathit{x},t)\in \mathsf{\Omega }\times J,\hfill \\ u(\mathit{x},0)=u_0(\mathit{x}),\hfill & \mathit{x}\in \mathsf{\Omega },\hfill \\ u(\mathit{x},t)=\mathsf{\Delta }u(\mathit{x},t)=0,\hfill & (\mathit{x},t)\in \partial \mathsf{\Omega }\times \overline{J},\hfill \end{array}\right.\end{array}$$

(1)

in which $\mathsf{\Omega }\subset R^d(d=2,3)$ is a bounded convex polygonal domain with boundary $\partial \mathsf{\Omega }$. $J=(0,T]$ is the time interval, where $0<T<\infty$. Coefficient $\gamma$ is a positive constant. $u_0(\mathit{x})$ is known initial function. $f(u)=u^3-u$, is the Lipschitz continuous function of u, that is, there exists a positive constant C such that

$$|f(u_1)-f(u_2)|⩽C|u_1-u_2|.$$

(2)

For simplicity, we assume that $u_0(\mathit{x})=0$ in the following theoretical analysis.

When $\gamma =0$ in (1), the second-order diffusion equation is obtained, which is called the standard Fisher-Kolmogorov equation. Coullet et al. [2] added a stabilizing fourth-order derivative term to obtain the equation (1), and called it as the extended Fisher-Kolmogorov equation. The EFK equation has very important physical background. It was applied to various kinds of physics and engineering. Because of the complexity of the equation, it is difficult to get its analytical solution. Therefore, it is very important to solve the equation numerically. And there have been a lot of research results on this equation. For instance, direct local boundary integral equation method [3], Fourier pseudo-spectral method [4], and modified cubic B-spline and polynomial based differential quadrature method [5,6] are applied to solve the EFK equation. Various four-order finite difference schemes [7,8,9,10] have been developed to solve the EFK equations, which have the four-order and second-order accuracy in space and time variables, respectively. He [11] used a three-level linearly implicit finite difference method to solve the 1D and 2D EFK equation, and the method is second-order convergent both in time and space direction. Sweilam et al. [12] studied stochastic EFK equation by combining the compact finite difference scheme and the semi-implicit Euler-Maruyama approach. Danumjaya and Pani [13] applied mixed finite element methods to solve the EFK equation with different types of boundary conditions. Doss and Nandini [14] employed a splitting technique to derive an H${}^1$-Galerkin mixed finite element scheme for the EFK equation. Wang et al. [15] used a new linearized CNMFE method to study the EFK equation, in which $\nabla u$ belongs to the weaker ${(L^2(\mathsf{\Omega }))}^2$ space instead of the classical $H(\mathrm{div};\mathsf{\Omega })$ space. In this paper, we establish the CNMFE scheme for the 2D EFK equation and analysis the uniqueness, stability and convergence of the CNMFE solutions.

However, in practical problems, using CNMFE method to solve EFK equations is a problem with a large number of unknowns, which makes running extremely slow, leads to an increase in rounding errors, and makes it difficult to obtain accurate numerical solutions. Therefore, it is crucial to lessen the number of unknowns in the CNMFE method of EFK equation to reduce the computational load, save CPU running time, and minimize the accumulation of rounding errors in the calculation process. The Proper orthogonal decomposition (POD) technology, reducing the amount of unknowns in numerical schemes, is an efficient dimensionality reduction method for solving partial differential equations (PDEs). It reduces the CPU running time and reduces the accumulation of round-off errors. POD technology can be applied to numerical methods of various PDEs, such as finite difference method [16,17,18], Schwarz domain decomposition method [19], finite spectral element method [20,21] etc. For reduced dimension FE and MFE models based on POD techniques, there are two dimensionality reduction method. One method is that the optimized models are established through reducing the dimension of the FE or MFE subspace. The parabolic equation [22], Sobolev equation [23], Burgers equation [24] and some other PDEs [25,26,27,28] are solved by the reduced-dimension method of constructing the POD subspace. The other is a new reduced-dimension method for the CNFE [29,30] and CNMFE [31,32,33,34] solution coefficient vectors proposed by Luo et al. in 2020. In [35], Luo and Yang combined unknown solution coefficient vector dimensionality reduction method and continuous space-time finite element method to study the 2D non-stationary incompressible Navier-Stokes equations.

Although there are existing reduce-order methods of the EFK equation, such as Literature [17], which combines the high-order compact finite difference scheme with the POD method to obtain the spatial sixth-order accuracy while greatly shortening the calculation time. Literature [33] uses the POD technique to reduce the dimensionality of the two-grid CNMFE solution coefficient vectors, in which, the matrix theory is used to deeply explore the stability and convergence of the TGRDECNMFE solutions. However, the reduce-dimension method proposed in this paper is theoretically significantly different from the above research methods. This paper integrates POD technology with the linearized CNMFE method, and proposes a new reduce-dimension method for 2D nonlinear EFK equations. In this study, we adopt the CN fully discrete format in the temporal direction to achieve second-order time accuracy. Additionally, we linearize the nonlinear term to significantly enhance computational speed and efficiency while preserving model accuracy. Unlike Literature [33], our work utilizes traditional finite element techniques to analyze uniqueness, stability, and convergence of the RDCNMFE solutions. Our proposed method combines the advantages of POD and linearized CNMFE methods for effective dimensionality reduction of the nonlinear EFK equation. This approach not only maintains accuracy but also greatly improves computational speed and efficiency, providing a novel insight and solution for solving the EFK equation.

The structure of this paper is as follows. In Section 2, we propose the CNMFE scheme and discuss the uniqueness, stability and convergence of the CNMFE solutions. In Section 3, we generate the POD basis to derive the reduce-dimension matrix model of the CNMFE solution coefficient vectors and analysis the uniqueness, stability and error estimate of the RDCNMFE solutions. In Section 4, we employ the numerical experiment of the 2D EFK equation to verify the effectiveness of RDCNMFE method and correctness of theoretical analysis. In Section 5, we summarize the major conclusions and findings.

2. The Crank-Nicolson Mixed Finite Element Method of the EFK Equation

2.1. The CNMFE Schemes

The Sobolev spaces and their corresponding norms mentioned in this paper adhere to standard definitions (see [1]). In order to develop the CNMFE scheme of the EFK equation (1), we initially introduce a diffusion term $\omega =-\mathsf{\Delta }u$, resulting in the following two lower-order equations

$$\left\{\begin{array}{c}{u}_t-\gamma \mathsf{\Delta }\omega +\omega +f(u)=0,(\mathit{x},t)\in \mathsf{\Omega }\times J,\hfill \\ -\mathsf{\Delta }u=\omega ,(\mathit{x},t)\in \mathsf{\Omega }\times J.\hfill \end{array}\right.$$

(3)

Let $V=H_0^1(\mathsf{\Omega })$. With the Green’s formula, the mixed functional form of (3) is as follows.

Problem 1. 

Find $\{u,\omega \}$: $[0,T]\mapsto V\times V$ satisfying

$$\begin{array}{c}\hfill \left\{\begin{array}{c}(u_t,v)+\gamma (\nabla \omega ,\nabla v)+(\omega ,v)+(f(u),v)=0,\forall v\in V,\hfill \\ (\omega ,g)-(\nabla u,\nabla g)=0,\forall g\in V,\hfill \\ u(\mathit{x},0)=0,\omega (\mathit{x},0)=0,\mathit{x}\in \mathsf{\Omega }.\hfill \end{array}\right.\end{array}$$

(4)

Let ${\Im }_h$ be the quasi-uniform triangulation subdivision on $\overline{\mathsf{\Omega }}$. The FE subspace $V_h$, spanned by the orthonormal basis ${\left\{{\zeta }_j(x)\right\}}_{j=1}^M$ under the inner product in $H^1(\mathsf{\Omega })$, is defined as follows:

$$\begin{array}{c}{V}_h=\{v_h\in V\cap C(\mathsf{\Omega });v_h{\mid }_K\in P_{M-1}(K),K\in \Im \}=\mathrm{span}\{{\zeta }_{\mathrm{j}}(\mathit{x}):1⩽\mathrm{j}⩽\mathrm{M}\},\hfill \end{array}$$

(5)

in which $P_{M-1}(K)$ is generated by $(M-1)$th degree polynomials on $K\in \Im$.

For given integer $N>0$, let $\mathsf{\Delta }t=T/N$, ${\varphi }^n=\varphi (t_n)$, ${\overline{\partial }}_t{\varphi }^n=({\varphi }^n-{\varphi }^{n-1})/\mathsf{\Delta }t$ and ${\varphi }^{n-{\displaystyle \frac{1}{2}}}=({\varphi }^n+{\varphi }^{n-1})/2$. Therefore, we can rewrite Problem 1 at time $t=t_{n-{\displaystyle \frac{1}{2}}}$ in the following form.

$$\begin{array}{c}\hfill \left\{\begin{array}{c}(u_t(t_{n-{\displaystyle \frac{1}{2}}}),v)+\gamma (\nabla \omega (t_{n-{\displaystyle \frac{1}{2}}}),\nabla v)+(\omega (t_{n-{\displaystyle \frac{1}{2}}}),v)+(f(u(t_{n-{\displaystyle \frac{1}{2}}})),v)=0,\forall v\in V,\hfill \\ (\omega (t_{n-{\displaystyle \frac{1}{2}}}),g)-(\nabla u(t_{n-{\displaystyle \frac{1}{2}}}),\nabla g)=0,\forall g\in V.\hfill \end{array}\right.\end{array}$$

(6)

Then, the equivalent expression of (6) is as follows.

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle (\frac{u^n-u^{n-1}}{\mathsf{\Delta }t},v)+\gamma (\nabla {\omega }^{n-\frac{1}{2}},\nabla v)+({\omega }^{n-\frac{1}{2}},v)+(f({\widehat{\widehat{u}}}^{n-\frac{1}{2}}),v)=(R_1^{n-\frac{1}{2}},v)+(R_2^{n-\frac{1}{2}},v),\forall v\in V,}\hfill \\ ({\omega }^{n-{\displaystyle \frac{1}{2}}},g)-(\nabla u^{n-{\displaystyle \frac{1}{2}}},\nabla g)=0,\forall g\in V,\hfill \end{array}\right.\end{array}$$

(7)

where

$$f({\widehat{\widehat{u}}}^{n-{\displaystyle \frac{1}{2}}})={\displaystyle \frac{3}{2}}f(u^{n-1})-{\displaystyle \frac{1}{2}}f(u^{n-2}),$$

(8)

$$R_1^{n-{\displaystyle \frac{1}{2}}}={\displaystyle \frac{u^n-u^{n-1}}{\mathsf{\Delta }t}}-u_t(t_{n-{\displaystyle \frac{1}{2}}}),$$

(9)

$$R_2^{n-{\displaystyle \frac{1}{2}}}=f({\widehat{\widehat{u}}}^{n-{\displaystyle \frac{1}{2}}})-f(u(t_{n-{\displaystyle \frac{1}{2}}})).$$

(10)

Let $\{u_h^n,{\omega }_h^n\}$ represent the CNMFE approximations of solutions $\{u,\omega \}$ for Problem 1 at $t_n=n\mathsf{\Delta }t$. Thus, the CNMFE model of Problem 1 can be constructed as follows.

Problem 2. 

Find $\{u_h^n,{\omega }_h^n\}\in V_h\times V_h(1⩽n⩽N)$ satisfying

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle ({\overline{\partial }}_t{u}_h^n,v_h)+\gamma (\nabla {\omega }_h^{n-\frac{1}{2}},\nabla v_h)+({\omega }_h^{n-\frac{1}{2}},v_h)+(f({\widehat{\widehat{u}}}_h^{n-\frac{1}{2}}),v_h)=0,\forall v_h\in V_h,}\hfill \\ ({\omega }_h^{n-{\displaystyle \frac{1}{2}}},g_h)-(\nabla u_h^{n-{\displaystyle \frac{1}{2}}},\nabla g_h)=0,\forall g_h\in V_h,\hfill \\ u_h^0(\mathit{x})=0,{\omega }_h^0(\mathit{x})=0,\mathit{x}\in \mathsf{\Omega }.\hfill \end{array}\right.\end{array}$$

(11)

Remark 1. 

By means of a linearized term

$$f({\widehat{\widehat{u}}}_h^{n-{\displaystyle \frac{1}{2}}})={\displaystyle \frac{3}{2}}f(u_h^{n-1})-{\displaystyle \frac{1}{2}}f(u_h^{n-2})={\displaystyle \frac{3}{2}}{(u_h^{n-1})}^3-{\displaystyle \frac{1}{2}}{(u_h^{n-2})}^3-{\displaystyle \frac{3}{2}}{u}_h^{n-1}+{\displaystyle \frac{1}{2}}{u}_h^{n-2},$$

(12)

we can observe that the formulation (11) is a linear scheme.

2.2. The Uniqueness, Stability and Convergence of the CNMFE Solutions

By employing the base functions of the FE space $V_h$, the solutions $\{u_h^n,{\omega }_h^n\}$ of Problem 2 can be formulated in the form of

$$\begin{array}{c}{\displaystyle u_h^n=\sum _{j=1}^M{U}_{hj}^n{\zeta }_j=\zeta ·{\mathit{U}}_h^n,{\omega }_h^n=\sum _{j=1}^M{W}_{hj}^n{\zeta }_j=\zeta ·{\mathit{W}}_h^n,}\hfill \end{array}$$

(13)

where $\zeta =({\zeta }_1,{\zeta }_2,\cdots ,{\zeta }_M)$ is the orthonormal base, satisfied

$$({\zeta }_i,{\zeta }_j)={\delta }_{ij},i,j=1,2,\cdots ,M.$$

${\mathit{U}}_h^n={(U_{h1}^n,U_{h2}^n,\cdots ,U_{hM}^n)}^T$ and ${\mathit{W}}_h^n={(W_{h1}^n,W_{h2}^n,\cdots ,W_{hM}^n)}^T$ are the unknown coefficient vectors of CNMFE solutions. From the expression of the solutions $\{u_h^n,{\omega }_h^n\}$ of (13), Problem 2 can be represented equivalently in the following matrix model.

Problem 3. 

Find $\{{\mathit{U}}_h^n,{\mathit{W}}_h^n\}\in R^M\times R^M$ and $\{u_h^n,{\omega }_h^n\}\in V_h\times V_h(1⩽n⩽N)$ satisfying

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle {\overline{\partial }}_t{\mathit{U}}_h^n+\gamma \mathit{B}{\mathit{W}}_h^{n-\frac{1}{2}}+{\mathit{W}}_h^{n-\frac{1}{2}}=\mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-\frac{1}{2}}),}\hfill \\ {\mathit{W}}_h^{n-{\displaystyle \frac{1}{2}}}-{\mathit{BU}}_h^{n-{\displaystyle \frac{1}{2}}}=0,\hfill \\ {\displaystyle u_h^n=\sum _{j=1}^M{U}_{hj}^n{\zeta }_j={\mathit{U}}_h^n·\zeta ,{\omega }_h^n=\sum _{j=1}^M{W}_{hj}^n{\zeta }_j={\mathit{W}}_h^n·\zeta ,}\hfill \end{array}\right.\end{array}$$

(14)

where

$$\mathit{B}={(b_{ij})}_{1⩽i,j⩽M},b_{ij}=(\nabla {\zeta }_j,\nabla {\zeta }_i),$$

and

$$\begin{array}{cc}{\displaystyle \mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-\frac{1}{2}})}\hfill & {\displaystyle =-(\frac{3}{2}\mathit{F}({\mathit{U}}_h^{n-1})-\frac{1}{2}\mathit{F}({\mathit{U}}_h^{n-2}))}\hfill \\ \hfill & {\displaystyle ={(-(\frac{3}{2}(f(\sum _{j=1}^M{U}_{hj}^{n-1}{\zeta }_j),{\zeta }_i)-\frac{1}{2}(f(\sum _{j=1}^M{U}_{hj}^{n-2}{\zeta }_j),{\zeta }_i)))}_{1⩽i⩽M.}}\hfill \end{array}$$

Theorem 1. 

The Problem 3 has the unique solutions $\{u_h^n,{\omega }_h^n\}\in V_h\times V_h$.

Proof of Theorem 1. 

Problem 3 can be written as

$$\left(\begin{array}{cc}\mathit{I}& {\displaystyle \frac{\gamma \mathsf{\Delta }t}{2}}\mathit{B}+{\displaystyle \frac{\mathsf{\Delta }t}{2}}\mathit{I}\\ -\mathit{B}& \mathit{I}\end{array}\right)\left(\begin{array}{c}{\mathit{U}}_h^n\\ {\mathit{W}}_h^n\end{array}\right)=\left(\begin{array}{cc}\mathit{I}& -{\displaystyle \frac{\gamma \mathsf{\Delta }t}{2}}\mathit{B}-{\displaystyle \frac{\mathsf{\Delta }t}{2}}\mathit{I}\\ \mathit{B}& -\mathit{I}\end{array}\right)\left(\begin{array}{c}{\mathit{U}}_h^{n-1}\\ {\mathit{W}}_h^{n-1}\end{array}\right)-\mathsf{\Delta }t\left(\begin{array}{c}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-{\displaystyle \frac{1}{2}}})\\ \tilde{\mathit{O}}\end{array}\right),$$

(15)

where $\mathit{I}$ stands for the $M\times M$ unit matrix, and $\tilde{\mathit{O}}$ denotes zero column vector of $M\times 1$.

From the definitions of matrices $\mathit{I}$ and $\mathit{B}$, we know that they are positive definite and thus invertible. From (15), we can conclude that the coefficient matrix is invertible as well. Hence, there exists a unique set of solutions $(u_h^n,{\omega }_h^n)\in V_h\times V_h$ for Problem 3.

□

In order to prove the stability of the CNMFE solutions, we need to introduce the properties of $\mathit{B}$ in Problem 3.

Lemma 1. 

[[36], Lemma 1.19] Matrix $\mathit{B}$ satisfies the following inequality

$${\parallel \mathit{B}\parallel }_{\infty }⩽C,{\parallel {\mathit{B}}^{-1}\parallel }_{\infty }⩽C.$$

(16)

Theorem 2. 

Both the CNMFE solutions $\{u_h^n,{\omega }_h^n\}$ and their solution coefficient vectors $\{{\mathit{U}}_h^n,{\mathit{W}}_h^n\}$ are unconditionally stable.

Proof of Theorem 2. 

We represent (14) as

$$\left\{\begin{array}{c}{\displaystyle {\mathit{B}}^{-1}{\overline{\partial }}_t{\mathit{U}}_h^n+\gamma {\mathit{W}}_h^{n-\frac{1}{2}}+{\mathit{B}}^{-1}{\mathit{W}}_h^{n-\frac{1}{2}}={\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-\frac{1}{2}}),}\hfill \\ {\mathit{W}}_h^{n-{\displaystyle \frac{1}{2}}}={\mathit{BU}}_h^{n-{\displaystyle \frac{1}{2}}}.\hfill \end{array}\right.$$

(17)

Substituting the second formula of (17) into the first formula, we obtain

$${\displaystyle {\mathit{B}}^{-1}{\overline{\partial }}_t{\mathit{U}}_h^n+\gamma {\mathit{BU}}_h^{n-\frac{1}{2}}+{\mathit{U}}_h^{n-\frac{1}{2}}={\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-\frac{1}{2}}).}$$

(18)

Taking the inner product of (18) and ${\overline{\partial }}_t{\mathit{U}}_h^n$, we have

$${\displaystyle ({\mathit{B}}^{-1}{\overline{\partial }}_t{\mathit{U}}_h^n,{\overline{\partial }}_t{\mathit{U}}_h^n)+\gamma ({\mathit{BU}}_h^{n-\frac{1}{2}},{\overline{\partial }}_t{\mathit{U}}_h^n)+({\mathit{U}}_h^{n-\frac{1}{2}},{\overline{\partial }}_t{\mathit{U}}_h^n)=({\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-\frac{1}{2}}),{\overline{\partial }}_t{\mathit{U}}_h^n).}$$

(19)

The left side of (19) is that

$$\begin{array}{cc}\hfill & {\displaystyle ({\mathit{B}}^{-1}{\overline{\partial }}_t{\mathit{U}}_h^n,{\overline{\partial }}_t{\mathit{U}}_h^n)+\gamma ({\mathit{BU}}_h^{n-\frac{1}{2}},{\overline{\partial }}_t{\mathit{U}}_h^n)+({\mathit{U}}_h^{n-\frac{1}{2}},{\overline{\partial }}_t{\mathit{U}}_h^n)}\hfill \\ =\hfill & {\displaystyle \parallel {({\mathit{B}}^{-1})}^{\frac{1}{2}}{\overline{\partial }}_t{\mathit{U}}_h^n{\parallel }^2+\frac{\gamma }{2\mathsf{\Delta }t}(\parallel {\mathit{B}}^{\frac{1}{2}}{\mathit{U}}_h^n{\parallel }^2-\parallel {\mathit{B}}^{\frac{1}{2}}{\mathit{U}}_h^{n-1}{\parallel }^2)+\frac{1}{2\mathsf{\Delta }t}(\parallel {\mathit{U}}_h^n{\parallel }^2-\parallel {\mathit{U}}_h^{n-1}{\parallel }^2),}\hfill \end{array}$$

(20)

and the right side of (19) is that

$${\displaystyle ({\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-\frac{1}{2}}),{\overline{\partial }}_t{\mathit{U}}_h^n)⩽C\parallel {({\mathit{B}}^{-1})}^{\frac{1}{2}}{\parallel }_{\infty }^2\parallel \mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-\frac{1}{2}})\parallel +\parallel {({\mathit{B}}^{-1})}^{\frac{1}{2}}{\overline{\partial }}_t{\mathit{U}}_h^n{\parallel }^2.}$$

(21)

Combining lemma 1 with the global Lipschitz-continuity of $f(u)$, the first term of (21) can be estimated as

$$\begin{array}{cc}\hfill & {\displaystyle \parallel {({\mathit{B}}^{-1})}^{\frac{1}{2}}{\parallel }_{\infty }^2{\parallel \mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-\frac{1}{2}})\parallel }^2}\hfill \\ ⩽\hfill & {\displaystyle C\parallel \frac{3}{2}\mathit{F}({\mathit{U}}_h^{n-1})-\frac{1}{2}\mathit{F}({\mathit{U}}_h^{n-2}){\parallel }^2}\hfill \\ ⩽\hfill & {\displaystyle C\left(\parallel \frac{1}{2}\mathit{F}({\mathit{U}}_h^{n-1})-\frac{1}{2}\mathit{F}({\mathit{U}}_h^{n-2}){\parallel }^2+\parallel \mathit{F}({\mathit{U}}_h^{n-1})-\mathit{F}({\mathit{U}}_h^0){\parallel }^2+{\parallel \mathit{F}({\mathit{U}}_h^0)\parallel }^2\right)}\hfill \\ ⩽\hfill & {\displaystyle C(\parallel {\mathit{U}}_h^{n-1}-{\mathit{U}}_h^{n-2}{\parallel }^2+\parallel {\mathit{U}}_h^{n-1}-{\mathit{U}}_h^0{\parallel }^2+\parallel \mathit{F}({\mathit{U}}_h^0){\parallel }^2)}\hfill \\ ⩽\hfill & {\displaystyle C(\parallel {\mathit{U}}_h^0{\parallel }^2+\parallel \mathit{F}({\mathit{U}}_h^0){\parallel }^2+\parallel {\mathit{U}}_h^{n-1}{\parallel }^2+\parallel {\mathit{U}}_h^{n-2}{\parallel }^2)}\hfill \\ ⩽\hfill & {\displaystyle C+C(\parallel {\mathit{U}}_h^{n-1}{\parallel }^2+\parallel {\mathit{U}}_h^{n-2}{\parallel }^2),2⩽n⩽N,}\hfill \end{array}$$

(22)

where $\parallel {\mathit{U}}_h^0\parallel$ and $\parallel \mathit{F}({\mathit{U}}_h^0)\parallel$ are bounded. Combining (20), (21) and (22) we have

$$\begin{array}{c}{\displaystyle \frac{\gamma }{2\mathsf{\Delta }t}(\parallel {\mathit{B}}^{\frac{1}{2}}{\mathit{U}}_h^n{\parallel }^2-\parallel {\mathit{B}}^{\frac{1}{2}}{\mathit{U}}_h^{n-1}{\parallel }^2)+\frac{1}{2\mathsf{\Delta }t}(\parallel {\mathit{U}}_h^n{\parallel }^2-\parallel {\mathit{U}}_h^{n-1}{\parallel }^2)⩽C+C(\parallel {\mathit{U}}_h^{n-1}{\parallel }^2+\parallel {\mathit{U}}_h^{n-2}{\parallel }^2).}\hfill \end{array}$$

(23)

Multiplying both sides of the inequality of (23) by $2\mathsf{\Delta }t$ and summating from 2 to $n(n⩽N)$, we have

$$\begin{array}{cc}\parallel {\mathit{U}}_h^n{\parallel }^2\hfill & {\displaystyle ⩽\parallel {\mathit{U}}_h^1{\parallel }^2+\gamma \parallel {\mathit{B}}^{\frac{1}{2}}{\mathit{U}}_h^1{\parallel }^2+2\mathsf{\Delta }t\sum _{i=2}^{n}C+C\mathsf{\Delta }t\sum _{i=0}^{n-1}{\parallel {\mathit{U}}_h^i\parallel }^2}\hfill \\ \hfill & {\displaystyle ⩽(1+\gamma \parallel {\mathit{B}}^{\frac{1}{2}}{\parallel }_{\infty })\parallel {\mathit{U}}_h^1{\parallel }^2+CT+C\mathsf{\Delta }t\sum _{i=0}^{n-1}{\parallel {\mathit{U}}_h^i\parallel }^2,}\hfill \end{array}$$

(24)

noting that $\parallel {\mathit{U}}_h^1\parallel ⩽C$, then

$${\displaystyle \parallel {\mathit{U}}_h^n{\parallel }^2⩽C+C\mathsf{\Delta }t\sum _{i=0}^{n-1}{\parallel {\mathit{U}}_h^i\parallel }^2.}$$

(25)

Applying the Gronwall Lemma to (25), we obtain

$$\parallel {\mathit{U}}_h^n{\parallel }^2⩽C{e}^{Cn\mathsf{\Delta }t}⩽C,1⩽n⩽N.$$

(26)

And because

$$\parallel {\mathit{W}}_h^n\parallel =\parallel \mathit{B}{\mathit{U}}_h^n{\parallel ⩽\parallel \mathit{B}\parallel }_{\infty }\parallel {\mathit{U}}_h^n\parallel ⩽C\parallel {\mathit{U}}_h^n\parallel ⩽C,1⩽n⩽N.$$

(27)

From (26) and (27), we have

$$\parallel {\mathit{U}}_h^n\parallel +\parallel {\mathit{W}}_h^n\parallel ⩽C,1⩽n⩽N.$$

(28)

Due to $\parallel \zeta \parallel ⩽C$, we can easily get

$$\parallel u_h^n\parallel +\parallel {\omega }_h^n\parallel =\parallel {\mathit{U}}_h^n·\zeta \parallel +\parallel {\mathit{W}}_h^n·\zeta \parallel ⩽C\parallel {\mathit{U}}_h^n\parallel \parallel \zeta \parallel +C\parallel {\mathit{W}}_h^n\parallel \parallel \zeta \parallel ⩽C,1⩽n⩽N.$$

(29)

(28) and (29) indicate that the CNMFE solutions and their coefficient vectors are unconditionally stable.

□

In order to analyze the convergence of the linearized CNMFE scheme, we first need to establish the projection operator $R_h$.

Lemma 2. 

Let $R_h:V\to V_h$ satisfy the following relation

$$(\nabla (\varphi -R_h\varphi ),\nabla \chi )=0,\forall \chi \in V_h,$$

(30)

and

$$\parallel \varphi -R_h{\varphi \parallel }_{L^2(\mathsf{\Omega })}+h\parallel \varphi -R_h{\varphi \parallel }_{H^1(\mathsf{\Omega })}⩽C{h}^2{\parallel \varphi \parallel }_{H^2(\mathsf{\Omega })},$$

(31)

$$\parallel {\varphi }_t-R_h{\varphi }_t{\parallel }_{L^2(\mathsf{\Omega })}⩽C{h}^2{(\parallel \varphi \parallel }_{H^2(\mathsf{\Omega })}+\parallel {\varphi }_t{\parallel }_{H^2(\mathsf{\Omega })}).$$

(32)

In order to facilitate theoretical analysis, the errors can be decomposed into

$$u(t_n)-u_h^n=u(t_n)-R_{h}u(t_n)+R_{h}u(t_n)-u_h^n={\tau }^n+{\theta }^n,$$

(33)

$$\omega (t_n)-{\omega }_h^n=\omega (t_n)-R_h\omega (t_n)+R_h\omega (t_n)-{\omega }_h^n={\pi }^n+{\kappa }^n.$$

(34)

Subtracting (11) from (7) and using (30) at $t=t_{n-{\displaystyle \frac{1}{2}}}$, we obtain the following error equations.

$$\begin{array}{c}{\displaystyle (\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t},v_h)+\gamma (\nabla {\kappa }^{n-\frac{1}{2}},\nabla v_h)+({\kappa }^{n-\frac{1}{2}},v_h)+(f({\widehat{\widehat{u}}}^{n-\frac{1}{2}})-f({\widehat{\widehat{u}}}_h^{n-\frac{1}{2}}),v_h)}\hfill \\ {\displaystyle =-(\frac{{\tau }^n-{\tau }^{n-1}}{\mathsf{\Delta }t},v_h)-({\pi }^{n-\frac{1}{2}},v_h)+(R_1^{n-\frac{1}{2}},v_h)+(R_2^{n-\frac{1}{2}},v_h),\forall v_h\in V_h,}\hfill \end{array}$$

(35)

$$({\kappa }^{n-{\displaystyle \frac{1}{2}}},g_h)-(\nabla {\theta }^{n-{\displaystyle \frac{1}{2}}},\nabla g_h)=-({\pi }^{n-{\displaystyle \frac{1}{2}}},g_h),\forall g_h\in V_h.$$

(36)

In order to obtain the error estimates, we introduce the following lemma, which is easily obtained by Taylor expansion.

Lemma 3. 

$R_1^{n-{\displaystyle \frac{1}{2}}}$ and $R_2^{n-{\displaystyle \frac{1}{2}}}$ satisfy the following estimates

$$\parallel R_1^{n-{\displaystyle \frac{1}{2}}}{\parallel }^2⩽C\mathsf{\Delta }{t}^3{\parallel u_{ttt}\parallel }_{L^{\infty }(L^2)},$$

(37)

$$\parallel R_2^{n-{\displaystyle \frac{1}{2}}}{\parallel }^2⩽C(u)\mathsf{\Delta }{t}^3{\parallel u_{ttt}\parallel }_{L^{\infty }(L^2)}.$$

(38)

Based on Lemmas 2 and 3, the following theorems for Crank-Nicolson fully discrete error estimates is obtained.

Theorem 3. 

With $u_h^1=R_{h}u(t_1)$, assume that the solution’s regular properties for (4) satisfy $u_t\in L^2(H^2)$, $u_{ttt}\in L^{\infty }(L^2)$, $u,\omega \in L^{\infty }(H^2)$. Then there exists a positive constant C independent of h and $\mathsf{\Delta }t$ such that

$$\begin{array}{c}\parallel u(t_J)-u_h^J\parallel ⩽C(h^2\parallel |*|\parallel +{(\mathsf{\Delta }t)}^2\parallel u_{ttt}{\parallel }_{L^{\infty }(L^2)}),\hfill \end{array}$$

(39)

where $\parallel |*|\parallel \triangleq {\parallel u\parallel }_{L^{\infty }(H^2)}+{\parallel \omega \parallel }_{L^{\infty }(H^2)}+{\parallel u_t\parallel }_{L^2(H^2)}$.

Proof of Theorem 3. 

Taking $v_h={\theta }^{n-{\displaystyle \frac{1}{2}}}$ and $g_h={\kappa }^{n-{\displaystyle \frac{1}{2}}}$ in (35) and (36), respectively, we obtain

$$\begin{array}{c}{\displaystyle \frac{\parallel {\theta }^n{\parallel }^2-{\parallel {\theta }^{n-1}\parallel }^2}{2\mathsf{\Delta }t}+\gamma (\nabla {\kappa }^{n-\frac{1}{2}},\nabla {\theta }^{n-\frac{1}{2}})+({\kappa }^{n-\frac{1}{2}},{\theta }^{n-\frac{1}{2}})+(f({\widehat{\widehat{u}}}^{n-\frac{1}{2}})-f({\widehat{\widehat{u}}}_h^{n-\frac{1}{2}}),{\theta }^{n-\frac{1}{2}})}\hfill \\ {\displaystyle =-(\frac{{\tau }^n-{\tau }^{n-1}}{\mathsf{\Delta }t},{\theta }^{n-\frac{1}{2}})-({\pi }^{n-\frac{1}{2}},{\theta }^{n-\frac{1}{2}})+(R_1^{n-\frac{1}{2}},{\theta }^{n-\frac{1}{2}})+(R_2^{n-\frac{1}{2}},{\theta }^{n-\frac{1}{2}}),}\hfill \end{array}$$

(40)

$$\parallel {\kappa }^{n-{\displaystyle \frac{1}{2}}}{\parallel }^2-(\nabla {\theta }^{n-{\displaystyle \frac{1}{2}}},\nabla {\kappa }^{n-{\displaystyle \frac{1}{2}}})=-({\pi }^{n-{\displaystyle \frac{1}{2}}},{\kappa }^{n-{\displaystyle \frac{1}{2}}}).$$

(41)

Multiplying (41) by $\gamma$ and adding it to (40) and using Cauchy-Schwarz inequality and Young inequality, we obtain

$$\begin{array}{cc}\hfill & {\displaystyle \gamma \parallel {\kappa }^{n-\frac{1}{2}}{\parallel }^2+\frac{1}{2\mathsf{\Delta }t}(\parallel {\theta }^n{\parallel }^2-\parallel {\theta }^{n-1}{\parallel }^2)}\hfill \\ =\hfill & {\displaystyle -({\kappa }^{n-\frac{1}{2}},{\theta }^{n-\frac{1}{2}})-(f({\widehat{\widehat{u}}}^{n-\frac{1}{2}})-f({\widehat{\widehat{u}}}_h^{n-\frac{1}{2}}),{\theta }^{n-\frac{1}{2}})-\left(\frac{{\tau }^n-{\tau }^{n-1}}{\mathsf{\Delta }t},{\theta }^{n-\frac{1}{2}}\right)-({\pi }^{n-\frac{1}{2}},{\theta }^{n-\frac{1}{2}})}\hfill \\ \hfill & +(R_1^{n-{\displaystyle \frac{1}{2}}},{\theta }^{n-{\displaystyle \frac{1}{2}}})+(R_2^{n-{\displaystyle \frac{1}{2}}},{\theta }^{n-{\displaystyle \frac{1}{2}}})-\gamma ({\pi }^{n-{\displaystyle \frac{1}{2}}},{\kappa }^{n-{\displaystyle \frac{1}{2}}})\hfill \\ ⩽\hfill & {\displaystyle \left(\parallel \frac{{\tau }^n-{\tau }^{n-1}}{\mathsf{\Delta }t}\parallel +\parallel R_1^{n-\frac{1}{2}}\parallel +\parallel R_2^{n-\frac{1}{2}}\parallel \right)\parallel {\theta }^{n-\frac{1}{2}}\parallel +\parallel {\kappa }^{n-\frac{1}{2}}\parallel \parallel {\theta }^{n-\frac{1}{2}}\parallel +\gamma \parallel {\pi }^{n-\frac{1}{2}}\parallel \parallel {\kappa }^{n-\frac{1}{2}}\parallel }\hfill \\ \hfill & +(\parallel f({\widehat{\widehat{u}}}^{n-{\displaystyle \frac{1}{2}}})-f({\widehat{\widehat{u}}}_h^{n-{\displaystyle \frac{1}{2}}})\parallel +\parallel {\pi }^{n-{\displaystyle \frac{1}{2}}}\parallel )\parallel {\theta }^{n-{\displaystyle \frac{1}{2}}}\parallel \hfill \\ ⩽\hfill & {\displaystyle C\left(\parallel \frac{{\tau }^n-{\tau }^{n-1}}{\mathsf{\Delta }t}{\parallel }^2+\parallel R_1^{n-\frac{1}{2}}{\parallel }^2+\parallel R_2^{n-\frac{1}{2}}{\parallel }^2+{\parallel {\pi }^{n-\frac{1}{2}}\parallel }^2\right)+\frac{\gamma }{2}{\parallel {\kappa }^{n-\frac{1}{2}}\parallel }^2}\hfill \\ \hfill & +C\parallel f({\widehat{\widehat{u}}}^{n-{\displaystyle \frac{1}{2}}})-f({\widehat{\widehat{u}}}_h^{n-{\displaystyle \frac{1}{2}}}){\parallel }^2+C{\parallel {\theta }^{n-{\displaystyle \frac{1}{2}}}\parallel }^2.\hfill \end{array}$$

(42)

For the nonlinear term $\parallel f({\widehat{\widehat{u}}}^{n-{\displaystyle \frac{1}{2}}})-f({\widehat{\widehat{u}}}_h^{n-{\displaystyle \frac{1}{2}}}){\parallel }^2$, Using Young inequality, we have

$$\begin{array}{cc}\hfill & \parallel f({\widehat{\widehat{u}}}^{n-{\displaystyle \frac{1}{2}}})-f({\widehat{\widehat{u}}}_h^{n-{\displaystyle \frac{1}{2}}}){\parallel }^2\hfill \\ =\hfill & {\displaystyle \parallel \frac{3}{2}f(u^{n-1})-\frac{1}{2}f(u^{n-2})-\frac{3}{2}f(u_h^{n-1})+\frac{1}{2}f(u_h^{n-2}){\parallel }^2}\hfill \\ =\hfill & {\displaystyle \parallel \frac{3}{2}{(u^{n-1})}^3-\frac{3}{2}{u}^{n-1}-\frac{1}{2}{(u^{n-2})}^3+\frac{1}{2}{u}^{n-2}-\frac{3}{2}{(u_h^{n-1})}^3+\frac{3}{2}{u}_h^{n-1}+\frac{1}{2}{(u_h^{n-2})}^3-\frac{1}{2}{u}_h^{n-2}{\parallel }^2}\hfill \\ =\hfill & {\displaystyle \parallel \frac{3}{2}[{(u^{n-1})}^3-{(u_h^{n-1})}^3]-\frac{1}{2}[{(u^{n-2})}^3-{(u_h^{n-2})}^3]-\frac{3}{2}(u^{n-1}-u_h^{n-1})+\frac{1}{2}(u^{n-2}-u_h^{n-2}){\parallel }^2}\hfill \\ ⩽\hfill & {\displaystyle \parallel \frac{3}{2}(u^{n-1}-u_h^{n-1})[{(u^{n-1})}^2+u^{n-1}{u}_h^{n-1}+{(u_h^{n-1})}^2]{\parallel }^2}\hfill \\ \hfill & {\displaystyle +\parallel \frac{1}{2}(u^{n-2}-u_h^{n-2})[{(u^{n-2})}^2+u^{n-2}{u}_h^{n-2}+{(u_h^{n-2})}^2]{\parallel }^2}\hfill \\ \hfill & {\displaystyle +\parallel \frac{3}{2}(u^{n-1}-u_h^{n-1}){\parallel }^2+\parallel \frac{1}{2}(u^{n-2}-u_h^{n-2}){\parallel }^2}\hfill \\ =\hfill & {\displaystyle \frac{9}{4}(C_u^2+1)\parallel {\theta }^{n-1}+{\tau }^{n-1}{\parallel }^2+\frac{1}{4}(C_u^2+1){\parallel {\theta }^{n-2}+{\tau }^{n-2}\parallel }^2}\hfill \\ ⩽\hfill & C(\parallel {\tau }^{n-1}{\parallel }^2+\parallel {\tau }^{n-2}{\parallel }^2+\parallel {\theta }^{n-1}{\parallel }^2+\parallel {\theta }^{n-2}{\parallel }^2).\hfill \end{array}$$

(43)

Noting that

$$\parallel {\displaystyle \frac{{\tau }^n-{\tau }^{n-1}}{\mathsf{\Delta }t}}{\parallel }^2⩽{\displaystyle \frac{1}{\mathsf{\Delta }t}}{\int }_{t_{n-1}}^{t_n}{\parallel {\tau }_t(s)\parallel }^{2}ds.$$

(44)

Substitute (43) and (44) into (42) to obtain

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\gamma }{2}\parallel {\kappa }^{n-\frac{1}{2}}{\parallel }^2+\frac{1}{2\mathsf{\Delta }t}(\parallel {\theta }^n{\parallel }^2-\parallel {\theta }^{n-1}{\parallel }^2)}\hfill \\ ⩽\hfill & {\displaystyle C\left(\frac{1}{\mathsf{\Delta }t}{\int }_{t_{n-1}}^{t_n}\parallel {\tau }_t{(s)\parallel }^{2}ds+\parallel R_1^{n-\frac{1}{2}}{\parallel }^2+\parallel R_2^{n-\frac{1}{2}}{\parallel }^2+{\parallel {\pi }^{n-\frac{1}{2}}\parallel }^2\right)}\hfill \\ \hfill & +C(\parallel {\tau }^{n-1}{\parallel }^2+\parallel {\tau }^{n-2}{\parallel }^2+\parallel {\theta }^n{\parallel }^2+\parallel {\theta }^{n-1}{\parallel }^2+\parallel {\theta }^{n-2}{\parallel }^2).\hfill \end{array}$$

(45)

Summing from $n=2,\cdots ,J$, the resulting inequality becomes

$$\begin{array}{cc}\hfill & {\displaystyle \gamma \mathsf{\Delta }t\sum _{n=2}^J\parallel {\kappa }^{n-\frac{1}{2}}{\parallel }^2+(1-C\mathsf{\Delta }t){\parallel {\theta }^J\parallel }^2}\hfill \\ ⩽\hfill & {\displaystyle \parallel {\theta }^1{\parallel }^2+C\mathsf{\Delta }t\sum _{n=2}^J(\parallel R_1^{n-\frac{1}{2}}{\parallel }^2+\parallel R_2^{n-\frac{1}{2}}{\parallel }^2+\parallel {\pi }^{n-\frac{1}{2}}{\parallel }^2)}\hfill \\ \hfill & {\displaystyle +C\mathsf{\Delta }t\sum _{n=0}^{J-1}(\parallel {\tau }^n{\parallel }^2+\parallel {\theta }^n{\parallel }^2)+C{\int }_0^{t_J}{\parallel {\tau }_t\parallel }^{2}ds.}\hfill \end{array}$$

(46)

Choosing $\mathsf{\Delta }{t}_0$, when $0<\mathsf{\Delta }t<\mathsf{\Delta }{t}_0$, $1-C\mathsf{\Delta }t>0$, and using discrete Gronwall Lemma and Lemma 3, we obtain

$$\begin{array}{cc}\hfill & {\displaystyle \gamma \mathsf{\Delta }t\sum _{n=2}^J\parallel {\kappa }^{n-\frac{1}{2}}{\parallel }^2+\frac{1}{2}{\parallel {\theta }^J\parallel }^2}\hfill \\ ⩽\hfill & {\displaystyle \parallel {\theta }^1{\parallel }^2+C\mathsf{\Delta }t\sum _{n=0}^{J-1}\parallel {\tau }^n{\parallel }^2+C\mathsf{\Delta }t\sum _{n=2}^J(\parallel R_1^{n-\frac{1}{2}}{\parallel }^2+\parallel R_2^{n-\frac{1}{2}}{\parallel }^2+\parallel {\pi }^{n-\frac{1}{2}}{\parallel }^2)+C{\int }_0^{t_J}{\parallel {\tau }_t\parallel }^{2}ds}\hfill \\ ⩽\hfill & {\displaystyle \parallel {\theta }^1{\parallel }^2+C(\parallel {\tau }^n{\parallel }_{L^{\infty }(L^2)}^2+\parallel {\pi }^{n-\frac{1}{2}}{\parallel }_{L^{\infty }(L^2)}^2+{\int }_0^{t_J}\parallel {\tau }_t{\parallel }^2{ds)+C(\mathsf{\Delta }t)}^4{\parallel u_{ttt}\parallel }_{L^{\infty }(L^2)}^2.}\hfill \end{array}$$

(47)

Combining Lemma 2, (47) and the triangle inequality, the proof of (39) is accomplished.

□

Theorem 4. 

With $u_h^1=R_{h}u(t_1)$ and ${\omega }_h^1=R_h\omega (t_1)$,assume that the solution’s regular properties for (4) satisfy $u_t,{\omega }_t\in L^2(H^2)$, $u_{ttt}\in L^{\infty }(L^2)$, $u,\omega \in L^{\infty }(H^2)$. Then there exists a positive constant C independent of h and $\mathsf{\Delta }t$ such that

$$\begin{array}{c}\parallel \omega (t_J)-{\omega }_h^J\parallel ⩽C(h^2\parallel |★|\parallel +{(\mathsf{\Delta }t)}^2\parallel u_{ttt}{\parallel }_{L^{\infty }(L^2)}),\hfill \end{array}$$

(48)

where $\parallel |★|\parallel \triangleq {\parallel u\parallel }_{L^{\infty }(H^2)}+{\parallel \omega \parallel }_{L^{\infty }(H^2)}+\parallel u_t{\parallel }_{L^2(H^2)}+{\parallel {\omega }_t\parallel }_{L^2(H^2)}$.

Proof of Theorem 4. 

From (36) we can get

$${\displaystyle \left(\frac{{\kappa }^n-{\kappa }^{n-1}}{\mathsf{\Delta }t},g_h\right)-\left(\frac{\nabla {\theta }^n-\nabla {\theta }^{n-1}}{\mathsf{\Delta }t},\nabla g_h\right)=-\left(\frac{{\pi }^n-{\pi }^{n-1}}{\mathsf{\Delta }t},g_h\right),\forall g_h\in V_h.}$$

(49)

Taking ${\displaystyle v_h=\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}}$ and $g_h={\kappa }^{n-{\displaystyle \frac{1}{2}}}$ in (35) and (49), respectively, we obtain

$$\begin{array}{cc}\hfill & {\displaystyle \parallel \frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}{\parallel }^2+\gamma \left(\nabla {\kappa }^{n-\frac{1}{2}},\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)+\left({\kappa }^{n-\frac{1}{2}},\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)}\hfill \\ \hfill & {\displaystyle +\left(f({\widehat{\widehat{u}}}^{n-\frac{1}{2}})-f({\widehat{\widehat{u}}}_h^{n-\frac{1}{2}}),\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)}\hfill \\ =\hfill & {\displaystyle -\left(\frac{{\tau }^n-{\tau }^{n-1}}{\mathsf{\Delta }t},\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)-\left({\pi }^{n-\frac{1}{2}},\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)}\hfill \\ \hfill & {\displaystyle +\left(R_1^{n-\frac{1}{2}},\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)+\left(R_2^{n-\frac{1}{2}},\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right),}\hfill \end{array}$$

(50)

$${\displaystyle \frac{\parallel {\kappa }^n{\parallel }^2-{\parallel {\kappa }^{n-1}\parallel }^2}{2\mathsf{\Delta }t}-\left(\frac{\nabla {\theta }^n-\nabla {\theta }^{n-1}}{\mathsf{\Delta }t},\nabla {\kappa }^{n-\frac{1}{2}}\right)=-\left(\frac{{\pi }^n-{\pi }^{n-1}}{\mathsf{\Delta }t},{\kappa }^{n-\frac{1}{2}}\right).}$$

(51)

Multiplying (51) by $\gamma$ and adding it to (50) and using Cauchy-Schwarz inequality and Young inequality, we have

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\gamma }{2\mathsf{\Delta }t}(\parallel {\kappa }^n{\parallel }^2-\parallel {\kappa }^{n-1}{\parallel }^2)+\parallel \frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}{\parallel }^2}\hfill \\ =\hfill & {\displaystyle -\left({\kappa }^{n-\frac{1}{2}},\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)-\left(f({\widehat{\widehat{u}}}^{n-\frac{1}{2}})-f({\widehat{\widehat{u}}}_h^{n-\frac{1}{2}}),\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)}\hfill \\ \hfill & {\displaystyle -\left(\frac{{\tau }^n-{\tau }^{n-1}}{\mathsf{\Delta }t},\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)-\left({\pi }^{n-\frac{1}{2}},\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)}\hfill \\ \hfill & {\displaystyle +\left(R_1^{n-\frac{1}{2}},\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)+\left(R_2^{n-\frac{1}{2}},\frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}\right)-\gamma \left(\frac{{\pi }^n-{\pi }^{n-1}}{\mathsf{\Delta }t},{\kappa }^{n-\frac{1}{2}}\right)}\hfill \\ ⩽\hfill & {\displaystyle C(\parallel f({\widehat{\widehat{u}}}^{n-\frac{1}{2}})-f({\widehat{\widehat{u}}}_h^{n-\frac{1}{2}}){\parallel }^2+\parallel R_1^{n-\frac{1}{2}}{\parallel }^2+\parallel R_2^{n-\frac{1}{2}}{\parallel }^2+\parallel \frac{{\tau }^n-{\tau }^{n-1}}{\mathsf{\Delta }t}{\parallel }^2+\parallel \frac{{\pi }^n-{\pi }^{n-1}}{\mathsf{\Delta }t}{\parallel }^2)}\hfill \\ \hfill & {\displaystyle +C(\parallel {\pi }^{n-\frac{1}{2}}{\parallel }^2+\parallel {\kappa }^{n-\frac{1}{2}}{\parallel }^2)+\frac{1}{2}\parallel \frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}{\parallel }^2.}\hfill \end{array}$$

(52)

Similar to (44),

$$\parallel {\displaystyle \frac{{\pi }^n-{\pi }^{n-1}}{\mathsf{\Delta }t}}{\parallel }^2⩽{\displaystyle \frac{1}{\mathsf{\Delta }t}}{\int }_{t_{n-1}}^{t_n}{\parallel {\pi }_t(s)\parallel }^{2}ds.$$

(53)

Multiplying (52) by $2\mathsf{\Delta }t$, summing from $n=2,\cdots ,J$ and using (43), (44) and (53) to obtain

$$\begin{array}{cc}\hfill & {\displaystyle \gamma \parallel {\kappa }^J{\parallel }^2+\mathsf{\Delta }t\sum _{n=2}^J\parallel \frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}{\parallel }^2}\hfill \\ ⩽\hfill & {\displaystyle \gamma \parallel {\kappa }^1{\parallel }^2+C\mathsf{\Delta }t\sum _{n=2}^J(\parallel {\tau }^{n-1}{\parallel }^2+\parallel {\tau }^{n-2}{\parallel }^2+\parallel {\theta }^{n-1}{\parallel }^2+\parallel {\theta }^{n-2}{\parallel }^2+\parallel R_1^{n-\frac{1}{2}}{\parallel }^2+\parallel R_2^{n-\frac{1}{2}}{\parallel }^2)}\hfill \\ \hfill & {\displaystyle +C\mathsf{\Delta }t\sum _{n=2}^J\parallel {\pi }^{n-\frac{1}{2}}{\parallel }^2+C{\int }_0^{t_J}(\parallel {\tau }_t{(s)\parallel }^2+\parallel {\pi }_t{(s)\parallel }^2)ds+C\mathsf{\Delta }t\sum _{n=2}^J{\parallel {\kappa }^{n-\frac{1}{2}}\parallel }^2.}\hfill \end{array}$$

(54)

Using the Gronwall Lemma, we have

$$\begin{array}{cc}\hfill & {\displaystyle \gamma \parallel {\kappa }^J{\parallel }^2+\mathsf{\Delta }t\sum _{n=2}^J\parallel \frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}{\parallel }^2}\hfill \\ ⩽\hfill & {\displaystyle \gamma \parallel {\kappa }^1{\parallel }^2+C\mathsf{\Delta }t\sum _{n=0}^{J-1}(\parallel {\tau }^n{\parallel }^2+\parallel {\theta }^n{\parallel }^2)+C\mathsf{\Delta }t\sum _{n=2}^J(\parallel R_1^{n-\frac{1}{2}}{\parallel }^2+\parallel R_2^{n-\frac{1}{2}}{\parallel }^2+\parallel {\pi }^{n-\frac{1}{2}}{\parallel }^2)}\hfill \\ \hfill & {\displaystyle +C{\int }_0^{t_J}(\parallel {\tau }_t{(s)\parallel }^2+\parallel {\pi }_t(s){\parallel }^2)ds.}\hfill \end{array}$$

(55)

Substitute (47) and Lemma 3 into (55), we get

$$\begin{array}{cc}\hfill & {\displaystyle \gamma \parallel {\kappa }^J{\parallel }^2+\mathsf{\Delta }t\sum _{n=2}^J\parallel \frac{{\theta }^n-{\theta }^{n-1}}{\mathsf{\Delta }t}{\parallel }^2}\hfill \\ ⩽\hfill & {\displaystyle \gamma \parallel {\kappa }^1{\parallel }^2+C{\int }_0^{t_J}(\parallel {\tau }_t{(s)\parallel }^2+\parallel {\pi }_t{(s)\parallel }^2)ds+C(\parallel {\tau }^n{\parallel }_{L^{\infty }(L^2)}^2+\parallel {\pi }^{n-\frac{1}{2}}{\parallel }_{L^{\infty }(L^2)}^2)}\hfill \\ \hfill & {\displaystyle +C{(\mathsf{\Delta }t)}^4{\parallel u_{ttt}\parallel }_{L^{\infty }(L^2)}^2}\hfill \end{array}$$

(56)

Combining Lemma 2, (56) with the triangle inequality, we complete the proof. □

3. The reduced-dimension method of the Crank-Nicolson mixed finite element solution coefficient vectors of the EFK equation

3.1. Construction of POD basis vectors

Firstly, we need to compute the initial K-step solution coefficient vectors ${\{{\mathit{U}}_h^n,{\mathit{W}}_h^n\}}_{n=1}^K$ of Problem 3, which will yield two $M\times K$ snapshot matrices ${\mathit{A}}_1=({\mathit{U}}_h^1,{\mathit{U}}_h^2,\cdots ,{\mathit{U}}_h^K)$ and ${\mathit{A}}_2=({\mathit{W}}_h^1,{\mathit{W}}_h^2,\cdots ,{\mathit{W}}_h^K)$. Secondly, we need to compute the positive eigenvalues and arrange them in descending order ${\mu }_{i,1}⩾{\mu }_{i,2}⩾\cdots ⩾{\mu }_{i,r_i}>0(r_i=\mathrm{rank}({\mathit{A}}_i))$ of ${\mathit{A}}_i{\mathit{A}}_i^T(i=1,2)$, and the associated eigenvectors to form the matrix ${\tilde{\mathsf{\Phi }}}_i=({\mathsf{\phi }}_{i,1},{\mathsf{\phi }}_{i,2},\cdots ,{\mathsf{\phi }}_{i,r_i})\in R^{M\times r_i}$. Finally, choosing the first d column vectors of ${\tilde{\mathsf{\Phi }}}_i$ as a set of POD bases ${\mathsf{\Phi }}_i=({\mathsf{\phi }}_{i,1},{\mathsf{\phi }}_{i,2},\cdots ,{\mathsf{\phi }}_{i,d})(d⩽r_i,i=1,2)$, they satisfy the following equalities

$$\parallel {\mathit{A}}_i-{\mathsf{\Phi }}_i{\mathsf{\Phi }}_i^T{\mathit{A}}_i{\parallel }_2=\sqrt{{\mu }_{i,d+1}},i=1,2,$$

(57)

where ${\displaystyle {\parallel \mathit{A}\parallel }_2=\underset{\mathit{v}\in R^M}{sup}\frac{\parallel \mathit{A}\mathit{v}\parallel }{\parallel \mathit{v}\parallel }}$ and ${\displaystyle \parallel \mathit{v}\parallel =\sqrt{{\displaystyle \sum _{i=1}^M}{\left|v_i\right|}^2}}$ of vector $\mathit{v}={(v_1,v_2,\cdots ,v_M)}^T$. When $n=1,2,\cdots ,K$, it follows that

$$\parallel {\mathit{U}}_h^n-{\mathsf{\Phi }}_1{\mathsf{\Phi }}_1^T{\mathit{U}}_h^n\parallel =\parallel ({\mathit{A}}_1-{\mathsf{\Phi }}_1{\mathsf{\Phi }}_1^T{\mathit{A}}_1){\mathit{e}}^n\parallel ⩽\parallel ({\mathit{A}}_1-{\mathsf{\Phi }}_1{\mathsf{\Phi }}_1^T{\mathit{A}}_1){\parallel }_2\parallel {\mathit{e}}^n\parallel ⩽\sqrt{{\mu }_{1,d+1}},$$

(58)

$$\parallel {\mathit{W}}_h^n-{\mathsf{\Phi }}_2{\mathsf{\Phi }}_2^T{\mathit{W}}_h^n\parallel =\parallel ({\mathit{A}}_2-{\mathsf{\Phi }}_2{\mathsf{\Phi }}_2^T{\mathit{A}}_2){\mathit{e}}^n\parallel ⩽\parallel ({\mathit{A}}_2-{\mathsf{\Phi }}_2{\mathsf{\Phi }}_2^T{\mathit{A}}_2){\parallel }_2\parallel {\mathit{e}}^n\parallel ⩽\sqrt{{\mu }_{2,d+1}},$$

(59)

Here, ${\mathit{e}}^n(1⩽n⩽K)$ denotes the unit vectors whose nth component is 1. Therefore, ${\mathsf{\Phi }}_i=({\mathsf{\phi }}_{i,1},{\mathsf{\phi }}_{i,2},\cdots ,{\mathsf{\phi }}_{i,d})$$(d⩽r_i,i=1,2)$ is a set of optimal POD bases.

Remark 2. 

Since the order $M\times M$ of matrix ${\mathit{A}}_i{\mathit{A}}_i^T(i=1,2)$ is much larger than the order $K\times K$ of matrix ${\mathit{A}}_i^T{\mathit{A}}_i$, and the eigenvalues of ${\mathit{A}}_i{\mathit{A}}_i^T$ and ${\mathit{A}}_i^T{\mathit{A}}_i$ are equal, thus we first calculate the initial d eigenvalues ${\mu }_{i,j}(1⩽j⩽d)$ of matrix ${\mathit{A}}_i^T{\mathit{A}}_i$ and their corresponding eigenvectors ${\mathsf{\varphi }}_{i,j}(1⩽j⩽d)$. Then we use the relation ${\mathsf{\phi }}_{i,j}={\mathit{A}}_i{\mathsf{\varphi }}_{i,j}/\sqrt{{\mu }_{i,j}}(1⩽j⩽d)$ to obtain the eigenvectors of ${\mathit{A}}_i{\mathit{A}}_i^T$, and generate a set of POD basis vectors ${\mathsf{\Phi }}_i=({\mathsf{\phi }}_{i,1},{\mathsf{\phi }}_{i,2},\cdots ,{\mathsf{\phi }}_{i,d})(d⩽r_i)$.

3.2. Formulation of the RDCNMFE scheme

Firstly, letting ${\mathsf{\alpha }}_d^n={({\mathsf{\alpha }}_1^n,{\mathsf{\alpha }}_2^n,\cdots ,{\mathsf{\alpha }}_d^n)}^T$, ${\mathsf{\beta }}_d^n={({\mathsf{\beta }}_1^n,{\mathsf{\beta }}_2^n,\cdots ,{\mathsf{\beta }}_d^n)}^T$, and ${\mathit{U}}_d^n=(U_{d1}^n,U_{d2}^n,\cdots ,$ $U_{dM}^n{)}^T$ and ${\mathit{W}}_d^n={(W_{d1}^n,W_{d2}^n,\cdots ,W_{dM}^n)}^T$ represent the RDCNMFE solution coefficient vectors. Secondly, we can instantly obtain the first K RDCNMFE solution vectors ${\mathit{U}}_d^n={\mathsf{\Phi }}_1{\mathsf{\Phi }}_1^T{\mathit{U}}_h^n=:{\mathsf{\Phi }}_1{\mathsf{\alpha }}_d^n$ and ${\mathit{W}}_d^n={\mathsf{\Phi }}_2{\mathsf{\Phi }}_2^T{\mathit{W}}_h^n=:{\mathsf{\Phi }}_2{\mathsf{\beta }}_d^n(1⩽n⩽K)$ from section 3.1. Finally, by using ${\mathit{U}}_d^n={\mathsf{\Phi }}_1{\mathsf{\alpha }}_d^n,{\mathit{W}}_d^n={\mathsf{\Phi }}_2{\mathsf{\beta }}_d^n(K+1⩽n⩽N)$ instead of $\{{\mathit{U}}_h^n,{\mathit{W}}_h^n\}$ in Problem 3, the following matrix-based RDCNMFE scheme can be established.

Problem 4. 

Find $\{{\mathit{U}}_d^n,{\mathit{W}}_d^n\}\in R^M\times R^M$ and $\{u_d^n,{\omega }_d^n\}\in V_h\times V_h(1⩽n⩽N)$ satisfying

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill & {\mathsf{\alpha }}_d^n={\mathsf{\Phi }}_1^T{\mathit{U}}_h^n,{\mathsf{\beta }}_d^n={\mathsf{\Phi }}_2^T{\mathit{W}}_h^n,1⩽n⩽K,\hfill \\ \hfill & {\mathsf{\Phi }}_1{\overline{\partial }}_t{\mathsf{\alpha }}_d^n+\gamma \mathit{B}{\mathsf{\Phi }}_2{\mathsf{\beta }}_d^{n-{\displaystyle \frac{1}{2}}}+{\mathsf{\Phi }}_2{\mathsf{\beta }}_d^{n-{\displaystyle \frac{1}{2}}}=\mathit{F}({\mathsf{\Phi }}_1{\widehat{\widehat{\mathsf{\alpha }}}}_d^{n-{\displaystyle \frac{1}{2}}}),K+1⩽n⩽N,\hfill \\ \hfill & {\mathsf{\Phi }}_2{\mathsf{\beta }}_d^{n-{\displaystyle \frac{1}{2}}}-\mathit{B}{\mathsf{\Phi }}_1{\mathsf{\alpha }}_d^{n-{\displaystyle \frac{1}{2}}}=0,K+1⩽n⩽N,\hfill \\ \hfill & {\displaystyle u_d^n=\sum _{j=1}^M{U}_{dj}^n{\zeta }_j={\mathit{U}}_d^n·\zeta ={\mathsf{\Phi }}_1{\mathsf{\alpha }}_d^n·\zeta ,{\omega }_d^n=\sum _{j=1}^M{W}_{dj}^n{\zeta }_j={\mathit{W}}_d^n·\zeta ={\mathsf{\Phi }}_2{\mathsf{\beta }}_d^n·\zeta ,1⩽n⩽N,}\hfill \end{array}\hfill \end{array}\right.$$

(60)

in which $\{{\mathit{U}}_h^n,{\mathit{W}}_h^n\}(n=1,2,\cdots ,K)$ are the first K solution coefficient vectors for Problem 3. The definitions of matrix $\mathit{B}$, vector $\mathit{F}$, as well as FE basis function vectors $\zeta =({\zeta }_1(x),{\zeta }_2(x),\cdots ,{\zeta }_M(x))$ are given in section 2.2.

3.3. The uniqueness, stability and error estimate of the RDCNMFE solutions

For the RDCNMFE solutions of Problem 4, the unique, stable and convergent properties are as follows.

Theorem 5. 

Under the same assumptions as Theorem 3 and 4, if $\{u^n,{\omega }^n\}\in V\times V$ are the weak solutions to Problem 1 and $\{u_d^n,{\omega }_d^n\}\in V_h\times V_h$ are the solutions to Problem 4, then for any $1⩽n⩽N$, the RDCNMFE solutions have uniqueness and unconditional stability, and the following error estimate is obtained.

$$\begin{array}{c}\hfill \parallel u^n-u_d^n\parallel +\parallel {\omega }^n-{\omega }_d^n\parallel ⩽C(h^2+\mathsf{\Delta }{t}^2+\sqrt{{\mu }_{1,d+1}}+\sqrt{{\mu }_{2,d+1}}).\end{array}$$

(61)

Proof of Theorem 5. 

(1) Discuss the uniqueness of the RDCNMFE solutions.

• For $1⩽n⩽K$, Theorem 1 guarantees the uniqueness of solutions $\{u_h^n,{\omega }_h^n\}$ to problem 3. Therefore, the solutions $\{u_d^n,{\omega }_d^n\}$ generated from first and fourth formulas in problem 4 must also be unique.
  For $K+1⩽n⩽N$, using ${\mathit{U}}_d^n={\mathsf{\Phi }}_1{\mathsf{\alpha }}_d^n$, ${\mathit{W}}_d^n={\mathsf{\Phi }}_2{\mathsf{\beta }}_d^n(K+1⩽n⩽N)$, the last three formulas of Problem 4 can be converted to

  $${\overline{\partial }}_t{\mathit{U}}_d^n+\gamma \mathit{B}{\mathit{W}}_d^{n-{\displaystyle \frac{1}{2}}}+{\mathit{W}}_d^{n-{\displaystyle \frac{1}{2}}}=\mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}),K+1⩽n⩽N,$$

  (62)

  $${\mathit{W}}_d^{n-{\displaystyle \frac{1}{2}}}-\mathit{B}{\mathit{U}}_d^{n-{\displaystyle \frac{1}{2}}}=0,K+1⩽n⩽N,$$

  (63)

  $${\displaystyle u_d^n=\sum _{j=1}^M{U}_{dj}^n{\zeta }_j={\mathit{U}}_d^n·\zeta ,{\omega }_d^n=\sum _{j=1}^M{W}_{dj}^n{\zeta }_j={\mathit{W}}_d^n·\zeta ,K+1⩽n⩽N.}$$

  (64)
  For $K+1⩽n⩽N$, the set of solutions ${\{u_h^n,{\omega }_h^n\}}_{n=K+1}^N$ for Problem 3 is unique. Since (62) – (64) follow the same formats as problem 3, the set of solutions ${\{u_d^n,{\omega }_d^n\}}_{n=K+1}^N$ for (62) – (64) is also unique.

(2)

Demonstrate the stability of the RDCNMFE solutions.

When $1⩽n⩽K$, since the vectors in ${\mathsf{\Phi }}_1$ and ${\mathsf{\Phi }}_2$ are orthonormal, using Theorem 2, we have

$$\begin{array}{cc}\parallel u_d^n\parallel +\parallel {\omega }_d^n\parallel \hfill & =\parallel {\mathit{U}}_d^n·\zeta \parallel +\parallel {\mathit{W}}_d^n·\zeta \parallel \hfill \\ \hfill & =\parallel {\mathsf{\Phi }}_1{\mathsf{\Phi }}_1^T{\mathit{U}}_h^n·\zeta \parallel +\parallel {\mathsf{\Phi }}_2{\mathsf{\Phi }}_2^T{\mathit{W}}_h^n·\zeta \parallel \hfill \\ \hfill & ⩽C(\parallel u_h^n\parallel +\parallel {\omega }_h^n\parallel )\hfill \\ \hfill & ⩽C,1⩽n⩽L.\hfill \end{array}$$

(65)

When $K+1⩽n⩽N$, since $\mathit{B}$ is a positive definite symmetric matrix, we can rewrite (62) as

$${\mathit{B}}^{-1}{\overline{\partial }}_t{\mathit{U}}_d^n+\gamma {\mathit{W}}_d^{n-{\displaystyle \frac{1}{2}}}+{\mathit{B}}^{-1}{\mathit{W}}_d^{n-{\displaystyle \frac{1}{2}}}={\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}),K+1⩽n⩽N.$$

(66)

Putting (63) into (66), we have

$${\mathit{B}}^{-1}{\overline{\partial }}_t{\mathit{U}}_d^n+\gamma \mathit{B}{\mathit{U}}_d^{n-{\displaystyle \frac{1}{2}}}+{\mathit{U}}_d^{n-{\displaystyle \frac{1}{2}}}={\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}),K+1⩽n⩽N.$$

(67)

Taking the inner product of (67) and ${\overline{\partial }}_t{\mathit{U}}_d^n$,

$$({\mathit{B}}^{-1}{\overline{\partial }}_t{\mathit{U}}_d^n,{\overline{\partial }}_t{\mathit{U}}_d^n)+\gamma (\mathit{B}{\mathit{U}}_d^{n-{\displaystyle \frac{1}{2}}},{\overline{\partial }}_t{\mathit{U}}_d^n)+({\mathit{U}}_d^{n-{\displaystyle \frac{1}{2}}},{\overline{\partial }}_t{\mathit{U}}_d^n)=({\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}),{\overline{\partial }}_t{\mathit{U}}_d^n),K+1⩽n⩽N.$$

(68)

Then the left side of (68) is that

$$\begin{array}{cc}\hfill & {\displaystyle ({\mathit{B}}^{-1}{\overline{\partial }}_t{\mathit{U}}_d^n,{\overline{\partial }}_t{\mathit{U}}_d^n)+\gamma (\mathit{B}{\mathit{U}}_d^{n-\frac{1}{2}},{\overline{\partial }}_t{\mathit{U}}_d^n)+({\mathit{U}}_d^{n-\frac{1}{2}},{\overline{\partial }}_t{\mathit{U}}_d^n)}\hfill \\ =\hfill & {\displaystyle \parallel {({\mathit{B}}^{-1})}^{\frac{1}{2}}{\overline{\partial }}_t{\mathit{U}}_d^n{\parallel }^2+\frac{\gamma }{2\mathsf{\Delta }t}(\parallel {\mathit{B}}^{\frac{1}{2}}{\mathit{U}}_d^n{\parallel }^2-\parallel {\mathit{B}}^{\frac{1}{2}}{\mathit{U}}_d^{n-1}{\parallel }^2)}\hfill \\ \hfill & {\displaystyle +\frac{1}{2\mathsf{\Delta }t}(\parallel {\mathit{U}}_d^n{\parallel }^2-\parallel {\mathit{U}}_d^{n-1}{\parallel }^2),K+1⩽n⩽N,}\hfill \end{array}$$

(69)

and the right side of (68) is that

$$({\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}),{\overline{\partial }}_t{\mathit{U}}_d^n)⩽C\parallel {({\mathit{B}}^{-1})}^{{\displaystyle \frac{1}{2}}}{\parallel }_{\infty }^2\parallel \mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}{\parallel }^2+{\parallel {({\mathit{B}}^{-1})}^{{\displaystyle \frac{1}{2}}}{\overline{\partial }}_t{\mathit{U}}_d^n\parallel }^2,K+1⩽n⩽N.$$

(70)

Using the same technique as (22), we have

$$\parallel {({\mathit{B}}^{-1})}^{{\displaystyle \frac{1}{2}}}{\parallel }_{\infty }^2\parallel \mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}{\parallel }^2⩽C+C(\parallel {\mathit{U}}_d^{n-1}{\parallel }^2+\parallel {\mathit{U}}_d^{n-2}{\parallel }^2),K+1⩽n⩽N.$$

(71)

Combining (69), (70) and (71), we have

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\gamma }{2\mathsf{\Delta }t}(\parallel {\mathit{B}}^{\frac{1}{2}}{\mathit{U}}_d^n{\parallel }^2-\parallel {\mathit{B}}^{\frac{1}{2}}{\mathit{U}}_d^{n-1}{\parallel }^2)+\frac{1}{2\mathsf{\Delta }t}(\parallel {\mathit{U}}_d^n{\parallel }^2-\parallel {\mathit{U}}_d^{n-1}{\parallel }^2)}\hfill \\ ⩽\hfill & {\displaystyle C+C(\parallel {\mathit{U}}_d^{n-1}{\parallel }^2+\parallel {\mathit{U}}_d^{n-2}{\parallel }^2),K+1⩽n⩽N.}\hfill \end{array}$$

(72)

Multiplying both sides of (72) by $2\mathsf{\Delta }t$, and summating from 2 to n, we obtain

$$\begin{array}{cc}{\displaystyle \parallel {\mathit{U}}_d^n{\parallel }^2}\hfill & {\displaystyle ⩽\parallel {\mathit{U}}_d^1{\parallel }^2+\gamma \parallel {\mathit{B}}^{\frac{1}{2}}{\mathit{U}}_d^1{\parallel }^2+2\mathsf{\Delta }t\sum _{i=2}^{n}C+C\mathsf{\Delta }t\sum _{i=0}^{n-1}{\parallel {\mathit{U}}_d^i\parallel }^2}\hfill \\ \hfill & {\displaystyle ⩽(1+\gamma \parallel {\mathit{B}}^{\frac{1}{2}}{\parallel }_{\infty }^2)\parallel {\mathit{U}}_d^1{\parallel }^2+CT+C\mathsf{\Delta }t\sum _{i=0}^{n-1}{\parallel {\mathit{U}}_d^i\parallel }^2,K+1⩽n⩽N.}\hfill \end{array}$$

(73)

Noting that

$$\parallel {\mathit{U}}_d^1{\parallel }^2=\parallel {\mathsf{\Phi }}_1{\mathsf{\Phi }}_1^T{\mathit{U}}_h^1{\parallel }^2⩽C{\parallel {\mathit{U}}_h^1\parallel }^2⩽C,$$

(74)

putting (74) into (73), we have

$$\parallel {\mathit{U}}_d^n{\parallel }^2⩽C+C\mathsf{\Delta }t\sum _{i=0}^{n-1}{\parallel {\mathit{U}}_d^i\parallel }^2,K+1⩽n⩽N.$$

(75)

Using the Gronwall Lemma for (75),

$$\parallel {\mathit{U}}_d^n{\parallel }^2⩽C{e}^{Cn\mathsf{\Delta }t}⩽C,K+1⩽n⩽N.$$

(76)

And

$$\parallel {\mathit{W}}_d^n\parallel =\parallel \mathit{B}{\mathit{U}}_d^n{\parallel ⩽\parallel \mathit{B}\parallel }_{\infty }\parallel {\mathit{U}}_d^n\parallel ⩽C\parallel {\mathit{U}}_d^n\parallel ⩽C,K+1⩽n⩽N.$$

(77)

So we have

$$\parallel {\mathit{U}}_d^n\parallel +\parallel {\mathit{W}}_d^n\parallel ⩽C,K+1⩽n⩽N.$$

(78)

Thus, noting that $\parallel \zeta \parallel ⩽C$, we get

$$\parallel u_d^n\parallel +\parallel {\omega }_d^n\parallel =\parallel {\mathit{U}}_d^n·\zeta \parallel +\parallel {\mathit{W}}_d^n·\zeta \parallel ⩽C\parallel {\mathit{U}}_d^n\parallel ·\parallel \zeta \parallel +C\parallel {\mathit{W}}_d^n\parallel ·\parallel \zeta \parallel ⩽C,K+1⩽n⩽N.$$

(79)

From (65) and (79), the RDCNMFE solutions $\{u_d^n,{\omega }_d^n\}(1⩽n⩽N)$ are unconditional stable.

(3)

Analyse the error estimate of the RDCNMFE solutions.

When $1⩽n⩽K$, from $\parallel \zeta \parallel ⩽C$, (58) and (59), we have

$$\begin{array}{cc}\parallel u_h^n-u_d^n\parallel +\parallel {\omega }_h^n-{\omega }_d^n\parallel \hfill & ⩽\parallel {\mathit{U}}_h^n-{\mathit{U}}_d^n{\parallel }_{\infty }\parallel \zeta \parallel +\parallel {\mathit{W}}_h^n-{\mathit{W}}_d^n{\parallel }_{\infty }\parallel \zeta \parallel \hfill \\ \hfill & ⩽C\parallel {\mathit{U}}_h^n-{\mathsf{\Phi }}_1{\mathsf{\Phi }}_1^T{\mathit{U}}_h^n\parallel +\parallel {\mathit{W}}_h^n-{\mathsf{\Phi }}_2{\mathsf{\Phi }}_2^T{\mathit{W}}_h^n\parallel \hfill \\ \hfill & ⩽C(\sqrt{{\mu }_{1,d+1}}+\sqrt{{\mu }_{2,d+1}}),1⩽n⩽K.\hfill \end{array}$$

(80)

When $K+1⩽n⩽N$, letting ${\delta }^n={\mathit{U}}_h^n-{\mathit{U}}_d^n$ and ${\rho }^n={\mathit{W}}_h^n-{\mathit{W}}_d^n$, and combining (17), (66) and (63), we obtain

$${\mathit{B}}^{-1}{\overline{\partial }}_t{\delta }^n+\gamma {\rho }^{n-{\displaystyle \frac{1}{2}}}+{\mathit{B}}^{-1}{\rho }^{n-{\displaystyle \frac{1}{2}}}={\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-{\displaystyle \frac{1}{2}}})-{\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}),K+1⩽n⩽N,$$

(81)

$${\rho }^{n-{\displaystyle \frac{1}{2}}}=\mathit{B}{\delta }^{n-{\displaystyle \frac{1}{2}}},K+1⩽n⩽N.$$

(82)

Putting (82) into(81), we have

$${\mathit{B}}^{-1}{\overline{\partial }}_t{\delta }^n+\gamma \mathit{B}{\delta }^{n-{\displaystyle \frac{1}{2}}}+{\delta }^{n-{\displaystyle \frac{1}{2}}}={\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-{\displaystyle \frac{1}{2}}})-{\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}),K+1⩽n⩽N.$$

(83)

Taking the inner product of (83) and ${\overline{\partial }}_t{\delta }^n$,

$$\begin{array}{cc}\hfill & ({\mathit{B}}^{-1}{\overline{\partial }}_t{\delta }^n,{\overline{\partial }}_t{\delta }^n)+\gamma (\mathit{B}{\delta }^{n-{\displaystyle \frac{1}{2}}},{\overline{\partial }}_t{\delta }^n)+({\delta }^{n-{\displaystyle \frac{1}{2}}},{\overline{\partial }}_t{\delta }^n)\hfill \\ =\hfill & ({\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-{\displaystyle \frac{1}{2}}})-{\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}),{\overline{\partial }}_t{\delta }^n),K+1⩽n⩽N.\hfill \end{array}$$

(84)

Then the left side of (84) is that

$$\begin{array}{cc}\hfill & {\displaystyle ({\mathit{B}}^{-1}{\overline{\partial }}_t{\delta }^n,{\overline{\partial }}_t{\delta }^n)+\gamma (\mathit{B}{\delta }^{n-\frac{1}{2}},{\overline{\partial }}_t{\delta }^n)+({\delta }^{n-\frac{1}{2}},{\overline{\partial }}_t{\delta }^n)}\hfill \\ =\hfill & {\displaystyle \parallel {({\mathit{B}}^{-1})}^{\frac{1}{2}}{\overline{\partial }}_t{\delta }^n{\parallel }^2+\frac{\gamma }{2\mathsf{\Delta }t}(\parallel {\mathit{B}}^{\frac{1}{2}}{\delta }^n{\parallel }^2-\parallel {\mathit{B}}^{\frac{1}{2}}{\delta }^{n-1}{\parallel }^2)}\hfill \\ \hfill & {\displaystyle +\frac{1}{2\mathsf{\Delta }t}(\parallel {\delta }^n{\parallel }^2-\parallel {\delta }^{n-1}{\parallel }^2),K+1⩽n⩽N,}\hfill \end{array}$$

(85)

and the right side of (84) is that

$$\begin{array}{cc}\hfill & ({\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-{\displaystyle \frac{1}{2}}})-{\mathit{B}}^{-1}\mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}),{\overline{\partial }}_t{\delta }^n)\hfill \\ ⩽\hfill & C\parallel {({\mathit{B}}^{-1})}^{{\displaystyle \frac{1}{2}}}{\parallel }_{\infty }^2\parallel \mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-{\displaystyle \frac{1}{2}}})-\mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-{\displaystyle \frac{1}{2}}}){\parallel }^2+{\parallel {({\mathit{B}}^{-1})}^{{\displaystyle \frac{1}{2}}}{\overline{\partial }}_t{\delta }^n\parallel }^2,K+1⩽n⩽N.\hfill \end{array}$$

(86)

Combining Lemma 1 with the global Lipschitz-continuity of $f(u)$, the first term of (86) can be estimated as

$$\begin{array}{cc}\hfill & {\displaystyle \parallel {({\mathit{B}}^{-1})}^{\frac{1}{2}}{\parallel }_{\infty }^2{\parallel \mathit{F}({\widehat{\widehat{\mathit{U}}}}_h^{n-\frac{1}{2}})-\mathit{F}({\widehat{\widehat{\mathit{U}}}}_d^{n-\frac{1}{2}})\parallel }^2}\hfill \\ ⩽\hfill & {\displaystyle C\parallel \left(\frac{3}{2}\mathit{F}({\mathit{U}}_h^{n-1})-\frac{1}{2}\mathit{F}({\mathit{U}}_h^{n-2})\right)-\left(\frac{3}{2}\mathit{F}({\mathit{U}}_d^{n-1})-\frac{1}{2}\mathit{F}({\mathit{U}}_d^{n-2})\right){\parallel }^2}\hfill \\ ⩽\hfill & {\displaystyle C\parallel \frac{3}{2}({\mathit{U}}_h^{n-1}-{\mathit{U}}_d^{n-1}){\parallel }^2+\parallel \frac{1}{2}({\mathit{U}}_h^{n-2}-{\mathit{U}}_d^{n-2}){\parallel }^2}\hfill \\ ⩽\hfill & {\displaystyle C(\parallel {\delta }^{n-1}{\parallel }^2+\parallel {\delta }^{n-2}{\parallel }^2),K+1⩽n⩽N.}\hfill \end{array}$$

(87)

Combining (85), (86) and (87), we have

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\gamma }{2\mathsf{\Delta }t}(\parallel {\mathit{B}}^{\frac{1}{2}}{\delta }^n{\parallel }^2-\parallel {\mathit{B}}^{\frac{1}{2}}{\delta }^{n-1}{\parallel }^2)+\frac{1}{2\mathsf{\Delta }t}(\parallel {\delta }^n{\parallel }^2-\parallel {\delta }^{n-1}{\parallel }^2)}\hfill \\ ⩽\hfill & C(\parallel {\delta }^{n-1}{\parallel }^2+\parallel {\delta }^{n-2}{\parallel }^2),K+1⩽n⩽N.\hfill \end{array}$$

(88)

Multiplying both sides of (88) by $2\mathsf{\Delta }t$, and summating from $K+1$ to $n(n⩽N)$, we obtain

$$\begin{array}{cc}\parallel {\delta }^n{\parallel }^2\hfill & {\displaystyle ⩽\parallel {\delta }^K{\parallel }^2+\parallel {\mathit{B}}^{\frac{1}{2}}{\delta }^K{\parallel }^2+C\mathsf{\Delta }t\sum _{i=K-1}^{n-1}{\parallel {\delta }^i\parallel }^2}\hfill \\ \hfill & {\displaystyle ⩽C\parallel {\delta }^K{\parallel }^2+C\mathsf{\Delta }t\sum _{i=K-1}^{n-1}{\parallel {\delta }^i\parallel }^2,K+1⩽n⩽N.}\hfill \end{array}$$

(89)

Noting that

$$\parallel {\delta }^K{\parallel }^2=\parallel {\mathit{U}}_h^K-{\mathit{U}}_d^K{\parallel }^2={\parallel {\mathit{U}}_h^K-{\mathsf{\Phi }}_1{\mathsf{\Phi }}_1^T{\mathit{U}}_h^K\parallel }^2.$$

(90)

Putting (90) into (89), from (58) and (59), we have

$${\displaystyle \parallel {\delta }^n{\parallel }^2⩽C{\mu }_{1,d+1}+C\mathsf{\Delta }t\sum _{i=K-1}^{n-1}{\parallel {\delta }^i\parallel }^2,K+1⩽n⩽N.}$$

(91)

Applying the Gronwall Lemma for (91),

$$\parallel {\delta }^n{\parallel }^2⩽C{\mu }_{1,d+1}{e}^{Cn\mathsf{\Delta }t}⩽C{\mu }_{1,d+1},K+1⩽n⩽N.$$

(92)

And

$$\parallel {\rho }^n\parallel =\parallel \mathit{B}{\delta }^n{\parallel ⩽C\parallel \mathit{B}\parallel }_{\infty }\parallel {\delta }^n\parallel ⩽C\parallel {\delta }^n\parallel ,K+1⩽n⩽N,$$

(93)

thus, we get

$$\parallel {\delta }^n\parallel +\parallel {\rho }^n\parallel ⩽C(\sqrt{{\mu }_{1,d+1}}+\sqrt{{\mu }_{2,d+1}}),K+1⩽n⩽N.$$

(94)

Further, from $\parallel \zeta \parallel ⩽C$, we obtain

$$\begin{array}{cc}\parallel u_h^n-u_d^n\parallel +\parallel {\omega }_h^n-{\omega }_d^n\parallel \hfill & ⩽\parallel {\mathit{U}}_h^n-{\mathit{U}}_d^n{\parallel }_{\infty }·\parallel \zeta \parallel +\parallel {\mathit{W}}_h^n-{\mathit{W}}_d^n{\parallel }_{\infty }·\parallel \zeta \parallel \hfill \\ \hfill & ⩽\parallel {\mathit{U}}_h^n-{\mathit{U}}_d^n\parallel ·\parallel \zeta \parallel +\parallel {\mathit{W}}_h^n-{\mathit{W}}_d^n\parallel ·\parallel \zeta \parallel \hfill \\ \hfill & ⩽C(\sqrt{{\mu }_{1,d+1}}+\sqrt{{\mu }_{2,d+1}}),K+1⩽n⩽N.\hfill \end{array}$$

(95)

Using the triangle inequality, Theorem 3 and 4, (80) and (95), we obtain

$$\begin{array}{cc}\parallel u^n-u_d^n\parallel +\parallel {\omega }^n-{\omega }_d^n\parallel \hfill & ⩽\parallel u^n-u_h^n\parallel +\parallel u_h^n-u_d^n\parallel +\parallel {\omega }^n-{\omega }_h^n\parallel +\parallel {\omega }_h^n-{\omega }_d^n\parallel \hfill \\ \hfill & ⩽C(h^2+\mathsf{\Delta }{t}^2+\sqrt{{\mu }_{1,d+1}}+\sqrt{{\mu }_{2,d+1}}),1⩽n⩽N.\hfill \end{array}$$

(96)

 □

4. The Numerical Examples for the EFK Equations

In order to verify the effectiveness of the RDCNMFE method, we can compare the errors and orders between the genuine, CNMFE and RDCNMFE solutions and CPU runtime on different time nodes. We here adopt the numerical experiment that the following 2D nonlinear EFK equation has a genuine solution. Generally speaking, it has no genuine solution.

$$\left\{\begin{array}{cc}{u}_t+{\mathsf{\Delta }}^{2}u-\mathsf{\Delta }u+u^3-u=g(x_1,x_2,t),\hfill & \hfill (x_1,x_2,t)\in {(0,1)}^2\times (0,T],\\ u(0,x_2,t)=u(1,x_2,t)=0,u(x_1,0,t)=u(x_1,1,t)=0,\hfill & \hfill t\in [0,T],\\ \mathsf{\Delta }u(0,x_2,t)=\mathsf{\Delta }u(1,x_2,t)=0,\mathsf{\Delta }u(x_1,0,t)=\mathsf{\Delta }u(x_1,1,t)=0,\hfill & \hfill t\in [0,T],\\ u(x_1,x_2,0)=\sin (\pi x_1)\sin (\pi x_2),\hfill & \hfill (x_1,x_2)\in {[0,1]}^2,\end{array}\right.$$

(97)

in which $g(x_1,x_2,t)=e^{-t}(-2+4{\pi }^4+2{\pi }^2)\sin (\pi x_1)\sin (\pi x_2)+e^{-3t}{\sin }^3(\pi x_1){\sin }^3(\pi x_2)$, we may obtain the genuine solution $u(x_1,x_2,t)={\mathrm{e}}^{-t}\sin (\pi x_1)\sin (\pi x_2)$ of (97). The genuine solution of the auxiliary variable $\omega$ is $\omega (x_1,x_2,t)=2{\pi }^2{\mathrm{e}}^{-t}\sin (\pi x_1)\sin (\pi x_2)$.

Firstly, the initial 20 CNMFE solution vectors $\{{\mathit{U}}_h^n,{\mathit{W}}_h^n\}(n=1,2,\cdots ,20)$ are obtained by calculating Problem 3 to establish the snapshot matrixes ${\mathit{A}}_1=({\mathit{U}}_h^1,{\mathit{U}}_h^2,\cdots ,{\mathit{U}}_h^{20})$ and ${\mathit{A}}_2=({\mathit{W}}_h^1,{\mathit{W}}_h^2,\cdots ,{\mathit{W}}_h^{20})$. Secondly, we calculate the eigenvalues ${\mu }_{i,j}(i=1,2,j=1,2,\cdots ,20)$ (arranged decreasingly) and the associated eigenvectors ${\mathsf{\varphi }}_{i,j}(i=1,2,j=1,2,\cdots ,20)$ of the matrixes ${\mathit{A}}_i^T{\mathit{A}}_i(i=1,2)$. Then, we find that $\sqrt{{\mu }_{i,7}}⩽h^2+\mathsf{\Delta }{t}^2$, so that we just need to extract the first 6 eigenvectors ${\mathsf{\varphi }}_{i,j}$ of the matrix ${\mathit{A}}_i^T{\mathit{A}}_i(i=1,2)$ to construct a set of POD basis vectors ${\mathsf{\Phi }}_i=({\mathsf{\phi }}_{i,1},{\mathsf{\phi }}_{i,2},\cdots ,{\mathsf{\phi }}_{i,6})$, where ${\mathsf{\phi }}_{i,j}={\mathit{A}}_i{\mathsf{\varphi }}_{i,j}/\sqrt{{\mu }_{i,j}}(i=1,2,j=1,2,\cdots ,6)$. Finally, we compute the RDCNMFE solutions $\{u_d^n,{\omega }_d^n\}$ of (97) by Problem 4 and the CNMFE solutions $\{u_h^n,{\omega }_h^n\}$ of (97) by Problem 3 when $h=\sqrt{2}/100$ and $\mathsf{\Delta }t=1/100$ at $t=1$. And they are compared with the genuine solution respectively, exhibited in Figure 1 and Figure 2. As evident from Figure 1 and Figure 2, the RDCNMFE and CNMFE solutions closely mirror the genuine solution.

For comparison, we give the errors and convergence orders of the CNMFE and RDCNMFE solutions for u and $\omega$ at $t=1$ in Table 1 and Table 2, respectively. As can be observed from Table 1 and Table 2, the errors and convergence rates of the RDCNMFE method are nearly identical to the CNMFE method, which aligns with the previous theoretical analysis. Figure 3 offers a more visual representation of the comparative results of the error convergence orders between $u_h^n$ and $u_d^n$ and between ${\omega }_h^n$ and ${\omega }_d^n$.

To provide additional evidence for the effectiveness of the RDCNMFE method, we recorded the errors and the CPU running time for computing the CNMFE and RDCNMFE solutions at $t=1.0, 1.5, 2.0, 2.5, 3.0$, respectively. These results are obtained by solving Problems 3 and 4 and shown in Table 3.

As demonstrated in Table 3, with the increase of time, the CNMFE method with $2\times 6561$ degrees of freedom requires considerable CPU runtime at each time node, where the running time increases by approximately $30s$ for every additional $0.5s$. While performing at the same time node, the RDCNMFE method, with only $2\times 6$ degrees of freedom, requires very minimal CPU runtime, with an increase in running time of only about $0.1s$ for every additional $0.5s$. The calculation time of the CNMFE method is approximately 70 times greater than that of the RDCNMFE method when $t=3.0$. This experiment demonstrates that the RDCNMFE method significantly reduces CPU running time. Furthermore, the results from the Table 3 indicate that the RDCNMFE method produces slightly less $L^2$ errors compared to the CNMFE method. In summary, it is clear that the RDCNMFE method surpasses the CNMFE method, proving it as an effective numerical method for solving the 2D nonlinear EFK equation.

5. Conclusions

In this study, we have focused on the reduced-dimension method of Crank-Nicolson mixed finite element solution coefficient vectors for 2D nonlinear EFK equations In order to give a full demonstration of the effectiveness of the reduced-dimension method, we have established the CNMFE scheme of the 2D nonlinear EFK equation and discussed the uniqueness, stability, and convergence of CNMFE solution. Then, we have generated the POD basis based on the initial K CNMFE solutions, constructed the RDCNMFE matrix model, and employed the traditional finite element techniques to analyze the uniqueness, stability and error estimation of RDCNMFE solution. Additionally, numerical experiments have been conducted to comparatively evaluate the performance of the two methods. Compared with the CNMFE method, the RDCNMFE method has fewer degrees of freedom at each time node. This characteristic enables the RDCNMFE method to greatly reduce computational load, running time, and error accumulation. In particular, this paper has simplified the actual calculation by linearizing the nonlinear term without iterative numerical calculation. Therefore, the RDCNMFE method constitutes a novel numerical technique for tackling nonlinear PDEs.

In future research, we intend to extend the scope of application of this method to PDEs with variable coefficients (in this paper, the coefficient $\gamma$ is a constant). Additionally, we believe this approach can be adapted to more complex engineering problems.