An Unchanged Basis Function and Preserving Accuracy Crank–Nicolson Finite Element Reduced-Dimension Method for Symmetric Tempered Fractional Diffusion Equation

Abstract

We herein mainly employ a proper orthogonal decomposition (POD) to study the reduced dimension of unknown solution coefficient vectors in the Crank–Nicolson finite element (FE) (CNFE) method for the symmetric tempered fractional diffusion equation so that we can build the reduced-dimension recursive CNFE (RDRCNFE) method. In this case, the RDRCNFE method keeps the same basic functions and accuracy as the CNFE method. Especially, we adopt the matrix analysis to discuss the stability and convergence of RDRCNFE solutions, resulting in the very laconic theoretical analysis. We also use some numerical simulations to confirm the correctness of theoretical results.

1. Introduction

The fractional partial differential equations (PDEs) are of significant physical meaning and can be used to depict many natural phenomena. In addition, the fractional PDEs frequently come out in biology and economics, too. Thereby, the fractional PDEs have been attracting a lot of attention. For instance, Ding [1] established a high-order finite difference (FD) numerical scheme for a two-dimensional (2D) time-space tempered fractional diffusion-wave equation. Du et al. [2] addressed a high-order FD algorithm for the fractional diffusion wave equation with the Caputo fractional derivative. Li and Ding [3] also posed the higher order FD scheme for the reaction and anomalous diffusion equation. Xing and Wen [4] developed a fourth-order FD algorithm for the 2D space-fractional diffusion equations. Zhou et al. [5,6] also built some FD schemes for space fractional diffusion equations. Luo and Wang [7] established a reduced-order FD scheme for the fractional-order parabolic-type sine-Gordon equations, and Zhou and Luo [8] founded an optimized FD algorithm for the fractional-order parabolic-type sine-Gordon equations. Especially, Çelik and Duman [9] constructed the Crank–Nicolson finite element (FE) (CNFE) method with the unconditionally stable second-order time accuracy for the symmetric tempered fractional diffusion equation (STFDE), which is one of the most effective FE numerical methods.

Herein, we mainly study the reduced-dimension of unknown solution coefficient vectors to the CNFE method of the following STFDE in [9].

Problem 1.

Seek $\varpi :\Omega \times J\to \mathbb{R}$ that satisfies

$$\left\{\begin{array}{cc}{\displaystyle {\partial }_t\varpi (x,t)+{\partial }_{\left|x\right|}^{\theta ,\gamma }\varpi (x,t)=\rho (x,t),}\hfill & x\in \Omega ,t\in J,\hfill \\ {\displaystyle \varpi (0,t)=\varpi (l,t)=0,}\hfill & t\in \overline{J},\hfill \\ {\displaystyle \varpi (x,0)={\varpi }^0\left(x\right),}\hfill & x\in \overline{\Omega },\hfill \end{array}\right.$$

(1)

in which $\Omega =(0,l)$, $l>0$ is a real numer, $J=(0,T)$, T is the time upper limit, ${\partial }_t=\partial /\partial t$, $1<\theta <2$, $\gamma >0$, $\rho (x,t)$ and ${\varpi }^0\left(x\right)$ are the known source and initial functions, respectively, ${\partial }_{\left|x\right|}^{\theta ,\gamma }\varpi (x,t)$ is defined as follows

$$\begin{array}{c}{\displaystyle {\partial }_{\left|x\right|}^{\theta ,\gamma }\varpi (x,t)=\frac{-1}{2\cos (\theta \pi /2)}[D_{+}^{\theta ,\gamma }\varpi (x,t)+D_{-}^{\theta ,\gamma }\varpi (x,t)-2{\gamma }^{\theta }\varpi (x,t)],}\hfill \\ {\displaystyle D_{+}^{\theta ,\gamma }\varpi (x,t)=\frac{e^{-\gamma x}}{\Gamma (n-\theta )}\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}{\int }_{-\infty }^x{(x-\xi )}^{n-\theta -1}{e}^{\gamma \xi }\varpi (\xi ,t)\mathrm{d}\xi ,}\hfill \\ {\displaystyle D_{-}^{\theta ,\gamma }\varpi (x,t)=\frac{{(-1)}^n{e}^{\gamma x}}{\Gamma (n-\theta )}\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}{\int }_{-\infty }^x{(\xi -x)}^{n-\theta -1}{e}^{-\gamma \xi }\varpi (\xi ,t)\mathrm{d}\xi ,}\hfill \\ {\displaystyle \Gamma \left(\theta \right)={\int }_0^{\infty }{x}^{\theta -1}{e}^{-x}\mathrm{d}x,n=\left[\theta \right]+1,}\hfill \end{array}$$

and $\left[\theta \right]$ is the integer part for θ.

Although the CNFE method with the second-order time accuracy in [9] is one of the most effective numerical methods, when the partition on the region $\overline{\Omega }$ is sufficiently refined, it also includes lots of unknowns, which would bring many difficulties in practical application. Hence, a key task is to reduce the unknowns of the CNFE method so as to lessen CPU runtime and rounding error amassing, and to mitigate the calculated load in the calculated process.

A lot of numerical simulations (see [4,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]) have shown that the POD method is one of the most effective approaches to reduce the unknowns of numerical models. It has played an important role in the order reduction of numerical models, such as the FD scheme, the Galerkin method, the FE method, the finite volume element (FVE) method, and the reduced basis (RB) method for time-space integer order derivative PDEs (see [4,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]). However, the above reduced-order FE, FVE, and RB methods all lower the dimension of subspaces of approximate solutions by the POD technique.

Unfortunately, at the moment, the POD technique has not been used to reduce the dimension of unknown solution coefficient vectors to the CNFE method of STFDE with the spatial fractional-order derivative. Hence, we herein adopt the POD technique to lower the dimension of unknown CNFE solution coefficient vectors of STFDE so as to build the reduced-dimension recursive CNFE (RDRCNFE) method with very few unknowns.

It is worth noting that the RDRCNFE method is absolutely different from the existed reduced-order FE and CNFE methods or the reduced-order FVE and RB methods both theoretically and technically, and the RDRCNFE method has at least the following two aspects of advantages.

Firstly, the existed reduced-order FE and CNFE methods or the reduced-order FVE and RB methods were built by replacing the finite dimensional subspaces with the subspaces spanned by the continuous POD basic functions, where their accuracy is impacted by the reduction order; unlike those, the RDRCNFE method is established by using the POD technique to lower the dimension of unknown solution coefficient vectors after the basic functions in the CNFE method are absorbed into the stiffness matrix, so that the RDRCNFE method possesses the same FE subspace and accuracy as the CNFE method. Although the unknowns of the RDRCNFE method are greatly reduced, the FE basis functions of the RDRCNFE method remain unchanged so that the accuracy of the RDRCNFE method is unchanged and maintains the same as that of the CNFE method. More specifically, if we assume that the finite dimensional subspace in the FE and CNFE methods or the FVE and RB methods is as follows

$$\begin{array}{c}\hfill {\displaystyle {\mathbb{W}}_h=\mathrm{span}\left\{{\zeta }_j\left(x\right):1⩽j⩽M\right\},}\end{array}$$

where ${\left\{{\zeta }_j\left(x\right)\right\}}_{j=1}^M$ is the set of basis functions, then the classical solution ${\omega }_h^n$ in the FE and CNFE methods or the FVE, CS, and RB methods should be expressed as

$$\begin{array}{c}{\displaystyle {\omega }_h^n={\zeta }_1\left(x\right)z_1^n+{\zeta }_2\left(x\right)z_2^n+\cdots +{\zeta }_M\left(x\right)z_M^n}\hfill \\ =({\zeta }_1\left(x\right),{\zeta }_2\left(x\right),\cdots ,{\zeta }_M\left(x\right))·(z_1^n,z_2^n,\cdots ,z_M^n)\in {\mathbb{W}}_h,1⩽n⩽N.\hfill \end{array}$$

Then, the existed reduced-order FE and CNFE methods or the reduced-order FVE and RB methods are to lower the dimension of ${\mathbb{W}}_h$, namely the dimension of basis function vector $({\zeta }_1\left(x\right),{\zeta }_2\left(x\right),\cdots ,{\zeta }_M\left(x\right))$; in other words, the subspace ${\mathbb{W}}_h$ is replaced with the d-dimensional subspace ${\mathbb{W}}_d=$ span$\{{\phi }_1,{\phi }_2,\cdots ,{\phi }_d\}$ generated by the main fewer d (usually d = 5∼7) POD basis ${\phi }_1,{\phi }_2,\cdots ,{\phi }_d$. In addition, the classical solutions ${\displaystyle {\omega }_h^n={\zeta }_1\left(x\right)z_1^n+{\zeta }_2\left(x\right)z_2^n+\cdots +{\zeta }_M\left(x\right)z_M^n}$ are approximated with ${\omega }_d^n={\phi }_1{\beta }_1^n+{\phi }_2{\beta }_2^n+\cdots +{\phi }_d{\beta }_d^n$, where ${\beta }_i$$(1⩽i⩽d)$ are unknowns. Whereas, the RDRCNFE method herein is to lower the dimension of unknown solution coefficient vectors $(z_1^n,z_2^n,\cdots ,z_M^n)$$(1⩽n⩽N)$, namely the classical solutions ${\displaystyle {\omega }_h^n=({\zeta }_1\left(x\right),{\zeta }_2\left(x\right),\cdots ,{\zeta }_M\left(x\right))·(z_1^n,z_2^n,\cdots ,z_M^n)}$ are approximated with the linear combination ${\varpi }_d^n=({\zeta }_1\left(x\right),{\zeta }_2\left(x\right),\cdots ,{\zeta }_M\left(x\right))·{\left[{\Psi }_{M\times d}(b_1^n,b_2^n,\cdots ,b_d^n)\right]}_{M\times 1}^T$ of the POD basis vectors ${\Psi }_{M\times d}$ generated by the first few solution coefficient vectors $(z_1^n,z_2^n,\cdots ,z_M^n)$$(1⩽n⩽L\ll N)$, where $b_i$$(1⩽i⩽d)$ are unknowns. Thus, the FE basis functions in the RDRCNFE solutions ${\varpi }_d^n$ are unchanged and are still ${\left\{{\zeta }_j\left(x\right)\right\}}_{j=1}^M$ so that the RDRCNFE solutions ${\varpi }_d^n$ have the same accuracy as the classical solutions ${\displaystyle {\omega }_h^n}$. Therefore, the reduced dimension of unknown solution coefficient vectors herein is completely different from the existed reduced dimension of finite dimensional subspace, i.e., basis function vector.

Secondly, we adopt the matrix analysis to analyze the stability and convergence of RDRCNFE solutions, resulting in the very laconic and readily understood theory method, but the stability and convergence of the existed reduced-order FE and CNFE solutions are discussed by functional analysis so that their theory methods are abstract and complicated. Thereby, the RDRCNFE method is completely new.

Although the reduced-order methods of the unknown solution coefficient vectors for the hyperbolic, parabolic, Sobolev, viscoelastic wave, and unsteady Stokes equations with time-space integer order derivatives have been built in [27,28,29,30,31], respectively, the STFDE with the spatial fractional-order derivative herein is more complicated than the above five types of equations, so that both the structure for the RDRCNFE model and the theoretical analysis of existence, stability, and errors for the RDRCNFE solutions need more techniques and have more difficulties than those in [27,28,29,30,31]. However, the RDRCNFE method for the STFDE with the spatial fractional order derivative has very important applications.

The rest of the content herein is arranged in the following four sections. In Section 2, we retrospect the CNFE method of STFDE in [9] and the associated theoretical results such as the existence as well as the stability together with the error estimates for the CNFE solutions, and the most key step is to rewrite the CNFE format of functional form into the matrix form. In Section 3, we first build the RDRCNFE method with a set of POD basic vectors generated by the initial several CNFE solution coefficient vectors, and then, we employ the matrix analysis to analyze the stability together with the convergence of the RDRCNFE solutions. In Section 4, we utilize some numerical simulations to confirm the rightness of theory results so as to reveal the superiority of the RDRCNFE method. We provide the main conclusions and discussion in Section 5.

2. Retrospect the CNFE Method for STFDE and Rewrite Matrix-Form

2.1. Retrospect the CNFE Method for STFDE

Herein, the integer and fractional order Sobolev spaces (see [32]) are simultaneously adopted.

For $\theta >0$ and $\gamma >0$, set

$$H^{\theta ,\gamma }\left(\mathbb{R}\right)=\{v\in L^2\left(\mathbb{R}\right):{\left({\gamma }^2+{\omega }^2\right)}^{\theta /2}\widehat{v}\in L^2\left(\mathbb{R}\right)\},$$

endowed with norm

$${\parallel v\parallel }_{H^{\theta ,\gamma }\left(\mathbb{R}\right)}={\parallel {({\gamma }^2+{\omega }^2)}^{\theta /2}\widehat{v}\parallel }_{L^2\left(\mathbb{R}\right)},$$

where $\widehat{v}$ is Fourier’s transformation of the function v, and $\omega$ is the variable in Fourier’s transformation.

Let $\mathbb{W}=H^1(0,T;H_0^{\frac{\theta }{2},\gamma }$$(\Omega ))$, $H_0^{\theta ,\gamma }(\Omega )=\left\{v\in H^{\theta ,\gamma }{(\Omega ):v|}_{\partial \Omega }=0\right\}$, and $H^{\theta ,\gamma }(\Omega )$ consist of the restrictions in $\Omega =(0,l)$ of functions in the space $H^{\theta ,\gamma }\left(R\right)$. Thereupon, with integration by parts, the weak form for STFDE is established in the following.

Problem 2.

Seek $u\in \mathbb{W}$ that satisfies

$$\begin{array}{c}({\varpi }_t,\upsilon )+a(\varpi ,\upsilon )=(\rho ,\upsilon ),\forall \upsilon \in H_0^{\frac{\theta }{2},\gamma }(\Omega ),\hfill \end{array}$$

(2)

in which $(·,·)$ indicates the inner product in $L^2(\Omega )$, $a(\varpi ,\upsilon )=C_{\theta }[2(\varpi ,\upsilon )-$$(D_{+}^{\frac{\theta }{2},\gamma }\varpi ,D_{-}^{\frac{\theta }{2},\gamma }\upsilon )-(D_{-}^{\frac{\theta }{2},\gamma }\varpi ,$$D_{+}^{\frac{\theta }{2},\gamma }\upsilon \left)\right]$, $C_{\theta }=-1/\left(2\cos (\theta \pi /2)\right)$, and ${\displaystyle \varpi (x,0)={\varpi }^0\left(x\right)(x\in \overline{\Omega })}$.

Noting that $a(·,·)$ is not coercive in $H_0^{\frac{\theta }{2},\gamma }(\Omega )$ as mentioned in [9], it is necessary to make Problem 1 into coerciveness so that it is uniquely solvable in $H_0^{\frac{\theta }{2},\gamma }(\Omega )$. Thereby, by setting $\varpi (x,t)=e^{\kappa t}U(x,t)$ ($\kappa >0$) we can change Problem 1 into the following system of equations with respect to the unknown function U.

$$\left\{\begin{array}{cc}{\partial }_{t}U(x,t)+{\partial }_{\left|x\right|}^{\theta ,\gamma ,\kappa }U(x,t)=\eta (x,t),\hfill & x\in \Omega ,t\in J,\hfill \\ U(0,t)=U(l,t)=0,\hfill & t\in \overline{J},\hfill \\ U(x,0)={\varpi }^0\left(x\right),\hfill & x\in \overline{\Omega },\hfill \end{array}\right.$$

(3)

where ${\partial }_{\left|x\right|}^{\theta ,\gamma ,\kappa }U(x,t)={\partial }_{\left|x\right|}^{\theta ,\gamma }U(x,t)+\kappa U(x,t)$ and $\eta (x,t)=e^{-\kappa t}\rho (x,t)$.

Thereupon, with integration by parts, the weak form for the system of Equations (3) is built as follows.

Problem 3.

Seek $U\in \mathbb{W}$ that satisfies

$$\begin{array}{}\mathrm{(4)}& (U_t,\upsilon )+a_{\kappa }(U,\upsilon )=(\eta ,\upsilon ),\forall \upsilon \in H_0^{\frac{\theta }{2},\gamma }(\Omega ),\mathrm{(5)}& U(x,0)={\varpi }^0\left(x\right),x\in \overline{\Omega },\end{array}$$

where $a_{\kappa }(U,\upsilon )=2C_{\theta }(U,\upsilon )$$+\kappa (U,\upsilon )$$-C_{\theta }(D_{+}^{\frac{\theta }{2},\gamma }U,D_{-}^{\frac{\theta }{2},\gamma }\upsilon )$$-C_{\theta }(D_{-}^{\frac{\theta }{2},\gamma }U,D_{+}^{\frac{\theta }{2},\gamma }\upsilon )$, and $C_{\theta }=-1/\left(2\cos \right(\theta \pi /2\left)\right)$.

The following result has proven in ([9], Lemma 3.3 and Theorems 4.3–4.5).

Theorem 1.

When ${\varpi }^0\left(x\right)\in H_0^{\frac{\theta }{2},\gamma }(\Omega )$ and $\rho (x,t)\in L^2(0,T;L^2(\Omega ))$, the bilinear functional $a_{\kappa }(·,·)$ is coercive and continuous in $H_0^{\frac{\theta }{2},\gamma }(\Omega )$ such that Problem 3 has a unique weak solution $U\in \mathbb{W}$.

Let ${\Im }_h$ be the regular subdivision onto $\overline{\Omega }=[0,l]$ such that $\overline{\Omega }=\{\cup K:K\in {\Im }_h\}$ and $h={\max }_{K\in {\Im }_h}\left\{\max \right|x_1-x_2|:x_1,x_2\in K\}$. The FE subspace is defined as:

$$\begin{array}{ll}{\mathbb{W}}_h& =\{{\upsilon }_h\in H_0^{\frac{\theta }{2},\gamma }(\Omega )\cap C(\Omega ):{\upsilon }_h{|}_K\in {\mathbb{P}}_m\left(K\right),K\in {\Im }_h\}\\ & =\mathrm{span}\left\{{\zeta }_j\left(x\right):1⩽j⩽M-1\right\},\end{array}$$

(6)

where ${\mathbb{P}}_m\left(K\right)$ consists of the mth-degree polynomials on $K\in {\Im }_h$, and ${\left\{{\zeta }_j\left(x\right)\right\}}_{j=1}^{M-1}$ is a set of orthonormal bases under the inner product in $L^2(\Omega )$ and may be obtained by the orthogonalization in ([32], Section 6.3 of Chapter 1).

For positive integer $N>0$, let $\tau =T/N$, $U_h^n$ be the CNFE solutions for Problem 3 at $t_n=n\tau$$(0⩽n⩽N)$, $\overline{\partial }{U}_h^n=(U_h^n-U_h^{n-1})/\tau$, ${\overline{U}}_h^n=(U_h^n+U_h^{n-1})/2$, and ${\eta }^{n-\frac{1}{2}}=\eta (x,t_{n-\frac{1}{2}})$. Thereupon, the CNFE model of Problem 3 is built as the following functional form.

Problem 4.

Seek $U_h^n\in {\mathbb{W}}_h$ ($1⩽n⩽N$) that satisfies

$$\begin{array}{}\mathrm{(7)}& (U_h^n,{\upsilon }_h)+\tau a_{\kappa }({\overline{U}}_h^n,{\upsilon }_h)=(U_h^{n-1},{\upsilon }_h)+\tau ({\eta }^{n-\frac{1}{2}},{\upsilon }_h),\forall {\upsilon }_h\in {\mathbb{W}}_h,1⩽n⩽N,\mathrm{(8)}& U_h^0\left(x\right)={\Pi }_h{\varpi }^0\left(x\right),\mathbf{x}\in \overline{\Omega },\end{array}$$

in which ${\Pi }_h:H_0^{\frac{\theta }{2},\gamma }(\Omega )\to {\mathbb{W}}_h$ is an interpolation operator.

The following result of the existence as well as the stability together with the error estimations for the CNFE solutions to Problem 4 was proved in ([9], Theorems 5.1 and 5.3 and Corollary 5.4).

Theorem 2.

Problem 4 has a unique set of unconditionally stabilized CNFE solutions ${\left\{U_h^n\right\}}_{n=1}^N\subset {\mathbb{W}}_h$. When the solution U to Problem 3 is sufficiently smooth, the CNFE solutions $U_h^n$ have the following error estimations

$$\begin{array}{c}{\displaystyle \parallel U\left(t_n\right)-U_h^n{\parallel }_0⩽\sigma \left({\tau }^2+h^{m+1}\right),n=1,2,\dots ,N.}\hfill \end{array}$$

(9)

where ${\parallel ·\parallel }_0$ is the norm in $L^2(\Omega )$ and σ represents a usually positive constant independent of τ and h, which may be unequal at different places. Furthermore, the CNFE solutions ${\varpi }_h^n=e^{\kappa t_n}{U}_h^n$$(n=1,2,...,N)$ are also uniquely existing and unconditionally stabilized, and they have the following error estimations

$$\begin{array}{c}{\displaystyle \parallel \varpi \left(t_n\right)-{\varpi }_h^n{\parallel }_0⩽\sigma \left({\tau }^2+h^{m+1}\right),n=1,2,\dots ,N.}\hfill \end{array}$$

(10)

2.2. Rewrite the CNFE Functional Form into Matrix Form

With the orthonormal bases ${\left\{{\zeta }_j\left(x\right)\right\}}_{j=1}^{M-1}$, the CNFE solutions to Problem 4 can be denoted by the following vector form:

$$U_h^n=\sum _{j=1}^{M-1}{z}_j^n{\zeta }_j\left(x\right)={\mathbf{{\rm Y}}}^n·\mathbf{\zeta },n=1,2,...,N,$$

where ${\mathbf{{\rm Y}}}^n={(z_1^n,z_2^n,...,z_{M-1}^n)}^T$ and $\mathbf{\zeta }={({\zeta }_1\left(x\right),{\zeta }_2\left(x\right),...,{\zeta }_{M-1}\left(x\right))}^T.$ Thereupon, Problem 4 can be rewritten as the following matrix form.

Problem 5.

Seek ${\mathbf{{\rm Y}}}^n\in R^{M-1}$, $U_h^n\in {\mathbb{W}}_h$, and ${\varpi }_h^n\in {\mathbb{W}}_h$$(1⩽n⩽N)$ that meet

$$\begin{array}{c}{\displaystyle \left(\mathit{I}+\frac{\tau }{2}\mathit{D}\right){\mathbf{{\rm Y}}}^n=\left(\mathit{I}-\frac{\tau }{2}\mathit{D}\right){\mathbf{{\rm Y}}}^{n-1}+\tau {\mathit{F}}^{n-\frac{1}{2}},n=1,2,\dots ,N,}\hfill \end{array}$$

(11)

$$\begin{array}{c}{\displaystyle U_h^n=\sum _{j=1}^{M-1}{z}_j^n{\zeta }_j\left(x\right)={\mathbf{{\rm Y}}}^n·\mathbf{\zeta },{\varpi }_h^n=e^{\kappa t_n}{U}_h^n,n=1,2,\dots ,N,}\hfill \end{array}$$

(12)

where the matrix $\mathit{D}={\left(a_{\kappa }({\zeta }_i,{\zeta }_j)\right)}_{(M-1)\times (M-1)}$ is symmetrical positive definite owing to the coerciveness of $a_{\kappa }(·,·)$, $\mathit{I}$ is an $(M-1)\times (M-1)$ identity matrix, ${\mathbf{{\rm Y}}}^0=({\varpi }^0\left(x_1\right),{\varpi }^0\left(x_2\right),\dots ,$${\varpi }^0\left(x_{M-1}\right){)}^T$, $x_j$’s are the inner nodes in ${\Im }_h$, and ${\mathit{F}}^{n-\frac{1}{2}}=(({\eta }^{n-\frac{1}{2}},{\zeta }_1),$$({\eta }^{n-\frac{1}{2}},{\zeta }_2),$$\dots ,({\eta }^{n-\frac{1}{2}},{\zeta }_{M-1}){)}^T$.

Remark 1.

When the FEs in ${\Im }_h$ need to be sufficiently refined, there will be many unknowns in Problem 5, resulting in that the rounding errors are accumulated rapidly in the calculation process and the CNFE solutions appear to have a large deviation. Hence, it is necessary to adopt the POD technique to lower the dimension of unknown solution coefficient vectors ${\mathbf{{\rm Y}}}^n$ in Problem 5.

3. The RDRCNFE Method for STFDE

3.1. Structure of POD Basic Vectors

We first seek the initial L solution coefficient vectors ${\mathbf{{\rm Y}}}^n$ ($n=1,2,\dots ,L$) of Problem 5 to make up the matrix $\mathit{E}=\left({\mathbf{{\rm Y}}}^1,{\mathbf{{\rm Y}}}^2,...,\right.$${\left.{\mathbf{{\rm Y}}}^L\right)}_{(M-1)\times L}$ and calculate the positive eigenvalues ${\chi }_j>0$$(j=1,$$2,$$\cdots ,\kappa =$ rank$\left(\mathit{E}\right))$ (degressively sequenced) and the relative orthonormal eigenvectors $\tilde{\Psi }=({\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,\cdots ,$${\mathbf{\varphi }}_{\kappa })\in R^{M\times \kappa }$ of $\mathit{E}{\mathit{E}}^T$. Thus, a set of POD bases $\Psi =({\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,\cdots ,$${\mathbf{\varphi }}_d)$, consisting of the initial d vectors in $\tilde{\Psi }$, has the following property (see [4]):

$$\begin{array}{c}{\displaystyle \parallel \mathit{E}-\Psi {\Psi }^T{\mathit{E}\parallel }_{2,2}=\sqrt{{\chi }_{d+1}},d⩽\kappa ,}\hfill \end{array}$$

(13)

where ${\parallel \mathit{E}\parallel }_{2,2}={\sup }_{\mathbf{y}\ne 0}\parallel \mathit{E}\mathbf{y}\parallel /\parallel \mathbf{y}\parallel$ and $\parallel \mathbf{y}\parallel$ are the Euclidian norm of vector $\mathbf{y}$. Thereupon, we obtain

$$\begin{array}{c}{\displaystyle \parallel {\mathbf{{\rm Y}}}^n-\Psi {\Psi }^T{\mathbf{{\rm Y}}}^n\parallel =\parallel (\mathit{E}-\Psi {\Psi }^T\mathit{E}){\mathit{e}}^n\parallel ⩽\parallel \mathit{E}-\Psi {\Psi }^T{\mathit{E}\parallel }_{2,2}\parallel {\mathit{e}}^n\parallel ⩽\sqrt{{\chi }_{d+1}},1⩽n⩽L,}\hfill \end{array}$$

(14)

in which ${\mathit{e}}^n$$(1⩽n⩽L)$ represents the identity vectors with the nth component 1.

Remark 2.

Owing to $(M-1)\gg L$, but both ${\mathit{E}}^T\mathit{E}$ and $\mathit{E}{\mathit{E}}^T$ have identical positive eigenvalues ${\chi }_j$$(1⩽j⩽r)$, we can firstly find that the main d eigenvectors ${\mathbf{\phi }}_j$$(1⩽j⩽d)$ correspond to most of the main eigenvalues ${\chi }_j$$(1⩽j⩽d)$ of ${\mathit{E}}^T\mathit{E}$, we then can lightly obtain most of the main d eigenvectors ${\mathbf{\varphi }}_j=\mathit{E}{\mathbf{\phi }}_j/\sqrt{{\chi }_j}$$(1⩽j⩽d)$ of $\mathit{E}{\mathit{E}}^T$ to make up a set of POD basic vectors ${\Psi }_i=({\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,\cdots ,$${\mathbf{\varphi }}_d)$ ($d⩽\kappa$).

3.2. Construction of RDRCNFE Method

If we assume that ${\mathit{b}}_d^n={(b_1^n,b_2^n,\dots ,b_d^n)}^T$, ${\mathbf{{\rm Y}}}_d^n$$=(z_{d1}^n,z_{d2}^n,...,$$z_{d(M-1)}^n{)}^T=\Psi {\mathit{b}}_d^n=\Psi {\Psi }^T{\mathbf{{\rm Y}}}^n$, and $U_d^n=\mathbf{\zeta }·{\mathbf{{\rm Y}}}_d^n$, we immediately gain the initial L RDRCNFE solutions $U_d^n=\mathbf{\zeta }·{\mathbf{{\rm Y}}}_d^n=\mathbf{\zeta }·(\Psi {\Psi }^T{\mathbf{{\rm Y}}}^n)$ ($1⩽n⩽L$). If the unknown solution vectors ${\mathbf{{\rm Y}}}^n$$(L+1⩽n⩽N)$ in Problem 5 are replaced with ${\mathbf{{\rm Y}}}_d^n=\Psi {\mathit{b}}_d^n$$(L+1⩽n⩽N)$, by the positive definiteness and invertibility of matrix $(\mathit{I}+\tau \mathit{D}/2)$, we could build the following RDRCNFE method.

Problem 6.

Seek ${\mathit{b}}_d^n\in {\mathbb{R}}^d$, $U_d^n\in {\mathbb{W}}_h$, and ${\varpi }_d^n\in {\mathbb{W}}_h$$(n=1,2,\dots ,N)$ that meet

$$\begin{array}{l}{\mathit{b}}_d^n={\Psi }^T{\mathbf{{\rm Y}}}^n,1⩽n⩽L;\\ {\mathit{b}}_d^n={\Psi }^T{\left(\mathit{I}+\frac{\tau }{2}\mathit{D}\right)}^{-1}\left(\mathit{I}-\frac{\tau }{2}\mathit{D}\right)\Psi {\mathit{b}}_d^{n-1}\end{array}$$

(15)

$$+\tau {\Psi }^T{\left(\mathit{I}+\frac{\tau }{2}\mathit{D}\right)}^{-1}{\mathit{F}}^{n-\frac{1}{2}},L+1⩽n⩽N,$$

(16)

$$U_d^n=\sum _{j=1}^{M-1}{z}_{dj}^n{\zeta }_j\left(x\right)=\mathbf{\zeta }·(\Psi {\mathit{b}}_d^n),{\varpi }_d^n=e^{\kappa t_n}{U}_d^n,1⩽n⩽N,$$

(17)

where ${\mathbf{{\rm Y}}}^n$$(1⩽n⩽L)$ are the initial L solution vectors to Problem 5, and the matrix $\mathit{D}$ and vectors ${\mathit{F}}^{n-\frac{1}{2}}$ are given in Problem 5.

Remark 3.

It is obvious that Problem 6 has a unique set of RDRCNFE solutions ${\left\{U_d^n\right\}}_{n=1}^N\subset {\mathbb{W}}_h$ so that the RDRCNFE solutions ${\left\{{\varpi }_d^n\right\}}_{n=1}^N={\left\{e^{\kappa t_n}{U}_d^n\right\}}_{n=1}^N\subset {\mathbb{W}}_h$ of Problem 1 are also uniquely existing. Specially, at each time node, Problem 5 has $(M-1)$ unknowns, while Problem 6 has only d unknowns $(d\ll M-1)$, but it has the same basis functions ${\left\{{\zeta }_i\left(x\right)\right\}}_{i=1}^M$ as Problem 5. Namely, although the unknowns of Problem 6 are greatly reduced, it keeps the basis functions unchanged so as to maintain the accuracy unchanged, too. Hence, Problem 6 is obviously superior to Problem 5. This signifies that Problem 6 can not only immensely salvage CPU runtime and slow down the rounding error amassing but also raise the accuracy for numerical solutions in the practical calculations.

3.3. Stability and Error Estimations of the RDRCNFE Solutions

The theoretical analysis for the stability together with the errors of the RDRCNFE solutions requires the following matrix properties (see [33], Theorems 1.4.1 and 1.4.2).

Lemma 1.

The symmetrical positive definite matrix $\mathit{D}$ in Problem 5 possesses the following properties:

$$\parallel {\left(\mathit{I}+0.5\tau \mathit{D}\right)}^{-1}\left(\mathit{I}-0.5\tau \mathit{D}\right){\parallel }_{2,2}⩽1;{\parallel {\left(\mathit{I}+0.5\tau \mathit{D}\right)}^{-1}\parallel }_{2,2}⩽1.$$

For the RDRCNFE solutions to Problem 6, we have the following result of stability together with error estimations.

Theorem 3.

Under the same conditions as Theorem 2, the RDRCNFE solutions $U_d^n$$(n=1,2,\dots ,N)$ to Problem 6 are unconditionally stabilized and have the following error estimations

$$\begin{array}{c}{\displaystyle \parallel U\left(t_n\right)-U_d^n{\parallel }_0⩽\sigma \left({\tau }^2+h^{m+1}+\sqrt{{\chi }_{d+1}}\right),1⩽n⩽N,}\hfill \end{array}$$

(18)

where $U\left(t_n\right)$ are the state of solutions of Problem 3 at $t_n=n\tau$$(1⩽n⩽N)$. Furthermore, the associated RDRCNFE solutions ${\varpi }_d^n=e^{\kappa t_n}{U}_d^n$$(1⩽n⩽N)$ are also unconditionally stabilized and have the following error estimations

$$\begin{array}{c}{\displaystyle \parallel \varpi \left(t_n\right)-{\varpi }_d^n{\parallel }_0⩽\sigma \left({\tau }^2+h^{m+1}+\sqrt{{\chi }_{d+1}}\right),1⩽n⩽N.}\hfill \end{array}$$

(19)

Proof. 

(i) Prove the stability of solutions of Problem 6.

While $n=1,2,\dots ,L$, based on the orthonormality of POD bases $\Psi$, we obtain

$$\begin{array}{c}{\displaystyle \parallel U_d^n{\parallel }_0=\parallel {\mathbf{{\rm Y}}}_d^n{·\mathbf{\zeta }\parallel }_0=\parallel \Psi {\Psi }^T{\mathbf{{\rm Y}}}^n{·\mathbf{\zeta }\parallel }_0⩽\sigma {\parallel U_h^n\parallel }_0.}\hfill \end{array}$$

(20)

Thus, with the unconditionally stability of $U_h^n$ given by Theorem 2, we immediately assert that ${\left\{U_d^n\right\}}_{n=1}^L$ is unconditionally stabilized.

While $n=L+1,L+2,\dots ,N$, by ${\mathbf{{\rm Y}}}_d^n=\Psi {\mathit{b}}_d^n$, we could, respectively, revert (16) and (17) into

$$\begin{array}{}\mathrm{(21)}& {\mathbf{{\rm Y}}}_d^n={\left(\mathit{I}+\frac{\tau }{2}\mathit{D}\right)}^{-1}\left(\mathit{I}-\frac{\tau }{2}\mathit{D}\right){\mathbf{{\rm Y}}}_d^{n-1}+\tau {\left(\mathit{I}+\frac{\tau }{2}\mathit{D}\right)}^{-1}{\mathit{F}}^{n-\frac{1}{2}},n=L+1,L+2,\dots ,N,\mathrm{(22)}& U_d^n=\mathbf{\zeta }·{\mathbf{{\rm Y}}}_d^n,n=L+1,L+2,\dots ,N.\end{array}$$

Using Lemma 1, from (21) we obtain

$$\begin{array}{c}{\displaystyle \parallel {\mathbf{{\rm Y}}}_d^n\parallel ⩽\parallel {\left(\mathit{I}+0.5\tau \mathit{D}\right)}^{-1}\left(\mathit{I}-0.5\tau \mathit{D}\right){\parallel }_{2,2}\parallel {\mathbf{{\rm Y}}}_d^{n-1}\parallel +\tau \parallel {\left(\mathit{I}+0.5\tau \mathit{D}\right)}^{-1}{\parallel }_{2,2}\parallel {\mathit{F}}^{n-\frac{1}{2}}\parallel }\hfill \\ {\displaystyle ⩽\parallel {\mathbf{{\rm Y}}}_d^{n-1}\parallel +\tau \parallel {\mathit{F}}^{n-\frac{1}{2}}\parallel ,n=L+1,L+2,\dots ,N.}\hfill \end{array}$$

(23)

Using (11) and Lemma 1, we obtain

$$\begin{array}{c}{\displaystyle \parallel {\mathbf{{\rm Y}}}^n\parallel ⩽\parallel {\left(\mathit{I}+0.5\tau \mathit{D}\right)}^{-1}\left(\mathit{I}-0.5\tau \mathit{D}\right){\mathbf{{\rm Y}}}^{n-1}\parallel +\tau \parallel {\left(\mathit{I}+0.5\tau \mathit{D}\right)}^{-1}{\mathit{F}}^{n-\frac{1}{2}}\parallel }\hfill \\ {\displaystyle ⩽\parallel {\mathbf{{\rm Y}}}^{n-1}\parallel +\tau \parallel {\mathit{F}}^{n-\frac{1}{2}}\parallel }\hfill \\ {\displaystyle ⩽\parallel {\mathbf{{\rm Y}}}^0\parallel +\tau \sum _{i=1}^n\parallel {\mathit{F}}^{i-\frac{1}{2}}\parallel }\hfill \\ {\displaystyle ⩽\sigma ,n=1,2,...,N.}\hfill \end{array}$$

(24)

Thus, summing for (23) from $L+1$ unto n, noting that $\parallel {\mathbf{{\rm Y}}}_d^L\parallel =\parallel \Psi {\Psi }^T{\mathbf{{\rm Y}}}^L\parallel$$⩽\sigma \parallel {\mathbf{{\rm Y}}}^L\parallel$, by (24) we obtain

$$\begin{array}{c}{\displaystyle \parallel {\mathbf{{\rm Y}}}_d^n\parallel ⩽\parallel {\mathbf{{\rm Y}}}_d^L\parallel +\tau \sum _{i=L+1}^n\parallel {\mathit{F}}^{i-\frac{1}{2}}\parallel }\hfill \\ {\displaystyle ⩽\sigma \parallel {\mathbf{{\rm Y}}}^L\parallel +\tau \sum _{i=L+1}^n\parallel {\mathit{F}}^{i-\frac{1}{2}}\parallel }\hfill \\ {\displaystyle ⩽\sigma ,n=L+1,L+2,\dots ,N.}\hfill \end{array}$$

(25)

Thus, noting that $U_d^n=\mathbf{\zeta }·{\mathbf{{\rm Y}}}_d^n$, ${\parallel \mathbf{\zeta }\parallel }_0⩽\sigma$, and ${\parallel \mathbf{y}\parallel }_{\infty }⩽\parallel \mathbf{y}\parallel$, we obtain

$$\begin{array}{c}{\displaystyle \parallel U_d^n{\parallel }_0⩽\parallel {\mathbf{{\rm Y}}}_d^n{·\mathbf{\zeta }\parallel }_0⩽\parallel {\mathbf{{\rm Y}}}_d^n{\parallel }_{\infty }{\parallel \mathbf{\zeta }\parallel }_0⩽\parallel {\mathbf{{\rm Y}}}_d^n{\parallel ·\parallel \mathbf{\zeta }\parallel }_0⩽\sigma ,n=L+1,L+2,...,N,}\hfill \end{array}$$

(26)

which implies that ${\left\{U_d^n\right\}}_{n=L+1}^N$ is also unconditionally stable. Hence, both ${\left\{U_d^n\right\}}_{n=1}^N$ and ${\left\{{\varpi }_d^n\right\}}_{n=1}^N={\left\{e^{\kappa t_n}{U}_d^n\right\}}_{n=1}^N$ are unconditionally stabilized.

(ii) Estimate the errors of ROECNFSE solutions.

While $1⩽n⩽L$, noting that $U_h^n=\mathbf{\zeta }·{\mathbf{{\rm Y}}}^n$, ${\parallel \mathbf{\zeta }\parallel }_0⩽\sigma$, and ${\parallel \mathbf{y}\parallel }_{\infty }⩽\parallel \mathbf{y}\parallel$ ($\forall \mathbf{y}\in {\mathbb{R}}^{(M-1)}$), by (14) we obtain

$$\begin{array}{c}{\displaystyle \parallel U_h^n-U_d^n{\parallel }_0⩽{\parallel \mathbf{\zeta }\parallel }_0\parallel {\mathbf{{\rm Y}}}^n-{\mathbf{{\rm Y}}}_d^n{\parallel }_{\infty }⩽\sigma \parallel {\mathbf{{\rm Y}}}^n-\Psi {\Psi }^T{\mathbf{{\rm Y}}}^n\parallel ⩽\sigma \sqrt{{\chi }_{d+1}}.}\hfill \end{array}$$

(27)

While $n=L+1,L+2,\dots ,N$, by (11) and (21), using Lemma 1, we obtain

$$\begin{array}{c}{\displaystyle \parallel {\mathbf{{\rm Y}}}^n-{\mathbf{{\rm Y}}}_d^n\parallel =\parallel {\left(\mathit{I}+0.5\tau \mathit{D}\right)}^{-1}\left(\mathit{I}-0.5\tau \mathit{D}\right)({\mathbf{{\rm Y}}}^{n-1}-{\mathbf{{\rm Y}}}_d^{n-1})\parallel }\hfill \\ {\displaystyle ⩽\parallel {\left(\mathit{I}+0.5\tau \mathit{D}\right)}^{-1}\left(\mathit{I}-0.5\tau \mathit{D}\right){\parallel }_{2,2}\parallel {\mathbf{{\rm Y}}}^{n-1}-{\mathbf{{\rm Y}}}_d^{n-1}\parallel }\hfill \\ {\displaystyle ⩽\parallel {\mathbf{{\rm Y}}}^{n-1}-{\mathbf{{\rm Y}}}_d^{n-1}\parallel .}\hfill \end{array}$$

(28)

Thus, from (28) and (27), we obtain

$$\begin{array}{c}{\displaystyle \parallel {\mathbf{{\rm Y}}}^n-{\mathbf{{\rm Y}}}_d^n\parallel ⩽\parallel {\mathbf{{\rm Y}}}^L-{\mathbf{{\rm Y}}}_d^L\parallel }\hfill \\ {\displaystyle ⩽\sigma \sqrt{{\chi }_{d+1}},n=L+1,L=2,\dots ,N.}\hfill \end{array}$$

(29)

Therefore, we have

$$\begin{array}{c}{\displaystyle \parallel U_h^n-U_d^n{\parallel }_0={\parallel \mathbf{\zeta }·({\mathbf{{\rm Y}}}^n-{\mathbf{{\rm Y}}}_d^n)\parallel }_0}\hfill \\ {\displaystyle ⩽{\parallel \mathbf{\zeta }\parallel }_0{\parallel {\mathbf{{\rm Y}}}^n-{\mathbf{{\rm Y}}}_d^n\parallel }_{\infty }}\hfill \\ {\displaystyle ⩽\sigma \parallel {\mathbf{{\rm Y}}}^n-{\mathbf{{\rm Y}}}_d^n\parallel }\hfill \\ {\displaystyle ⩽\sigma \sqrt{{\chi }_{d+1}},n=L+1,L=2,\dots ,N.}\hfill \end{array}$$

(30)

Combining (27) and (30) together with Theorem 2, we obtain (18). By $\varpi (x,t_n)=e^{\kappa t_n}U(x,t_n)$, ${\varpi }_h(x,t_h)=e^{\kappa t_n}{U}_h^n$, and ${\varpi }_d^n=e^{\kappa t_n}{U}_d^n$ ($n=1,2,\dots ,N$), from (18) and Theorem 2, we immediately obtain (19). Theorem 3 is proved. □

Remark 4.

Although the error estimates in Theorem 3 have one more term $\sqrt{{\chi }_{d+1}}$ than those in Theorem 2 due to the dimensionality reduction for the CNFE method, it can be used as a criterion for selecting the POD basis vectors. As long as the number d of POD bases meets $\sqrt{{\chi }_{d+1}}⩽{\tau }^2+h^{m+1}$, there is little effect on the total errors. A lot of numerical simulations (see [4,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]) have shown that the series $\left\{{\chi }_j\right\}$ of eigenvalues would rapidly decrease to 0. In a general way, while d = 5~7, $\sqrt{{\chi }_{d+1}}$ is already extremely small. Especially, if the RDRCNFE solution $U_d^{n_0+1}$ obtained by Problem 6 at some time node $t_{n_0+1}$ cannot satisfy the accuracy requirement, but $U_d^n$ at the time nodes $t_n$$(n⩽n_0)$ can still satisfy, we may choose a new set of solution vectors $({\mathbf{{\rm Y}}}^{n_0+1-L},{\mathbf{{\rm Y}}}^{n_0+2-L},\dots ,{\mathbf{{\rm Y}}}^{n_0-1},{\mathbf{{\rm Y}}}^{n_0})$ to construct a new set of POD basis vectors and to build the new RDRCNFE method, and then, we can calculate out the RDRCNFE solutions that satisfy accuracy requirements. In this way, we can calculate out the RDRCNFE solutions at an arbitrary time node. This is incomparable to the traditional CNFE method.

4. Some Numerical Simulations

Here, the correctness of theory results and the superiority of the RDRCNFE method are verified by the means of some numerical simulations. Problem 1 has an analytical solution, but it usually has no analytic solution when the source term and initial functions are complicated.

If the initial function ${\varpi }^0\left(x\right)=2{\gamma }^{6-\theta }{x}^3{(x-1)}^3\Gamma (-\theta )\cos (\theta \pi /2)$ and the source function $\rho (x,t)$ is denoted by

$$\begin{array}{ccc}& & {\displaystyle \rho (x,t)=2{\gamma }^{6-\theta }\cos (\pi \theta /2)\Gamma (-\theta )e^{-t}{x}^3{(1-x)}^3}\hfill \\ & & {\displaystyle +e^{-t}\left[3{\gamma }^5{(x-1)}^2{x}^2(2x-1)\nu (1-\theta ,x\gamma )\right.}\hfill \\ & & {\displaystyle +\gamma \left({\gamma }^5{x}^3(3x^2-x^3-3x+1)\nu (-\theta ,x\gamma )\right.}\hfill \\ & & {\displaystyle +3{\gamma }^{3}x(10x^2-5x^3-6x+3)\nu (2-\theta ,\gamma -\gamma x)}\hfill \\ & & {\displaystyle -{\gamma }^5{(x-1)}^3{x}^3\nu (-\theta ,\gamma -\gamma x)+{\gamma }^2(20x^3-30x^2+12x-1)\nu (3-\theta ,\gamma x)}\hfill \\ & & {\displaystyle +{\gamma }^2(1-12x+30x^2-20x^3)\nu (3-\theta ,\gamma -\gamma x)}\hfill \\ & & {\displaystyle -3{\gamma }^4{(x-1)}^2(2x-1)x^2\nu (1-\theta ,\gamma -\gamma x)}\hfill \\ & & {\displaystyle -3\gamma (5x^2-5x+1)(\nu (4-\theta ,\gamma x)+\nu (4-\theta ,\gamma -\gamma x))}\hfill \\ & & {\displaystyle +3(x-1)\left(\nu \right(5-\theta ,\gamma x)-\nu (5-\theta ,\gamma -\gamma x\left)\right)}\hfill \\ & & {\displaystyle \left.-3{\gamma }^{3}x(x-1)(5x^2-5x+1)\nu (2-\theta ,\gamma x)\right)}\hfill \\ & & {\displaystyle +2\Gamma (2-\theta )\left({\theta }^4-14{\theta }^3+71{\theta }^2-154\theta +120\right.}\hfill \\ & & {\displaystyle \left.+3(3-\theta )(2-\theta ){\gamma }^2(5x^2-5x+1)+3{\gamma }^{4}x(x-1)(5x^2-5x+1)\right)}\hfill \\ & & {\displaystyle \left.-\nu (6-\theta ,\gamma x)-6\nu (6-\theta ,\gamma -\gamma x)\right],}\hfill \end{array}$$

then Problem 1 has an analytic solution

$$\begin{array}{ccc}& & {\displaystyle \varpi (x,t)=2{\gamma }^{6-\theta }{x}^3{(x-1)}^3{e}^{-t}\cos (\theta \pi /2)\Gamma (-\theta ).}\hfill \end{array}$$

When $\Omega =J=(0,1)$, $\tau =h=1/1000$, and ${\mathbb{P}}_1\left(K\right)$ consists of linear polynomials (i.e., $m=1$), we firstly calculate out the 20 initial CNFE solutions ${\mathbf{{\rm Y}}}^n$ ($n=1,2,...,20$) under various cases of $\theta =1.1$, $1.5$, and $1.9$, and $\gamma =0.5$, $1.0$, and $2.0$ by Problem 5 and make up the matrix $\mathit{E}=({\mathbf{{\rm Y}}}^1,{\mathbf{{\rm Y}}}^2,...,{\mathbf{{\rm Y}}}^{20})$, respectively. We then calculate out the eigenvalues ${\chi }_i$ arrayed degressively and eigenvectors ${\mathbf{\phi }}_i$$(i=1,2,...,20)$ of the matrix ${\mathit{E}}^T\mathit{E}$ corresponding to the various cases of $\theta =1.1$, $1.5$, and $1.9$, and $\gamma =0.5$, $1.0$, and $2.0$, respectively. By reckoning, we find that $\sqrt{{\chi }_7}⩽3\times {10}^{-6}$. Thus, we just have to take the foremost six eigenvectors ${\mathbf{\varphi }}_i$$(i=1,2,...,6)$ and produce a set of POD basis vectors $\Psi =({\mathbf{\phi }}_1,{\mathbf{\phi }}_2,...,{\mathbf{\phi }}_6)$ by ${\mathbf{\phi }}_i=\mathit{E}{\mathbf{\varphi }}_i/\sqrt{{\chi }_i}$$(1⩽i⩽6)$. Finally, by using Matlab on a laptop (Think-Pad E530) to solve Problem 6, we calculate out the RDRCNFE solutions ${\varpi }_d^n=e^{t_n}{U}_d^n$ ($n=1000$, i.e., at $t=1$) under the various cases of $\theta =1.1$, $1.5$, and $1.9$, and $\gamma =0.5$, $1.0$, and $2.0$, respectively.

In order to explain the superiority of the RDRCNFE method by using the same Matlab on the same laptop to solve the CNFE method, we also calculate out the corresponding CNFE solutions ${\varpi }_h^n$ (when $n=1000$, i.e., at $t=1$) under the various cases of $\theta =1.1$, $1.5$, and $1.9$, and $\gamma =0.5$, $1.0$, and $2.0$, respectively, and record the $L^2$ errors $\parallel \varpi \left(t_n\right)-{\varpi }_h^n{\parallel }_0$ between the analytic solutions $\varpi \left(t_n\right)$ and the CNFE solution ${\varpi }_h^n$, the errors $\parallel \varpi \left(t_n\right)-{\varpi }_d^n{\parallel }_0$ between the analytic solutions $\varpi \left(t_n\right)$ and the RDRCNFE solution ${\varpi }_d^n$, and the CPU runtime for the CNFE method and the RDRCNFE method when $n=1000$, i.e., $t=1$ under the various cases of $\theta =1.1$, $1.5$, and $1.9$, and $\gamma =0.5$, $1.0$, and $2.0$, as listed in

Table 1. CPU runtime and errors between the analytic and the CNFE as well as RDRCNFE solutions when $h=\tau ={10}^{-3}$.

		CNFE Method		RDRCNFE Method
θ	γ	$\parallel \mathit{\varpi }\left({\mathit{t}}_{\mathit{n}}\right)-{\mathit{\varpi }}_{\mathit{h}}^{\mathit{n}}{\parallel }_0$	CPU Runtime	$\parallel \mathit{\varpi }\left({\mathit{t}}_{\mathit{n}}\right)-{\mathit{\varpi }}_{\mathit{d}}^{\mathit{n}}{\parallel }_0$	CPU Runtime
1.1	0.5	1.010356 × 10−6	43.568 S	4.150523 × 10−6	1.623 S
	1.0	1.012083 × 10−6	43.865 S	4.250732 × 10−6	1.665 S
	2.0	1.125338 × 10−6	43.914 S	5.071732 × 10−6	1.673 S
1.5	0.5	1.315376 × 10−6	43.931 S	4.352762 × 10−6	1.692 S
	1.0	1.414376 × 10−6	43.982 S	4.651718 × 10−6	1.713 S
	2.0	1.534283 × 10−6	44.173 S	5.052123 × 10−6	1.721 S
1.9	0.5	1.541232 × 10−6	43.842 S	4.356431 × 10−6	1.676 S
	1.0	1.562183 × 10−6	43.874 S	4.672762 × 10−6	1.813 S
	2.0	1.612386 × 10−6	44.187 S	5.131753 × 10−6	1.925 S

.

Table 1 explains that the CPU runtime of the RDRCNFE method is about 26 times that of the CNFE method. Thus, the RDRCNFE method can greatly lessen the CPU runtime. Especially, at each time node, the RDRCNFE method only includes six unknowns, but the CNFE method has a thousand unknowns. When the CNFE method is applied to large-scale numerical calculations in the real world, it has more than millions of unknowns. Hence, the RDRCNFE method can not only greatly lessen the CPU runtime and decrease the rounding error amassing but also raise the calculating accuracy so that the RDRCNFE method is far superior to the CNFE method.

When $\sqrt{{\chi }_7}=O\left({10}^{-6}\right)$ and $\tau =h=1/1000$, by Theorems 2 and 3, we can obtain that the theory errors between the analytic solutions $\varpi \left(t_n\right)$ and the CNFE solution ${\varpi }_h^n$ together with the theory errors between the analytic solutions $\varpi \left(t_n\right)$ and the RDRCNFE solutions ${\varpi }_d^n$ are about $O\left({10}^{-6}\right)$, but the numerical calculating errors in Table 1 are also $O\left({10}^{-6}\right)$. It is shown that the numerical results are in accord with the theoretical results.

To obtain the intuition, we exhibit the RDRCNFE solutions at $t=0.25$, $0.50$, $0.75$, and $1.00$ for $\gamma =1.0$ and $\theta =1.1$, $1.5$, and $1.9$ in Figure 1a, Figure 2a, Figure 3a, respectively. To compare with the RDRCNFE solutions, we also exhibit the CNFE solutions at $t=0.25$, $0.50$, $0.75$, and $1.00$ for $\gamma =1.0$ and $\theta =1.1$, $1.5$, and $1.9$ in Figure 1b, Figure 2b, Figure 3b, respectively. By comparison, we find that each pair in Figure 1, Figure 2 and Figure 3 is highly similar, which implies that the RDRCNFE solutions have the same accuracy as the the CNFE solutions and the RDRCNFE method is very effective for settling STFDE (i.e., Problem 1) even if it only employs six POD basis vectors.

5. Conclusions and Discussions

Herein, we have dealt with the reduced dimension for the unknown solution coefficient vectors to CNFE method of STFDE. We have made use of the POD technique to build the RDRCNFE method for STFDE, adopted the matrix analysis to analyze the stability and convergence of RDRCNFE solutions, and used some numerical simulations to confirm the correctness of theoretical results and the superiority of the RDRCNFE method. In particular, the RDRCNFE method is proposed for the first time and is absolutely different from the existed reduced-order methods such as the reduced-order methods in [4,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]. Therefore, the RDRCNFE method is a new development over the existing numerical methods.

Although we have only dealt with the reduced dimension for the unknown solution coefficient vectors to the CNFE method of STFDE, the method can be directly generalized to the two-dimensional (2D) or three-dimensional (3D) cases and even the more intricate real-world engineering numerical simulations. In fact, when the FE models for any 2D and 3D steady PDEs including the 2D and 3D STFDE are used to perform the numerical simulations on the computer, they need to express as the following matrix form

$${\mathit{U}}^{n+1}=\mathit{B}({\mathit{U}}^{n+1},{\mathit{U}}^n,{\mathit{U}}^{n-1}),n=1,2,\cdots ,$$

(31)

where ${\mathit{U}}^n$ are very high-dimension vectors, which can reach millions or even tens of millions of dimensions in the actual engineering computation. Thus, we may first compute the first L steps solutions ${\left\{{\mathit{U}}^n\right\}}_{n=1}^L$ (empirical value $L=20$), which are also formed from the observation values of experiments on all grid points at the L moments, to form the snapshot $\mathit{E}=({\mathit{U}}^1,{\mathit{U}}^2,\cdots ,{\mathit{U}}^L)$. Thereupon, we may employ the discrete POD method in Section 3.1 to find the POD basis vectors $\Psi =({\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,\cdots ,{\mathbf{\varphi }}_d)$, which are the eigenvectors ${\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,\cdots ,{\mathbf{\varphi }}_d$ corresponding to the main eigenvalues ${\lambda }_1⩾{\lambda }_2⩾\cdots {\lambda }_d>0$ of $\mathit{E}{\mathit{E}}^T$.

Let ${\mathit{U}}_d^n=\Psi {\mathit{b}}^n$ and ${\mathit{b}}^n={(b_1^n,b_2^n,\cdots ,b_d^n)}^T$ be the d-dimensional unknown vectors. Substituting ${\mathit{U}}_d^n=\Psi {\mathit{b}}^n$ into ${\mathit{U}}^n$ in the above large-scale linear or nonlinear algebra Equation (31), we obtain the following reduced-dimension system of equations that only includes d unknowns:

$$\{\begin{array}{l}{\mathit{b}}_d^n={\Psi }^T{\mathit{U}}^n,n=1,2,\cdots ,L;\\ \Psi {\mathit{b}}^{n+1}=\mathit{B}(\Psi {\mathit{b}}^{n+1},\Psi {\mathit{b}}^n,\Psi {\mathit{b}}^{n-1}),n=L,L+1,\cdots ,\\ {\mathit{U}}_d^n=\Psi {\mathit{b}}^n,n=1,2,\cdots ,L,L+1,L+2,\cdots ,\end{array}$$

(32)

where ${\mathit{U}}^n$$(n=1,2,\cdots ,L)$ are the known solution vectors for the system of Equation (31) or observation values of experiments on all grid points at the L moments. If the reduced-dimension equations are linear, they can be directly solved iteratively; if the reduced-dimension equations are nonlinear, they can be solved by the Newton method or other methods of nonlinear algebra equations. As a consequence, the RDRCNFE method possesses very extensive applying foreground.