Block-Centered Finite-Difference Methods for Time-Fractional Fourth-Order Parabolic Equations

Abstract

The block-centered finite-difference method has many advantages, and the time-fractional fourth-order equation is widely used in physics and engineering science. In this paper, we consider variable-coefficient fourth-order parabolic equations of fractional-order time derivatives with Neumann boundary conditions. The fractional-order time derivatives are approximated by $L1$ interpolation. We propose the block-centered finite-difference scheme for fourth-order parabolic equations with fractional-order time derivatives. We prove the stability of the block-centered finite-difference scheme and the second-order convergence of the discrete $L2$ norms of the approximate solution and its derivatives of every order. Numerical examples are provided to verify the effectiveness of the block-centered finite-difference scheme.

1. Introduction

The block-centered finite-difference method was first applied to the simulation of oil reservoirs [1]. Russell and Wheeler [1] proved that the block-centered finite-difference method is equivalent to the mixed finite-element method with a special numerical quadrature formula. Based on this equivalence, it is easier to discuss the stability and convergence of the block-centered finite-difference method. In addition, the block-centered finite-difference method can simultaneously approximate the exact solution of the original problem and its derivatives, thus preserving the local conservation of the problem. Moreover, for problems with Neumann boundary conditions, the numerical solution of nodes near the boundary does not need to be considered separately. On this basis, Weiser and Wheeler [2] introduced the block-centered finite-difference method for linear self-adjoint and non-self-adjoint elliptic and parabolic problems with Neumann boundary conditions in a rectangular area. They proved that the errors of the discrete L2 norms for the solution and the first derivative are both of the second order. Refs. [3,4,5] considered the block-centered finite-difference method for the nonlinear Darcy–Forchheimer equation. The block-centered finite-difference method on non-uniform grids has been discussed in [6,7,8,9,10,11]. In [12], Ren and Zhang studied the Crank–Nicolson block-centered difference method for solving linear parabolic equations in bounded domains. Li and Rui [13] introduced and analyzed the block-centered finite-difference method for distributed-order time-fractional diffusion-wave equations with Neumann boundary conditions. In addition, [14,15,16,17,18] discussed the two-grid and parallel block-centered finite-difference schemes for parabolic equations and diffusion equations with fractional-order time derivatives. The resulting schemes have a second-order accuracy in space and a $(2-\alpha )$-order accuracy in time, and the unconditional stability and convergence have been proved theoretically. In [19], Shi and Xie derived and analyzed the fourth-order compact block-centered finite-difference schemes for one-dimensional and two-dimensional variable-coefficient elliptic and parabolic problems. They demonstrated the stability of the solution and flux and performed optimal fourth-order error estimation.

The fourth-order parabolic problem has important practical significance in science and engineering. It can be used to describe bistable phenomena encountered in various fields [20], such as the competition and spatial sorting of biological populations, migration of riverbeds, charge-density distribution of quantum semiconductors, etc. [21,22]. Since the exact solution of the fourth-order equation is difficult to obtain, the numerical method of the fourth-order parabolic equation has attracted extensive attention from researchers in recent years. In [22], Jüngel studied the positivity-preserving numerical scheme for a class of fourth-order nonlinear parabolic systems in quantum semiconductor modeling and performed transient calculations using a macroscopic quantum model for the first time. The time-fractional derivative is especially good at describing dynamic processes with history dependence; therefore, the time-fractional differential equation can be used to depict physical problems with time variables with great accuracy such as in [23,24]. Fractional Caputo derivatives can be used to study the dynamics of plankton–fish models in the presence of toxic compounds produced by harmful algal blooms [25]. The time-fractional fourth-order parabolic equation can better describe the propagation of waves in intense laser beams and the charge-density distribution of quantum semiconductors. Currently, many researchers are dedicated to the study of fractional-order differential equations. Aziz [26] studied two inverse source problems of fourth-order parabolic equations with fractional time derivatives. Li and Liao [27] used a class of $L1$-Galerkin finite-element methods to study the numerical solution of time-fractional nonlinear parabolic problems. They provided the optimal error estimates of several fully discrete linearized Galerkin finite-element methods for solving nonlinear problems. The authors of [28] established a fully discrete weak Galerkin finite-element method for the initial boundary value problems of two-dimensional sub-diffusion equations with Caputo fractional time derivatives. In [29], Liu and Du proposed and discussed the finite-difference/finite-element method for solving nonlinear time-fractional fourth-order reaction and diffusion problems. A new implicit compact difference scheme for fourth-order fractional diffuse wave systems was constructed in [30]. In addition, Ji and Sun [31] studied the compact algorithm for a class of fourth-order fractional diffusion equations with first-order Dirichlet boundary conditions.

So far, no block-centered finite-difference methods for fourth-order parabolic equations with fractional-order time derivatives have been published in the literature. For the Neumann boundary conditions, which provide the boundary charge density, the time fractional fourth-order parabolic equation is more suitable to be solved using the block-centered finite-difference method, without separately considering the numerical solution of the nodes near the boundary. Therefore, it is of great theoretical and practical significance to propose and develop a block-centered finite-difference method for time-fractional fourth-order parabolic equations. This paper discusses the block-centered finite-difference method [30] for variable-coefficient fourth-order parabolic equations of fractional-order time derivatives with Neumann boundary conditions. In this method, the mixed finite-element method is used for theoretical analysis, which gives the error analysis a certain regularity. The fractional-order time derivatives are approximated by $L1$ interpolation. The block-centered finite-difference scheme is established, and the error estimations of the discrete $L2$ norm of the approximate solution and its derivatives are provided. Numerical examples are presented to verify the effectiveness of the block-centered finite-difference scheme.

This paper is organized as follows. Section 2 introduces the notations used in this paper. Section 3 presents the block-centered finite-difference scheme and error estimation for fourth-order ordinary differential equations. Section 4 establishes the block-centered finite-difference schemes for the fractional-order time derivatives and proves the stability and convergence of the schemes. In Section 5, numerical examples are provided to verify the convergence of the proposed schemes.

2. Notations

We first introduce some notations and definitions used in this paper, which will help with the following analysis. We use notations similar to those in [2]. Define the partition ${\mathcal{T}}_x$ of $\mathsf{\Omega }=(x_L,x_R)$ as

$${\mathcal{T}}_x:x_L=x_{1/2}<x_{3/2}<\cdots <x_{N-1/2}<x_{N+1/2}=x_R.$$

For each $i=1$ to N, define

$$\begin{array}{c}\begin{array}{cc}\hfill & x_i=\frac{1}{2}(x_{i+1/2}+x_{i-1/2}),\hfill \\ \hfill & h_i=x_{i+1/2}-x_{i-1/2},\hfill \\ \hfill & h=\underset{i}{max}{h}_i,\hfill \\ \hfill & h_{i+1/2}=x_{i+1}-x_i=\frac{1}{2}(h_i+h_{i+1}),\hfill \\ \hfill & {\mathsf{\Omega }}_i=(x_{i-1/2},x_{i+1/2}).\hfill \end{array}\hfill \end{array}$$

The block-centered dual partition grids are defined as ${\mathcal{T}}_x^{*}=\left\{x_i\right\}$.

Take a positive integer J, and let $\tau =\widehat{T}/J$, $t_n=n\tau$$(0\le n\le J)$.

For any function $g\left(x\right)$, let $g_i$, $g_{i+1/2}$, $g_i^n$ donate $g\left(x_i\right)$, $g\left(x_{i+1/2}\right)$, $g(x_i,t_n)$.

Define the following notations

$${\left[d_{x}g\right]}_{i+1/2}=\frac{g_{i+1}-g_i}{h_{i+1/2}},$$

$${\left[D_{x}g\right]}_i=\frac{g_{i+1/2}-g_{i-1/2}}{h_i}.$$

For functions F and G, define the midpoint quadrature formula and trapezoidal quadrature formula on ${\mathsf{\Omega }}_i$ as

$$\begin{array}{c}\hfill {(F,G)}_{M,{\mathsf{\Omega }}_i}=h_{i}F\left(x_i\right)G\left(x_i\right),\end{array}$$

and

$$\begin{array}{c}\hfill {(F,G)}_{T,{\mathsf{\Omega }}_i}=\frac{h_i}{2}[F\left(x_{i-1/2}\right)G\left(x_{i-1/2}\right)+F\left(x_{i+1/2}\right)G\left(x_{i+1/2}\right)].\end{array}$$

Given functions $f\left(x\right)$ and $g\left(x\right)$, define the $L2$ inner product and norm

$$\begin{array}{cc}\hfill & (f,g)={\int }_{\mathsf{\Omega }}f\left(x\right)g\left(x\right)dx,\hfill \\ \hfill & {\parallel f\parallel }^2=(f,f),\hfill \end{array}$$

and the discrete $L2$ inner products and norms

$$\begin{array}{cc}\hfill & {(f,g)}_M=\sum _{i=1}^N{h}_i{f}_i{g}_i,\hfill \\ \hfill & {\parallel f\parallel }_M^2={(f,f)}_M,\hfill \\ \hfill & {(f,g)}_T=\sum _{i=1}^N{h}_{i+1/2}{f}_{i+1/2}{g}_{i+1/2},\hfill \\ \hfill & {\parallel f\parallel }_T^2={(f,f)}_T,\hfill \\ \hfill & {\parallel f\parallel }_{L^{\infty }\left(\mathsf{\Omega }\right)}=\underset{1\le i\le N}{max}\left|f\left(x_i\right)\right|.\hfill \end{array}$$

Define $S_c^d\left({\mathcal{T}}_x\right)$ as the finite-dimensional space of one-dimensional functions that have c continuous derivatives on $(x_L,x_R)$ and are piecewise polynomials of degree d in each interval ${\mathsf{\Omega }}_i$. When $c=-1$, the functions themselves may be discontinuous.

The notation $\parallel g\parallel =O\left(h^k\right),k>0$ means that there exists a constant C such that $\parallel g\left(x\right)\parallel \le C{h}^k$ as h approaches zero.

3. Fourth-Order Ordinary Differential Equation

In order to discuss the block-centered difference method for the time-fractional fourth-order parabolic equation, we first consider the block-centered difference scheme for the fourth-order ordinary differential equation.

We consider fourth-order variable-coefficient ordinary differential equations with Neumann boundary conditions

$$\begin{array}{c}\left\{\begin{array}{c}{\left(a\left(x\right)w^{″}\right)}^{″}=f\left(x\right),x\in (x_L,x_R),\hfill \\ w^{\prime }{{|}_{x=x_L}=0,w^{\prime }|}_{x=x_R}=0,\hfill \\ {\left(a\left(x\right)w^{″}\right)}^{\prime }{{|}_{x=x_L}=0,{\left(a\left(x\right)w^{″}\right)}^{\prime }|}_{x=x_R}=0.\hfill \end{array}\right.\hfill \end{array}$$

(1)

where $f\left(x\right)$ is a known smooth function.

Let

$$v\left(x\right)=-w^{\prime }\left(x\right),p\left(x\right)=a\left(x\right)v^{\prime }\left(x\right),u\left(x\right)=-p^{\prime }\left(x\right),$$

then, we have

$$\begin{array}{cc}\hfill u^{\prime }\left(x\right)=f\left(x\right),& x\in (x_L,x_R),\hfill \\ \hfill u\left(x\right)=& -p^{\prime }\left(x\right),\hfill \\ \hfill a\left(x\right)v^{\prime }\left(x\right)=& p\left(x\right),\hfill \\ \hfill v\left(x\right)=& -w^{\prime }\left(x\right),\hfill \\ \hfill u\left(x\right)=0,v\left(x\right)=0,& atx=x_L,x=x_R.\hfill \end{array}$$

(2)

The block-centered finite-difference approximations $U_{i+1/2}$, $V_{i+1/2}$, $W_i$, and $P_i$ to $u\left(x_{i+1/2}\right)$, $v\left(x_{i+1/2}\right)$, $w\left(x_i\right)$, and $p\left(x_i\right)(i=1,2,\cdots ,N)$, respectively, satisfy the following

$$\begin{array}{cc}\hfill {\left[D_{x}U\right]}_i=f_i,& U_{1/2}=0,U_{N+1/2}=0,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill U_{i+1/2}=& {[-d_{x}P]}_{i+1/2},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\left[D_{x}V\right]}_i=\frac{P_i}{a_i},& V_{1/2}=0,V_{N+1/2}=0,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill V_{i+1/2}& ={[-d_{x}W]}_{i+1/2},\hfill \end{array}$$

(6)

which approximate the original Equation (2). The above block-centered finite-difference scheme can be written as a mixed finite-element scheme with approximate integration

$$\begin{array}{c}\hfill {({\left(U\right)}^{\prime },1)}_{{\mathsf{\Omega }}_i}={(f,1)}_{M,{\mathsf{\Omega }}_i},\end{array}$$

(7)

$$\begin{array}{c}\hfill {(U,\chi )}_T-(P,{\chi }^{\prime })=0,\chi \in \tilde{S},\end{array}$$

(8)

$$\begin{array}{c}\hfill {({\left(V\right)}^{\prime },1)}_{{\mathsf{\Omega }}_i}={(\frac{P}{a},1)}_{M,{\mathsf{\Omega }}_i},\end{array}$$

(9)

$$\begin{array}{c}\hfill {(V,\chi )}_T-(W,{\chi }^{\prime })=0,\chi \in \tilde{S}.\end{array}$$

(10)

where $\tilde{S}=S_0^1\left({\mathcal{T}}_x\right)\bigcap \{\chi :\chi \left(x_L\right)=\chi \left(x_R\right)=0\}$, U and V are in $\tilde{S}$, and P and W are in $S_{-1}^0\left({\mathcal{T}}_x\right)$.

 Lemma 1 

([2]). If $w^{\left(5\right)}\left(x\right)$ is continuous and $f^{‴}\left(x\right)$ is in $L^1\left({\mathsf{\Omega }}_i\right)$ for all i,

$$\begin{array}{c}\hfill {\parallel U-u\parallel }_T=O\left(h^2\right),{\parallel P-p\parallel }_M=O\left(h^2\right),\\ \hfill {\parallel V-v\parallel }_T=O\left(h^2\right),{\parallel W-w\parallel }_M=O\left(h^2\right).\end{array}$$

(11)

 Proof.

By using Equations (2)–(6) and the Taylor expansion, we can obtain

$$\begin{array}{cc}\hfill u_{i+1/2}-U_{i+1/2}& ={\int }_{x_{1/2}}^{x_{i+1/2}}(u^{\prime }-U^{\prime })dx={\int }_{x_{1/2}}^{x_{i+1/2}}(f\left(x\right)-f\left(x_i\right))dx=\sum _{j=1}^i{\int }_{{\mathsf{\Omega }}_j}(f\left(x\right)-f\left(x_j\right))dx\hfill \\ \hfill & =\sum _{j=1}^i{\int }_{{\mathsf{\Omega }}_j}((x-x_j)f_j^{\prime }+\frac{{(x-x_j)}^2}{2}{f}_j^{″}+{\int }_{x_j}^x\frac{{(x-s)}^2}{2}{f}^{‴}\left(s\right)ds)dx\hfill \\ \hfill & =\sum _{j=1}^i(\frac{1}{24}{h}_j^3{f}_j^{″}+{\int }_{{\mathsf{\Omega }}_j}{\int }_{x_j}^x\frac{{(x-s)}^2}{2}{f}^{‴}\left(s\right)dsdx)\hfill \\ \hfill & \le \frac{1}{24}\sum _{j=1}^i{h}_j^3{f}_j^{″}+|\sum _{j=1}^i{\int }_{{\mathsf{\Omega }}_j}{\int }_{x_j}^x\frac{{(x-s)}^2}{2}{f}^{‴}\left(s\right)dsdx|\hfill \\ \hfill & \le \frac{1}{24}\sum _{j=1}^i{h}_j^3{f}_j^{″}+\sum _{j=1}^i\frac{h_j^3}{8}{\int }_{{\mathsf{\Omega }}_j}\left|f^{‴}\right|dx.\hfill \end{array}$$

(12)

Then, we have

$$|u_{i+1/2}-U_{i+1/2}|\le \sum _{j=1}^i\frac{h_j^2}{2}{\parallel f^{″}\parallel }_{L^1}.$$

$$\begin{array}{cc}\hfill p_{i+1}-P_{i+1}& =-\sum _{j=1}^i{\int }_{x_j}^{x_{j+1}}(u-U_{j+1/2})dx-{\int }_0^{x_1}udx\hfill \\ \hfill & =-{\int }_0^{x_1}udx+\sum _{j=1}^i{\int }_{x_j}^{x_{j+1}}{(U-u)}_{j+1/2}dx+\sum _{j=1}^i{\int }_{x_j}^{x_{j+1}}(u_{j+1/2}-u)dx.\hfill \end{array}$$

According to Equation (12), the first two terms on the right side of the equation are $O\left(h^2\right)$. Now, we estimate the third term on the right side.

$$\begin{array}{cc}\hfill \sum _{j=1}^i{\int }_{x_j}^{x_{j+1}}(p_{j+1/2}^{\prime }-p^{\prime })dx& =-\sum _{j=1}^i{\int }_{x_j}^{x_{j+1}}((x-x_{j+1/2})p_{j+1/2}^{\prime \prime }+{\int }_{x_{j+1/2}}^x(x-s)p^{\prime \prime \prime }\left(s\right)ds)dx\hfill \\ \hfill & \le |\sum _{j=2}^i\frac{h_j^2}{8}|p_{j+1/2}^{\prime \prime }-p_{j-1/2}^{\prime \prime }|+\frac{h_1^2}{8}|p_{3/2}^{\prime \prime }|+\frac{h_{i+1}^2}{8}|p_{i+1/2}^{\prime \prime }|\hfill \\ \hfill & +\sum _{j=1}^i\frac{h_j^2}{2}{\int }_{x_j}^{x_{j+1}}\left|p^{\prime \prime \prime }\right|ds\hfill \\ \hfill & =O\left(h^2\right).\hfill \end{array}$$

We can obtain

$$|p_{i+1}-P_{i+1}|=O\left(h^2\right).$$

Similarly, we can reach the same conclusion for V and W by employing the Taylor expansion. □

The second-order error estimation for a block-centered difference scheme applied to fourth-order ordinary differential equations has been derived.

4. Time-Fractional Fourth-Order Parabolic Equation

In this section, we consider the block-centered finite-difference method for a time-fractional fourth-order parabolic equation when $0<\alpha <1$.

We consider the following variable-coefficient fractional fourth-order parabolic problem with initial and boundary value conditions

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{}_0^c{D}_t^{\alpha }w(x,t)+\mu \frac{{\partial }^2}{\partial x^2}(a(x)\frac{{\partial }^{2}w(x,t)}{\partial x^2})=f(x,t),x\in (x_L,x_R),t\in (0,\widehat{T}],\hfill \\ w(x,0)=\phi \left(x\right),x\in (x_L,x_R),t\in (0,\widehat{T}],\hfill \\ \frac{\partial w(x,t)}{\partial x}{{|}_{x=x_L}=0,\frac{\partial w(x,t)}{\partial x}|}_{x=x_R}=0,t\in (0,\widehat{T}],\hfill \\ \frac{\partial }{\partial x}(a(x)\frac{{\partial }^{2}w(x,t)}{\partial x^2}){{|}_{x=x_L}=0,\frac{\partial }{\partial x}\left(a\left(x\right)\frac{{\partial }^{2}w(x,t)}{\partial x^2}\right)|}_{x=x_R}=0,t\in (0,\widehat{T}].\hfill \end{array}\right.\end{array}$$

(13)

where $\mu$ is a constant; $a\left(x\right)$, $f(x,t)$, and $\phi \left(x\right)$ are known smooth functions; and it is assumed that $0<a_0\le a\left(x\right)\le a_1$.

We consider the case $0<\alpha <1$. ${}_0^c{D}_t^{\alpha }w$ in (13) is defined as the Caputo fractional derivative of $\alpha$, which is given by

$$\begin{array}{c}\hfill {}_0^c{D}_t^{\alpha }w(x,t)=\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^t{(t-\tau )}^{-\alpha }\frac{\partial w(x,\tau )}{\partial \tau }d\tau .\end{array}$$

4.1. Block-Centered Finite-Difference Scheme

In this subsection, we provide the block-centered difference scheme for a time-fractional fourth-order parabolic equation.

 Lemma 2 

([32]). Suppose $f\left(t\right)\in C^2[0,\widehat{T}]$, $0<\alpha <1$,

$$\begin{array}{c}\hfill {}_0^c{D}_t^{\alpha }f\left(t_n\right)=D_{\tau }^{\alpha }f\left(t_n\right)+R\left(f\left(t_n\right)\right)=\frac{{\tau }^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}[a_0^{\left(\alpha \right)}f\left(t_n\right)-\sum _{k=1}^{n-1}(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)})f\left(t_k\right)-a_{n-1}^{\left(\alpha \right)}f\left(t_0\right)]\\ \hfill +R\left(f\left(t_n\right)\right),\end{array}$$

where

$$\begin{array}{cc}\hfill & a_k^{\left(\alpha \right)}={(k+1)}^{1-\alpha }-k^{1-\alpha },\hfill \\ \hfill & |R\left(f\left(t_n\right)\right)|\le \frac{1}{2\mathsf{\Gamma }(1-\alpha )}[\frac{1}{4}+\frac{\alpha }{(1-\alpha )(2-\alpha )}]\underset{t_0\le t\le t_n}{max}\left|f^{\prime \prime }\left(t\right)\right|{\tau }^{2-\alpha }.\hfill \end{array}$$

 Lemma 3 

([32]). Given that $0<\alpha <1$, we have $\underset{l\to +\infty }{lim}{a}_l=0$,

$$\begin{array}{c}\hfill 1=a_0^{\left(\alpha \right)}>a_1^{\left(\alpha \right)}>a_2^{\left(\alpha \right)}>\cdots >a_l^{\left(\alpha \right)}>0.\end{array}$$

The block-centered finite-difference method (13) defines ${\left\{U_{i+1/2}^n\right\}}_{n=1}^M$, ${\left\{V_{i+1/2}^n\right\}}_{n=1}^M$, ${\left\{W_i^n\right\}}_{n=1}^M$ and ${\left\{P_i^n\right\}}_{n=1}^M(i=1,2,\cdots ,N)$, satisfying

$$\begin{array}{cc}\hfill {\left[D_{\tau }^{\alpha }W\right]}_i^n+\mu {\left[D_{x}U\right]}_i^n=f_i^n,& \alpha \in (0,1),U_{1/2}^n=0,U_{N+1/2}^n=0,\hfill \end{array}$$

(14)

$$\begin{array}{c}\hfill U_{i+1/2}^n={[-d_{x}P]}_{i+1/2}^n,\end{array}$$

(15)

$$\begin{array}{cc}\hfill {\left[D_{x}V\right]}_i^n=\frac{P_i^n}{a_i},& V_{1/2}^n=0,V_{N+1/2}^n=0,\hfill \end{array}$$

(16)

$$\begin{array}{c}\hfill V_{i+1/2}^n={[-d_{x}W]}_{i+1/2}^n,\end{array}$$

(17)

where $U^0=E_h{u}^0$, $P^0=E_h{p}^0$, $V^0=E_h{v}^0$, and $W^0=E_h{w}^0$. Here, $U_{i+1/2}^n$, $P_i^n$, $V_{i+1/2}^n$, and $W_i^n$ are approximations to $u_{i+1/2}^n$, $p_i^n$, $v_{i+1/2}^n$, and $w_i^n$, respectively, and $E_h{u}^n$, $E_h{p}^n$, $E_h{v}^n$, and $E_h{w}^n$ are their corresponding elliptic projections.

The above block-centered finite-difference scheme can be written as a mixed finite-element scheme with approximate integration

$$\begin{array}{cc}\hfill {(D_{\tau }^{\alpha }{W}^n,1)}_{{\mathsf{\Omega }}_i}+\mu {({\left(U^n\right)}^{\prime },1)}_{{\mathsf{\Omega }}_i}={(f^n,1)}_{M,{\mathsf{\Omega }}_i},& \alpha \in (0,1),\hfill \end{array}$$

(18)

$$\begin{array}{c}\hfill {(U^n,\chi )}_T-(P^n,{\chi }^{\prime })=0,\chi \in \tilde{S},\end{array}$$

(19)

$$\begin{array}{c}\hfill {({\left(V^n\right)}^{\prime },1)}_{{\mathsf{\Omega }}_i}={(\frac{P^n}{a},1)}_{M,{\mathsf{\Omega }}_i},\end{array}$$

(20)

$$\begin{array}{c}\hfill {(V^n,\chi )}_T-(W^n,{\chi }^{\prime })=0,\chi \in \tilde{S}.\end{array}$$

(21)

4.2. Stability Analysis

In this subsection, we prove the stability of the scheme (14)–(17) when $0<\alpha <1$.

 Theorem 1.

For the block-centered difference scheme, the following stable inequality holds unconditionally for sufficiently small τ

$$\parallel W^n{\parallel }_M\le C_1\parallel W^0{\parallel }_M+C_2\underset{1\le k\le J}{max}{\parallel f^k\parallel }_M.$$

 Proof. 

For $n=1$, by multiplying (14) by $h_i{W}_i^1$ and summing on i from 1 to N, we obtain ${(D_{\tau }^{\alpha }{W}^1,W^1)}_M+\mu {(W_{xx}^1,W_{xx}^1)}_M={(f^1,W^1)}_M$. So,

$${(a_0^{\left(\alpha \right)}{W}^1,W^1)}_M\le {(a_0^{\left(\alpha \right)}{W}^0,W^1)}_M+\rho {(f^1,W^1)}_M,$$

where $\rho ={\tau }^{\alpha }\mathsf{\Gamma }(2-\alpha )$. Using the Cauchy–Schwarz inequality and Young inequality, we can obtain

$$(1-\frac{\rho }{2})\parallel W^1{\parallel }_M^2\le \parallel W^0{\parallel }_M^2+2\rho {\parallel f^1\parallel }_M^2.$$

So, we have

$$\parallel W^1{\parallel }_M^2\le \frac{2}{2-\rho }\parallel W^0{\parallel }_M^2+\frac{4\rho }{2-\rho }{\parallel f^1\parallel }_M^2.$$

Thus, we can obtain

$$\parallel W^1{\parallel }_M\le C_1\parallel W^0{\parallel }_M+C_2{\parallel f^1\parallel }_M.$$

For $n\ge 2$, we suppose that the stability conclusion of the difference scheme is valid when $k\le n-1$.

Then, by multiplying (14) by $h_i{W}_i^n$ and summing on i from 1 to N, we can obtain ${(D_{\tau }^{\alpha }{W}^n,W^n)}_M$$\le {(f^n,W^n)}_M$. So,

$${(a_0^{\left(\alpha \right)}{W}^n,W^n)}_M\le {(\sum _{k=1}^{n-1}(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)})W^k,W^n)}_M+{(a_{n-1}^{\left(\alpha \right)}{W}^0,W^n)}_M+\rho {(f^n,W^n)}_M.$$

Using the Cauchy–Schwarz inequality and Young inequality, we obtain

$$\parallel W^n{\parallel }_M^2\le \frac{2}{2-\rho }\sum _{k=1}^{n-1}(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)})\parallel W^k{\parallel }_M^2+\frac{2}{2-\rho }{a}_{n-1}^{\left(\alpha \right)}\parallel W^0{\parallel }_M^2+\frac{4\rho }{2-\rho }{\parallel f^n\parallel }_M^2.$$

Through mathematical induction and using the relation of coefficient $a_k$, we have

$$\begin{array}{cc}\hfill \parallel W^n{\parallel }_M^2& \le \frac{2}{2-\rho }\sum _{k=1}^{n-1}(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)})(C_1^{\prime }\parallel W^0{\parallel }_M^2+C_2^{\prime }\underset{1\le k\le J}{max}\parallel f^k{\parallel }_M^2)+\frac{2}{2-\rho }{a}_{n-1}^{\left(\alpha \right)}{\parallel W^0\parallel }_M^2\hfill \\ & +\frac{4\rho }{2-\rho }{\parallel f^n\parallel }_M^2\hfill \\ \hfill & \le \frac{2}{2-\rho }[(C_1^{\prime }(a_0^{\left(\alpha \right)}-a_{n-1}^{\left(\alpha \right)})+a_{n-1}^{\left(\alpha \right)})\parallel W^0{\parallel }_M^2+(C_2^{\prime }(a_0^{\left(\alpha \right)}-a_{n-1}^{\left(\alpha \right)})+2\rho )\underset{1\le k\le J}{max}\parallel f^k{\parallel }_M^2].\hfill \end{array}$$

There are constants $C_3$ and $C_4$ that make

$$\parallel W^n{\parallel }_M^2\le C_3\parallel W^0{\parallel }_M^2+C_4\underset{1\le k\le J}{max}{\parallel f^k\parallel }_M^2.$$

So, we can obtain

$$\parallel W^n{\parallel }_M\le C_1\parallel W^0{\parallel }_M+C_2\underset{1\le k\le J}{max}{\parallel f^k\parallel }_M.$$

We complete the proof. □

4.3. Error Analysis

The error analysis of the block-centered difference scheme for the time-fractional parabolic equation is performed.

Error estimates for the finite-difference scheme of (14)–(17) are derived using a technique of mixed finite-element methods for parabolic partial differential equations.

$$\begin{array}{cc}\hfill \mu \frac{{\partial }^2}{\partial x^2}(a\left(x\right)\frac{{\partial }^2{w}^n}{\partial x^2})=\varphi =f^n{-}_0^c{D}_t^{\alpha }{w}^n,& \alpha \in (0,1),x\in (x_L,x_R),\hfill \end{array}$$

(22)

$$\begin{array}{c}\hfill \frac{\partial w^n}{\partial x}{{|}_{x=x_L}=0,\frac{\partial w^n}{\partial x}|}_{x=x_R}=0,\end{array}$$

(23)

$$\begin{array}{c}\hfill \frac{\partial }{\partial x}(a\left(x\right)\frac{{\partial }^2{w}^n}{\partial x^2}){{|}_{x=x_L}=0,\frac{\partial }{\partial x}(a\left(x\right)\frac{{\partial }^2{w}^n}{\partial x^2})|}_{x=x_R}=0.\end{array}$$

(24)

For fixed n, let $E_h{u}_{i+1/2}^n,E_h{v}_{i+1/2}^n,E_h{w}_i^n$, and $E_h{p}_i^n$ be defined by

$$\begin{array}{cc}\hfill \mu {\left[D_x{E}_{h}u\right]}_i^n=f_i^n-{{(}_0^c{D}_t^{\alpha }w)}_i^n,& E_h{u}_{1/2}^n=0,E_h{u}_{N+1/2}^n=0,\hfill \end{array}$$

(25)

$$\begin{array}{c}\hfill E_h{u}_{i+1/2}^n={[-d_x{E}_{h}p]}_{i+1/2}^n,\end{array}$$

(26)

$$\begin{array}{cc}\hfill {\left[D_x{E}_{h}v\right]}_i^n=\frac{E_h{p}_i^n}{a_i},& E_h{v}_{1/2}^n=0,E_h{v}_{N+1/2}^n=0,\hfill \end{array}$$

(27)

$$\begin{array}{c}\hfill E_h{v}_{i+1/2}^n={[-d_x{E}_{h}w]}_{i+1/2}^n,\end{array}$$

(28)

where $E_h{u}^0={\left(a\left(x\right){\phi }^{\prime \prime }\left(x\right)\right)}^{\prime }$, $E_h{p}^0=-a\left(x\right){\phi }^{\prime \prime }\left(x\right)$, $E_h{v}^0=-{\phi }^{\prime }\left(x\right)$, and $E_h{w}^0=\phi \left(x\right)$.

Equations (25)–(28) can be written as a mixed finite-element method with approximate integration

$$\begin{array}{c}\hfill \mu {({\left(E_h{u}^n\right)}^{\prime },1)}_{{\mathsf{\Omega }}_i}={(f^n,1)}_{M,{\mathsf{\Omega }}_i}-{{(}_0^c{D}_t^{\alpha }{w}^n,1)}_{M,{\mathsf{\Omega }}_i},\end{array}$$

(29)

$$\begin{array}{c}\hfill {(E_h{u}^n,\chi )}_T-(E_h{p}^n,{\chi }^{\prime })=0,\chi \in \tilde{S},\end{array}$$

(30)

$$\begin{array}{c}\hfill {({\left(E_h{v}^n\right)}^{\prime },1)}_{{\mathsf{\Omega }}_i}={(\frac{E_h{p}^n}{a},1)}_{M,{\mathsf{\Omega }}_i},\end{array}$$

(31)

$$\begin{array}{c}\hfill {(E_h{v}^n,\chi )}_T-(E_h{w}^n,{\chi }^{\prime })=0,\chi \in \tilde{S}.\end{array}$$

(32)

By the error of the ellipse projection, we have

$$\begin{array}{c}\hfill \parallel E_h{u}^n-u^n{\parallel }_T=O\left(h^2\right),{\parallel E_h{p}^n-p^n\parallel }_M=O\left(h^2\right),\\ \hfill \parallel E_h{v}^n-v^n{\parallel }_T=O\left(h^2\right),{\parallel E_h{w}^n-w^n\parallel }_M=O\left(h^2\right),\end{array}$$

(33)

which hold for sufficiently smooth w.

By differentiating t in Equations (25)–(28), we can obtain the following estimation

$$\begin{array}{c}\hfill \parallel E_h{w}_t^n-w_t^n{\parallel }_M=O\left(h^2\right).\end{array}$$

(34)

Set ${\xi }_{i+1/2}^n=U_{i+1/2}^n-E_h{u}_{i+1/2}^n$, ${\eta }_i^n=P_i^n-E_h{p}_i^n$, ${\theta }_{i+1/2}^n=V_{i+1/2}^n-E_h{v}_{i+1/2}^n$, ${\zeta }_i^n=W_i^n-E_h{w}_i^n$, and ${\sigma }_i^n=E_h{w}_i^n-w_i^n.$

By subtracting (25) from (14), we obtain

$$\begin{array}{cc}\hfill D_{\tau }^{\alpha }{\zeta }_i^n+\mu {\left[D_x\xi \right]}_i^n=R\left(w_i^n\right)-D_{\tau }^{\alpha }{\sigma }_i^n,& \alpha \in (0,1),{\xi }_{1/2}^n=0,{\xi }_{N+1/2}^n=0.\hfill \end{array}$$

(35)

By subtracting (26), (27) and (28) from (15), (16) and (17), respectively, we obtain

$$\begin{array}{c}\hfill {\xi }_{i+1/2}^n={[-d_x\eta ]}_{i+1/2}^n,\end{array}$$

(36)

$$\begin{array}{cc}\hfill {\eta }_i^n=a_i{\left[D_x\theta \right]}_i^n,& {\theta }_{1/2}^n=0,{\theta }_{N+1/2}^n=0,\hfill \end{array}$$

(37)

$$\begin{array}{c}\hfill {\theta }_{i+1/2}^n={[-d_x\zeta ]}_{i+1/2}^n.\end{array}$$

(38)

By multiplying (35) by $h_i{\zeta }_i^n$ and summing on i from 1 to N, we deduce that

$$\begin{array}{c}\hfill {(D_{\tau }^{\alpha }{\zeta }^n,{\zeta }^n)}_M+\mu {(D_x{\xi }^n,{\zeta }^n)}_M={(R\left(w^n\right),{\zeta }^n)}_M-{(D_{\tau }^{\alpha }{\sigma }^n,{\zeta }^n)}_M.\end{array}$$

(39)

By (36)–(38), and Lemma 2, we have

$$\begin{array}{cc}\hfill {(D_x{\xi }^n,{\zeta }^n)}_M& ={(D_x(-d_x{\eta }^n),{\zeta }^n)}_M\hfill \\ \hfill & ={(d_x{\eta }^n,d_x{\zeta }^n)}_T\hfill \\ \hfill & ={({\eta }^n,D_x{\theta }^n)}_M\hfill \\ \hfill & =\frac{1}{a}{({\eta }^n,{\eta }^n)}_M.\hfill \end{array}$$

(40)

Now,

$$\begin{array}{c}\hfill {(D_{\tau }^{\alpha }{\zeta }^n,{\zeta }^n)}_M=\frac{{\tau }^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}[a_0^{\left(\alpha \right)}{({\zeta }^n,{\zeta }^n)}_M-\sum _{k=1}^{n-1}(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}){({\zeta }^k,{\zeta }^n)}_M-a_{n-1}^{\left(\alpha \right)}{({\zeta }^0,{\zeta }^n)}_M].\end{array}$$

(41)

Let $\rho ={\tau }^{\alpha }\mathsf{\Gamma }(2-\alpha )$, and (39) can be written as

$$\begin{array}{c}\hfill \parallel {\zeta }^n{\parallel }_M^2+\frac{\mu \rho }{a_1}{\parallel {\eta }^n\parallel }_M^2\le \sum _{k=1}^{n-1}(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}){({\zeta }^k,{\zeta }^n)}_M+a_{n-1}^{\left(\alpha \right)}{({\zeta }^0,{\zeta }^n)}_M\\ \hfill +\rho |{(R\left(w^n\right),{\zeta }^n)}_M|+\rho |{(D_{\tau }^{\alpha }{\sigma }^n,{\zeta }^n)}_M|.\end{array}$$

(42)

By the Cauchy–Schwarz inequality and Young inequality, we have

$$\begin{array}{c}\hfill (1-\frac{\rho }{2})\parallel {\zeta }^n{\parallel }_M^2+\frac{2}{a_1}\mu \rho \parallel {\eta }^n{\parallel }_M^2\le \sum _{k=1}^{n-1}(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)})\parallel {\zeta }^k{\parallel }_M^2+a_{n-1}^{\left(\alpha \right)}{\parallel {\zeta }^0\parallel }_M^2\\ \hfill +4\rho \parallel D_{\tau }^{\alpha }{\sigma }^n{\parallel }_M^2+4\rho {\parallel R\left(w^n\right)\parallel }_M^2.\end{array}$$

(43)

According to the definition of the fractional derivative and Equation (34),

$$\begin{array}{cc}\hfill \parallel D_{\tau }^{\alpha }{\sigma }^n{\parallel }_M& =\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\parallel \sum _{k=1}^n\frac{{\sigma }^k-{\sigma }^{k-1}}{\tau }{\int }_{t_{k-1}}^{t_k}{(t_n-t)}^{-\alpha }dt\parallel }_M\hfill \\ \hfill & =\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\parallel \sum _{k=1}^n\frac{1}{\tau }{\int }_{t_{k-1}}^{t_k}{\sigma }_{t}dt{\int }_{t_{k-1}}^{t_k}{(t_n-t)}^{-\alpha }dt\parallel }_M\hfill \\ \hfill & \le \frac{1}{\mathsf{\Gamma }(1-\alpha )}\underset{0\le t\le t_n}{max}{\parallel {\sigma }_t\parallel }_M\sum _{k=1}^n{\int }_{t_{k-1}}^{t_k}{(t_n-t)}^{-\alpha }dt\hfill \\ \hfill & \le \frac{T^{1-\alpha }}{\mathsf{\Gamma }(2-\alpha )}\underset{0\le t\le t_n}{max}{\parallel {\sigma }_t\parallel }_M\hfill \\ \hfill & \le C{h}^2.\hfill \end{array}$$

(44)

Notice that ${\zeta }^0=0$. Using Theorem 1 and the inductive hypothesis, we deduce that

$$\begin{array}{c}\hfill (1-\frac{\rho }{2})\parallel {\zeta }^n{\parallel }_M^2+\frac{2}{a_1}\mu \rho {\parallel {\eta }^n\parallel }_M^2\le C({\tau }^{4-2\alpha }+h^4).\end{array}$$

(45)

Thus, we have

$$\begin{array}{c}\hfill \parallel {\zeta }^n{\parallel }_M^2\le C({\tau }^{4-2\alpha }+h^4),\end{array}$$

(46)

$$\begin{array}{c}\hfill \parallel {\eta }^n{\parallel }_M^2\le C({\tau }^{4-2\alpha }+h^4).\end{array}$$

(47)

By multiplying (35) by $h_i{\eta }_i^n$ and summing on i from 1 to N, we can obtain

$$\begin{array}{c}\hfill {(D_{\tau }^{\alpha }{\zeta }^n,{\eta }^n)}_M+\mu {(D_x{\xi }^n,{\eta }^n)}_M={(R\left(w^n\right),{\eta }^n)}_M-{(D_{\tau }^{\alpha }{\sigma }^n,{\eta }^n)}_M.\end{array}$$

(48)

Now,

$$\begin{array}{c}\hfill {(D_{\tau }^{\alpha }{\zeta }^n,{\eta }^n)}_M=\frac{{\tau }^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}[a_0^{\left(\alpha \right)}{({\zeta }^n,{\eta }^n)}_M-\sum _{k=1}^{n-1}(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}){({\zeta }^k,{\eta }^n)}_M-a_{n-1}^{\left(\alpha \right)}{({\zeta }^0,{\eta }^n)}_M].\end{array}$$

(49)

By using (36)–(38), and Lemma 2, we derive

$$\begin{array}{c}\hfill {({\zeta }^n,{\eta }^n)}_M=a{({\theta }^n,{\theta }^n)}_T,\end{array}$$

(50)

$$\begin{array}{c}\hfill {(D_x{\xi }^n,{\eta }^n)}_M={({\xi }^n,{\xi }^n)}_T.\end{array}$$

(51)

By substituting (49)–(51) into (48), (48) can be deformable to

$$\begin{array}{c}\hfill a_0\parallel {\theta }^n{\parallel }_T^2+\mu \rho {\parallel {\xi }^n\parallel }_T^2\le a_1\sum _{k=1}^{n-1}(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)}){({\theta }^k,{\theta }^n)}_T+a_1{a}_{n-1}^{\left(\alpha \right)}{({\theta }^0,{\theta }^n)}_T\\ \hfill +\rho |{(R\left(w^n\right),{\eta }^n)}_M|+\rho |{(D_{\tau }^{\alpha }{\sigma }^n,{\eta }^n)}_M|.\end{array}$$

(52)

By the Cauchy–Schwarz inequality, we derive

$$\begin{array}{c}\hfill (2a_0-a_1)\parallel {\theta }^n{\parallel }_T^2+2\mu \rho \parallel {\xi }^n{\parallel }_T^2\le a_1\sum _{k=1}^{n-1}(a_{n-k-1}^{\left(\alpha \right)}-a_{n-k}^{\left(\alpha \right)})\parallel {\theta }^k{\parallel }_T^2+a_1{a}_{n-1}^{\left(\alpha \right)}\parallel {\theta }^0{\parallel }_T^2+2\rho {\parallel {\eta }^n\parallel }_M^2+\\ \hfill \rho \parallel D_{\tau }^{\alpha }{\sigma }^n{\parallel }_M^2+\rho {\parallel R\left(w^n\right)\parallel }_M^2.\end{array}$$

(53)

Similarly, when $2a_0>a_1$, using the above mathematical induction and substituting (44) and (47) into (53), we can obtain

$$\begin{array}{c}\hfill (2a_0-a_1)\parallel {\theta }^n{\parallel }_T^2+2\mu \rho {\parallel {\xi }^n\parallel }_T^2\le C({\tau }^{4-2\alpha }+h^4).\end{array}$$

(54)

So, we can obtain

$$\parallel {\theta }^n{\parallel }_T^2\le C({\tau }^{4-2\alpha }+h^4),$$

$$\parallel {\xi }^n{\parallel }_T^2\le C({\tau }^{4-2\alpha }+h^4).$$

By (33) and the triangle inequality,

$$\parallel U^n-u^n{\parallel }_T^2\le C({\tau }^{4-2\alpha }+h^4),{\parallel P^n-p^n\parallel }_M^2\le C({\tau }^{4-2\alpha }+h^4),$$

$$\parallel V^n-v^n{\parallel }_T^2\le C({\tau }^{4-2\alpha }+h^4),{\parallel W^n-w^n\parallel }_M^2\le C({\tau }^{4-2\alpha }+h^4),$$

which hold for sufficiently smooth w.

We can draw the following conclusion.

 Theorem 2.

Let w be sufficiently smooth and satisfy (13). If U, P, V, and W satisfy (14)–(17), for all $n,n=1,2,\cdots ,J$,

$$\parallel U^n-u^n{\parallel }_T=O({\tau }^{2-\alpha }+h^2),{\parallel P^n-p^n\parallel }_M=O({\tau }^{2-\alpha }+h^2),$$

$$\parallel V^n-v^n{\parallel }_T=O({\tau }^{2-\alpha }+h^2),{\parallel W^n-w^n\parallel }_M=O({\tau }^{2-\alpha }+h^2).$$

5. Numerical Experiments

In this section, we present three numerical examples to verify the effectiveness and convergence of the block-centered finite-difference method.

 Example 1.

Consider the following fourth-order ordinary differential equation

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{w}^{″″}=f\left(x\right),x\in (0,1),\hfill \\ w^{\prime }{{|}_{x=0}=0,w^{\prime }|}_{x=1}=0,x\in (0,1),\hfill \\ w^{\prime \prime \prime }{{|}_{x=0}=0,w^{\prime \prime \prime }|}_{x=1}=0,x\in (0,1).\hfill \end{array}\right.\end{array}$$

(55)

We provide the exact solution to the problem $w\left(x\right)=cos\left(\pi x\right)$ and a source term $f\left(x\right)={\pi }^{4}cos\left(\pi x\right).$

The original equation is discretized by the block-centered difference scheme, resulting in Equations (3)–(6). Equations (4)–(6) are substituted into Equation (3) to form a system of linear equations about the unknown function W. The approximate solution of W is obtained using MATLAB software, whereas the approximate solutions of V, P, and U are obtained using Equations (4)–(6).

Take the spatial step $h=1/N.$ The space errors and convergence orders of W, V, P, and U are shown in

Table 1. The computational errors and convergence orders of W in space.

h	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
1/8	$2.5567\times {10}^{-2}$	-	$1.8434\times {10}^{-2}$	-
1/16	$6.4172\times {10}^{-3}$	1.9943	$4.5596\times {10}^{-3}$	2.0154
1/32	$1.6059\times {10}^{-3}$	1.9986	$1.1369\times {10}^{-3}$	2.0038
1/64	$4.0156\times {10}^{-4}$	1.9997	$2.8403\times {10}^{-4}$	2.0001

,

Table 2. The computational errors and convergence orders of V in space.

h	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
1/8	$6.1226\times {10}^{-2}$	-	$4.3293\times {10}^{-2}$	-
1/16	$1.5181\times {10}^{-2}$	2.0119	$1.0735\times {10}^{-2}$	2.0019
1/32	$3.7875\times {10}^{-3}$	2.0030	$2.6782\times {10}^{-3}$	2.0030
1/64	$9.4640\times {10}^{-4}$	2.0007	$6.6920\times {10}^{-4}$	2.0007

,

Table 3. The computational errors and convergence orders of P in space.

h	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
1/8	$1.2536\times {10}^{-1}$	-	$9.0382\times {10}^{-2}$	-
1/16	$3.1617\times {10}^{-2}$	1.9873	$2.2465\times {10}^{-2}$	2.0084
1/32	$7.9214\times {10}^{-3}$	1.9969	$5.6081\times {10}^{-3}$	2.0021
1/64	$1.9814\times {10}^{-3}$	1.9992	$1.4015\times {10}^{-3}$	2.0005

and

Table 4. The computational errors and convergence orders of U in space.

h	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
1/8	$2.0013\times {10}^{-1}$	-	$1.4151\times {10}^{-1}$	-
1/16	$4.9864\times {10}^{-2}$	2.0049	$3.5260\times {10}^{-2}$	2.0049
1/32	$1.2455\times {10}^{-2}$	2.0012	$8.8074\times {10}^{-3}$	2.0012
1/64	$3.1132\times {10}^{-3}$	2.0003	$2.2014\times {10}^{-3}$	2.0003

, respectively. The second and fourth columns in Table 1, Table 2, Table 3 and Table 4 show the maximum norm error and discrete L2 norm error, respectively, whereas the third and fifth columns show their corresponding spatial convergence orders, respectively. It can be seen from the table that the order of spatial convergence reaches the second order for both the maximum norm and discrete L2 error, which verifies the convergence of the block-centered finite-difference scheme.

For the functions $S\left(x\right)$ and $s\left(x\right)$, $L^{\infty }-Error$ and $L^2-Error$, as defined below, are used

$$\begin{array}{cc}\hfill & {\parallel S-s\parallel }_{L^{\infty }\left(\mathsf{\Omega }\right)}=\underset{1\le i\le N}{max}|S\left(x_i\right)-s\left(x_i\right)|,\hfill \\ \hfill & {\parallel S-s\parallel }_{L^2\left(\mathsf{\Omega }\right)}=(\sum _{i=1}^N|S\left(x_i\right)-s\left(x_i\right){{|}^{2}h)}^{1/2}.\hfill \end{array}$$

Figure 1 illustrates the numerical solution and the exact solution when $L=3$. It can be seen that the numerical solution can accurately fit the exact solution.

 Example 2.

In Equation (13), let $\widehat{T}$ and $a\left(x\right)$ be equal to 1 and μ equal to 10. Consider the following initial boundary value problem

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{}_0^c{D}_t^{\alpha }w(x,t)+10\frac{{\partial }^{4}w(x,t)}{\partial x^4}=f(x,t),x\in (0,1),t\in (0,1],\hfill \\ w(x,0)=0,x\in (0,1),\hfill \\ \frac{\partial w(x,t)}{\partial x}{{|}_{x=0}=0,\frac{\partial w(x,t)}{\partial x}|}_{x=1}=0,x\in (0,1),t\in (0,1],\hfill \\ \frac{{\partial }^{3}w(x,t)}{\partial x^3}{{|}_{x=0}=0,\frac{{\partial }^{3}w(x,t)}{\partial x^3}|}_{x=1}=0,x\in (0,1),t\in (0,1].\hfill \end{array}\right.\end{array}$$

(56)

The exact solution is $w(x,t)=t^{2}cos\left(\pi x\right)$. The spatial step is $h=1/N$ and the time step is $\tau =1/J$.

Table 5. The computational errors and convergence orders of W in space.

$\mathit{\alpha }$	h	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
0.4	1/4	$8.7432\times {10}^{-2}$	-	$4.3716\times {10}^{-2}$	-
	1/8	$2.2442\times {10}^{-2}$	1.9619	$7.9346\times {10}^{-3}$	2.4619
	1/16	$5.6459\times {10}^{-3}$	1.9909	$1.4115\times {10}^{-3}$	2.4909
	1/32	$1.4138\times {10}^{-3}$	1.9976	$2.4993\times {10}^{-4}$	2.4976
0.6	1/4	$8.4477\times {10}^{-2}$	-	$4.2238\times {10}^{-2}$	-
	1/8	$2.1739\times {10}^{-2}$	1.9583	$7.6860\times {10}^{-3}$	2.4583
	1/16	$5.4745\times {10}^{-3}$	1.9895	$1.3686\times {10}^{-3}$	2.4895
	1/32	$1.3734\times {10}^{-3}$	1.9950	$2.4279\times {10}^{-3}$	2.4950
0.8	1/4	$8.2380\times {10}^{-2}$	-	$4.1190\times {10}^{-2}$	-
	1/8	$2.1252\times {10}^{-2}$	1.9547	$7.5139\times {10}^{-3}$	2.4547
	1/16	$5.3659\times {10}^{-3}$	1.9857	$1.3415\times {10}^{-3}$	2.4857
	1/32	$1.3578\times {10}^{-3}$	1.9825	$2.4003\times {10}^{-3}$	2.4825

,

Table 6. The computational errors and convergence orders of V in space.

$\mathit{\alpha }$	h	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
0.4	1/4	$2.0975\times {10}^{-1}$	-	$1.4832\times {10}^{-1}$	-
	1/8	$5.1317\times {10}^{-2}$	2.0312	$3.6287\times {10}^{-2}$	2.0312
	1/16	$1.2760\times {10}^{-2}$	2.0078	$9.0228\times {10}^{-3}$	2.0078
	1/32	$3.1862\times {10}^{-3}$	2.0017	$2.2530\times {10}^{-3}$	2.0017
0.6	1/4	$1.9996\times {10}^{-1}$	-	$1.4139\times {10}^{-1}$	-
	1/8	$2.9079\times {10}^{-2}$	2.0265	$3.4704\times {10}^{-2}$	2.0265
	1/16	$1.2220\times {10}^{-2}$	2.0059	$8.6407\times {10}^{-3}$	2.0059
	1/32	$3.0592\times {10}^{-3}$	1.9980	$2.1632\times {10}^{-3}$	1.9980
0.8	1/4	$1.9302\times {10}^{-1}$	-	$1.3648\times {10}^{-1}$	-
	1/8	$4.7530\times {10}^{-2}$	2.0218	$3.3609\times {10}^{-2}$	2.0218
	1/16	$1.1878\times {10}^{-2}$	2.0006	$8.3987\times {10}^{-3}$	2.0006
	1/32	$3.0101\times {10}^{-3}$	1.9803	$2.1285\times {10}^{-3}$	1.9803

,

Table 7. The computational errors and convergence orders of P in space.

$\mathit{\alpha }$	h	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
0.4	1/4	$3.6117\times {10}^{-1}$	-	$2.7643\times {10}^{-1}$	-
	1/8	$9.5148\times {10}^{-2}$	1.9244	$6.8598\times {10}^{-2}$	2.0107
	1/16	$2.4090\times {10}^{-2}$	1.9817	$1.7117\times {10}^{-2}$	2.0027
	1/32	$6.0431\times {10}^{-3}$	1.9951	$4.2782\times {10}^{-3}$	2.0003
0.6	1/4	$3.3347\times {10}^{-1}$	-	$2.5523\times {10}^{-1}$	-
	1/8	$8.8296\times {10}^{-2}$	1.9171	$6.3658\times {10}^{-2}$	2.0034
	1/16	$2.2404\times {10}^{-2}$	1.9786	$1.5919\times {10}^{-2}$	1.9996
	1/32	$5.6446\times {10}^{-3}$	1.9888	$3.9962\times {10}^{-3}$	1.9940
0.8	1/4	$3.1382\times {10}^{-1}$	-	$2.4019\times {10}^{-1}$	-
	1/8	$8.3554\times {10}^{-2}$	1.9092	$6.0239\times {10}^{-2}$	1.9954
	1/16	$2.1336\times {10}^{-2}$	1.9694	$1.5160\times {10}^{-2}$	1.9905
	1/32	$5.4909\times {10}^{-3}$	1.9582	$3.8873\times {10}^{-3}$	1.9634

,

Table 8. The computational errors and convergence orders of U in space.

$\mathit{\alpha }$	h	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
0.4	1/4	$4.0755\times {10}^{-1}$	-	$2.8818\times {10}^{-1}$	-
	1/8	$1.0436\times {10}^{-1}$	1.9654	$7.3792\times {10}^{-2}$	1.9654
	1/16	$2.6240\times {10}^{-2}$	1.9917	$1.8554\times {10}^{-2}$	1.9917
	1/32	$6.5740\times {10}^{-3}$	1.9969	$4.6485\times {10}^{-3}$	1.9969
0.6	1/4	$3.1577\times {10}^{-1}$	-	$2.2328\times {10}^{-1}$	-
	1/8	$8.2550\times {10}^{-2}$	1.9355	$5.8372\times {10}^{-2}$	1.9355
	1/16	$2.0924\times {10}^{-2}$	1.9801	$1.4796\times {10}^{-2}$	1.9801
	1/32	$5.3212\times {10}^{-3}$	1.9754	$3.7626\times {10}^{-3}$	1.9754
0.8	1/4	$2.5066\times {10}^{-1}$	-	$1.7724\times {10}^{-1}$	-
	1/8	$6.7458\times {10}^{-2}$	1.8937	$4.7700\times {10}^{-2}$	1.8937
	1/16	$1.7557\times {10}^{-2}$	1.9419	$1.2415\times {10}^{-2}$	1.9419
	1/32	$4.8378\times {10}^{-3}$	1.8596	$3.4208\times {10}^{-3}$	1.8596

,

Table 9. The computational errors and convergence orders of W in time.

$\mathit{\alpha }$	$\mathit{\tau }$	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
0.4	1/8	$1.1526\times {10}^{-3}$	-	$3.6020\times {10}^{-5}$	-
	1/16	$3.9529\times {10}^{-4}$	1.5440	$1.2353\times {10}^{-5}$	1.5440
	1/32	$1.3445\times {10}^{-4}$	1.5558	$4.2017\times {10}^{-6}$	1.5558
	1/64	$4.4648\times {10}^{-5}$	1.5904	$1.3953\times {10}^{-6}$	1.5904
0.6	1/8	$2.7928\times {10}^{-3}$	-	$8.7275\times {10}^{-5}$	-
	1/16	$1.0864\times {10}^{-3}$	1.3621	$3.3950\times {10}^{-5}$	1.3621
	1/32	$4.1776\times {10}^{-4}$	1.3788	$1.3055\times {10}^{-5}$	1.3788
	1/64	$1.5565\times {10}^{-4}$	1.4244	$4.8640\times {10}^{-6}$	1.4244
0.8	1/8	$6.1370\times {10}^{-3}$	-	$1.9178\times {10}^{-4}$	-
	1/16	$2.7131\times {10}^{-3}$	1.1776	$8.4783\times {10}^{-5}$	1.1776
	1/32	$1.1895\times {10}^{-3}$	1.1896	$3.7171\times {10}^{-5}$	1.1895
	1/64	$5.2209\times {10}^{-4}$	1.1879	$1.6315\times {10}^{-5}$	1.1879

,

Table 10. The computational errors and convergence orders of V in time.

$\mathit{\alpha }$	$\mathit{\tau }$	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
0.4	1/8	$3.6206\times {10}^{-3}$	-	$2.5601\times {10}^{-3}$	-
	1/16	$1.2416\times {10}^{-3}$	1.5440	$8.7793\times {10}^{-4}$	1.5440
	1/32	$4.2024\times {10}^{-4}$	1.5629	$2.9716\times {10}^{-4}$	1.5629
	1/64	$1.3963\times {10}^{-4}$	1.5896	$9.8734\times {10}^{-5}$	1.5896
0.6	1/8	$8.7728\times {10}^{-3}$	-	$6.2033\times {10}^{-3}$	-
	1/16	$3.4129\times {10}^{-3}$	1.3620	$2.4133\times {10}^{-3}$	1.3620
	1/32	$1.3117\times {10}^{-3}$	1.3795	$9.2756\times {10}^{-4}$	1.3795
	1/64	$4.8812\times {10}^{-4}$	1.4262	$3.4515\times {10}^{-4}$	1.4262
0.8	1/8	$1.9279\times {10}^{-2}$	-	$1.3632\times {10}^{-2}$	-
	1/16	$8.5226\times {10}^{-3}$	1.1776	$6.0264\times {10}^{-3}$	1.1776
	1/32	$3.7346\times {10}^{-3}$	1.1903	$2.6408\times {10}^{-3}$	1.1903
	1/64	$1.6400\times {10}^{-3}$	1.1873	$1.1596\times {10}^{-3}$	1.1873

,

Table 11. The computational errors and convergence orders of P in time.

$\mathit{\alpha }$	$\mathit{\tau }$	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
0.4	1/8	$1.1372\times {10}^{-2}$	-	$8.0408\times {10}^{-3}$	-
	1/16	$3.8972\times {10}^{-3}$	1.5449	$2.7557\times {10}^{-3}$	1.5449
	1/32	$1.3166\times {10}^{-3}$	1.5656	$9.3098\times {10}^{-4}$	1.5656
	1/64	$4.3497\times {10}^{-4}$	1.5979	$3.0753\times {10}^{-4}$	1.5980
0.6	1/8	$2.7558\times {10}^{-2}$	-	$1.9486\times {10}^{-2}$	-
	1/16	$1.0719\times {10}^{-2}$	1.3623	$7.5792\times {10}^{-3}$	1.3623
	1/32	$4.1174\times {10}^{-3}$	1.3804	$2.9114\times {10}^{-3}$	1.3803
	1/64	$1.5301\times {10}^{-3}$	1.4281	$1.0817\times {10}^{-3}$	1.4285
0.8	1/8	$6.0563\times {10}^{-2}$	-	$4.2825\times {10}^{-2}$	-
	1/16	$2.6771\times {10}^{-2}$	1.1778	$1.8930\times {10}^{-2}$	1.1778
	1/32	$1.1729\times {10}^{-2}$	1.1906	$8.2936\times {10}^{-3}$	1.1906
	1/64	$5.1483\times {10}^{-3}$	1.1879	$3.6405\times {10}^{-3}$	1.1879

and

Table 12. The computational errors and convergence orders of U in time.

$\mathit{\alpha }$	$\mathit{\tau }$	${\mathit{l}}^{\mathit{\infty }}-\mathit{Error}$	$\mathbf{Order}$	${\mathit{l}}^2-\mathit{Error}$	$\mathbf{Order}$
0.4	1/8	$3.5715\times {10}^{-2}$	-	$2.5254\times {10}^{-2}$	-
	1/16	$1.2232\times {10}^{-2}$	1.5458	$8.6497\times {10}^{-3}$	1.5458
	1/32	$4.1246\times {10}^{-3}$	1.5684	$2.9167\times {10}^{-3}$	1.5683
	1/64	$1.3548\times {10}^{-3}$	1.6061	$9.5780\times {10}^{-4}$	1.6065
0.6	1/8	$8.6565\times {10}^{-2}$	-	$6.1211\times {10}^{-2}$	-
	1/16	$3.3662\times {10}^{-2}$	1.3626	$2.3803\times {10}^{-2}$	1.3626
	1/32	$1.2924\times {10}^{-2}$	1.3811	$9.1385\times {10}^{-3}$	1.3811
	1/64	$4.7941\times {10}^{-3}$	1.4308	$3.3898\times {10}^{-3}$	1.4307
0.8	1/8	$1.9026\times {10}^{-1}$	-	$1.3453\times {10}^{-1}$	-
	1/16	$8.4093\times {10}^{-2}$	1.1779	$5.9463\times {10}^{-2}$	1.1779
	1/32	$3.6837\times {10}^{-2}$	1.1908	$2.6047\times {10}^{-2}$	1.1908
	1/64	$1.6163\times {10}^{-2}$	1.1885	$1.1429\times {10}^{-2}$	1.1885

show the time and space errors, as well as the convergence orders of W, V, P, and U, when α is $0.4,0.6$, and $0.8$. Table 5, Table 6, Table 7 and Table 8 show the maximum norm errors, discrete L2 norm errors, and convergence orders of W, V, P, and U as the mesh size h is reduced with a fixed $\tau =1/2^{10}$. Table 9, Table 10, Table 11 and Table 12 show the maximum norm errors, discrete L2 norm errors, and convergence orders of W, V, P, and U as the mesh size τ is reduced with a fixed $h=1/2^{10}$. It can be seen from the tables that for the maximum norm and discrete L2 norm errors, the space convergence order has reached the second order and the time convergence order has reached the $2-\alpha$ order. Therefore, the validity of the block-centered finite-difference scheme is verified.

Figure 2 and Figure 3 show time and space images of the numerical solution W and the analytic solution w when the mesh is divided into $N=J^2=2^8$. Figure 4, Figure 5 and Figure 6 illustrate comparisons of the numerical solution W and the analytic solution w when α is $0.4,0.6$, and $0.8$, respectively. It can be seen that the numerical solution W closely matches the analytic solution w.

6. Conclusions

This paper discusses the block-centered finite-difference method for solving fourth-order parabolic equations of fractional-order time derivatives. By introducing intermediate variables, the fourth-order differential equation is transformed into a system of first-order differential equations. Based on the equivalence of the block-centered finite-difference scheme and the mixed finite-element method with a special numerical quadrature formula, the stability and convergence of the block-centered finite-difference scheme have been proved. The effectiveness of the block-centered finite-difference scheme is verified through numerical examples. The block-centered difference method can be extended to other fourth-order differential equations with Neumann boundary conditions. In future work, we will study the solutions of other types of fourth-order partial differential equations using the block-centered finite-difference method and numerical examples on non-uniform grids.