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

Preprint

Article

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

Altmetrics

Downloads

143

Views

Comments

A peer-reviewed article of this preprint also exists.

Taixiu Zhang,Zhe Yin,

Ailing Zhu^*

Taixiu Zhang,Zhe Yin,

Ailing Zhu^*

This version is not peer-reviewed

Submitted:

17 May 2023

Posted:

18 May 2023

You are already at the latest version

Alerts

Abstract

In this paper, we consider the fourth-order parabolic equations with integer and fractional order time derivatives with Neumann boundary conditions. The integer order time derivatives are approximated by backward Euler difference quotients, and the fractional order time derivatives are approximated by L1 interpolation. We propose the block-centered finite difference scheme for fourth-order parabolic equations of integer and 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 given to verify the effectiveness of the block-centered finite difference scheme.

Keywords:

Subject: Computer Science and Mathematics - Computational Mathematics

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, hence maintaining 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] discussed 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 and proved that the errors of the discrete L2 norms for the solution and the first derivative are both of the second order. Rui and Pan considered in [3,4,5], respectively, the block-centered finite-difference method for the nonlinear Darcy-Forchheimer equation and parabolic equations with time-dependent coefficient. The block-centered finite-difference method on non-uniform grids was 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. Xu [13] proposed two block-centered finite difference schemes for parabolic problems. [14] considered the block-centered finite difference method for a class of parabolic integral-differential equations with Neumann boundary conditions in a rectangular region. Dupont [15] provided a domain decomposition method for solving parabolic equations in the block-centered difference scheme. Li and Rui [16] introduced and analyzed the block-centered finite difference method for distributed-order time-fractional diffusion-wave equations with Neumann boundary conditions. In addition, [17,18,19,20,21] discussed the two-grid and the 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 - α)

-order-accuracy in time, and the unconditional stability and convergence have been proven theoretically. In [22], 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, demonstrating the stability of the solution and flux and performed the 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 [23], such as the competition and spatial sorting of biological populations, the migration of riverbeds, the charge density distribution of quantum semiconductors, etc [24,25]. In recent years, the numerical method of the fourth-order parabolic equation has attracted extensive attention of researchers. In [25], 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. Kaur [26] gave a three-level implicit compact difference scheme for the generalized form of the fourth-order parabolic equation. [27,28,29] discussed the finite element method for the initial boundary value problems of fourth-order parabolic equations. 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. Currently, many researchers are dedicated to the study of fractional order differential equations. Aziz [30] studied two inverse source problems of fourth-order parabolic equations with fractional time derivatives. Li and Liao [31] used a class of

L 1

-Galerkin finite element methods to study the numerical solution of time-fractional nonlinear parabolic problems and gave the optimal error estimations of several fully discrete linearized Galerkin finite element methods for solving nonlinear problems. [32] 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 [33], Liu and Du proposed and discussed the finite difference/finite element method for solving non-linear time-fractional fourth-order reaction and diffusion problems. A new implicit compact difference scheme for fourth-order fractional diffuse wave systems was constructed in [34]. In addition, Ji and Sun [35] 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 have been present in the literature. Therefore, it is of great theoretical and practical significance to propose and develop a block-centered finite difference method for fourth-order parabolic equations. This paper discusses the block-centered finite difference method [34] for fourth-order parabolic equations of integer and fractional order time derivatives with Neumann boundary conditions. The integer order time derivatives are approximated by backward Euler difference quotients, and the fractional order time derivatives are approximated by

L 1

interpolation. The block-centered finite difference schemes for the two types of equations are established, and the error estimations of the discrete

L 2

norm of the approximate solution and its derivatives are given. Numerical examples are provided to verify the effectiveness of the block-centered finite difference scheme.

This paper is structured as follows: Section 2 introduces the notations used in this paper. Section 3 gives the block-centered finite difference scheme and error estimation for fourth-order ordinary differential equations. Section 4 and 5, respectively, establish the block-centered finite difference schemes for the fourth-order parabolic equations of integer and fractional order time derivatives and prove the stability and convergence of the schemes. In Section 6, numerical examples are given to verify the convergence of the proposed schemes.

2. Notation

We first introduce some notations and definitions used in this paper, which will help the following analysis. At the same time, we use notations similar to those in [2]. Define the partition

T_{x}

Ω = (x_{L}, x_{R})

T_{x} : x_{L} = x_{1 / 2} < x_{3 / 2} < \dots < x_{N - 1 / 2} < x_{N + 1 / 2} = x_{R} .

For each

i = 1

to N, define

\begin{matrix} \begin{matrix} x_{i} = \frac{1}{2} (x_{i + 1 / 2} + x_{i - 1 / 2}), \\ h_{i} = x_{i + 1 / 2} - x_{i - 1 / 2}, \\ h = max_{i} h_{i}, \\ h_{i + 1 / 2} = x_{i + 1} - x_{i} = \frac{1}{2} (h_{i} + h_{i + 1}), \\ Ω_{i} = (x_{i - 1 / 2}, x_{i + 1 / 2}) . \end{matrix} \end{matrix}

The block-centered dual partition grids are defined as

T_{x}^{*} = {x_{i}}

Take positive integers M, let

τ = \hat{T} / J

t_{n} = n τ

(0 \leq n \leq J)

For any function

g (x)

, let

g_{i}

g_{i + 1 / 2}

g_{i}^{n}

donate

g (x_{i})

g (x_{i + 1 / 2})

g (x_{i}, t_{n})

Define the following notations

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

{[D_{x} g]}_{i} = \frac{g_{i + 1 / 2} - g_{i - 1 / 2}}{h_{i}},

d_{t} g_{i}^{n} = \frac{g_{i}^{n} - g_{i}^{n - 1}}{τ} .

For functions F and G, define the midpoint quadrature formula and trapezoidal quadrature formula on

Ω_{i}

\begin{matrix} {(F, G)}_{M, Ω_{i}} = h_{i} F (x_{i}) G (x_{i}), \end{matrix}

and

\begin{matrix} {(F, G)}_{T, Ω_{i}} = \frac{h_{i}}{2} [F (x_{i - 1 / 2}) G (x_{i - 1 / 2}) + F (x_{i + 1 / 2}) G (x_{i + 1 / 2})] . \end{matrix}

Given functions

f (x)

and

g (x)

, define

L 2

inner product and norm

\begin{matrix} (f, g) = \int_{Ω} f (x) g (x) d x, \\ {∥ f ∥}^{2} = (f, f), \end{matrix}

and the discrete inner products and norms

\begin{matrix} {(f, g)}_{M} = \sum_{i = 1}^{N} h_{i} f_{i} g_{i}, \\ {∥ f ∥}_{M}^{2} = {(f, f)}_{M}, \\ {(f, g)}_{T} = \sum_{i = 1}^{N} h_{i + 1 / 2} f_{i + 1 / 2} g_{i + 1 / 2}, \\ {∥ f ∥}_{T}^{2} = {(f, f)}_{T} . \end{matrix}

Define

S_{c}^{d} (T_{x})

to be the finite-dimensional space of one-dimensional functions which have c continuous derivatives on

(x_{L}, x_{R})

, and which are piecewise polynomials of degree d in each interval

Ω_{i}

. When

c = - 1

, the functions themselves may discontinuous.

The notation

∥ g ∥ = O (h^{k}), k > 0

means there exists a constant C such that

∥ g (x) ∥ \leq C h^{k},

as h approaches zero.

3. Fourth-order ordinary differential equation

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

We consider fourth-order ordinary differential equations with Neumann boundary condition

\begin{matrix} \{\begin{matrix} {(a (x) w^{″})}^{″} = f (x), x \in (x_{L}, x_{R}), \\ w^{'} {|_{x = x_{L}} = 0, w^{'} |}_{x = x_{R}} = 0, \\ {(a (x) w^{″})}^{'} {|_{x = x_{L}} = 0, {(a (x) w^{″})}^{'} |}_{x = x_{R}} = 0 . \end{matrix} \end{matrix}

(1)

Where

f (x)

is a known smooth function.

Let

v (x) = - w^{'} (x), p (x) = a (x) v^{'} (x), u (x) = - p^{'} (x),

then, we have

\begin{matrix} u^{'} (x) = f (x), & x \in (x_{L}, x_{R}), \\ u (x) = & - p^{'} (x), \\ a (x) v^{'} (x) = & p (x), \\ v (x) = & - w^{'} (x), \\ u (x) = 0, v (x) = 0, & a t x = x_{L}, x = x_{R} . \end{matrix}

(2)

The block-centered finite difference approximations

U_{i + 1 / 2}

V_{i + 1 / 2}

W_{i}

and

P_{i}

u (x_{i + 1 / 2})

v (x_{i + 1 / 2})

w (x_{i})

and

p (x_{i}) (i = 1, 2, \dots, N)

, respectively, are satisfied

\begin{matrix} {[D_{x} U]}_{i} = f_{i}, & U_{1 / 2} = 0, U_{N + 1 / 2} = 0, \end{matrix}

(3)

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

(4)

\begin{matrix} {[D_{x} V]}_{i} = \frac{P_{i}}{a_{i}}, & V_{1 / 2} = 0, V_{N + 1 / 2} = 0, \end{matrix}

(5)

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

(6)

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

\begin{matrix} {({(U)}^{'}, 1)}_{Ω_{i}} = {(f, 1)}_{M, Ω_{i}}, \end{matrix}

(7)

\begin{matrix} {(U, χ)}_{T} - (P, χ^{'}) = 0, χ \in \tilde{S}, \end{matrix}

(8)

\begin{matrix} {({(V)}^{'}, 1)}_{Ω_{i}} = {(\frac{P}{a}, 1)}_{M, Ω_{i}}, \end{matrix}

(9)

\begin{matrix} {(V, χ)}_{T} - (W, χ^{'}) = 0, χ \in \tilde{S} . \end{matrix}

(10)

Where

\tilde{S} = S_{0}^{1} (T_{x}) ⋂ {χ : χ (x_{L}) = χ (x_{R}) = 0}

, U and V are in

\tilde{S}

, and P and W are in

S_{- 1}^{0} (T_{x})

According to the reference [2], through Taylor’s expansion, the following conclusions can be obtained.

Lemma 1

([2]). If

w^{(5)} (x)

is continuous and

f^{″} (x)

is in

L^{1} (Ω_{i})

for all i,

\begin{matrix} {∥ U - u ∥}_{T} = O (h^{2}), {∥ P - p ∥}_{M} = O (h^{2}), \\ {∥ V - v ∥}_{T} = O (h^{2}), {∥ W - w ∥}_{M} = O (h^{2}) . \end{matrix}

(11)

The second-order error estimate of block-centered difference scheme for fourth-order ordinary differential equations is obtained.

4. Time integer fourth-order parabolic equation

In this section, we discuss the block-centered finite difference scheme for fourth-order parabolic equations of integer order, and prove the stability and convergence of its discrete scheme.

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

\begin{matrix} \{\begin{matrix} w_{t} (x, t) + μ \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, \hat{T}], \\ w (x, 0) = φ (x), x \in (x_{L}, x_{R}), \\ \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, \hat{T}], \\ \frac{\partial}{\partial x} (a (x) \frac{\partial^{2} w (x, t)}{\partial x^{2}}) {|_{x = x_{L}} = 0, \frac{\partial}{\partial x} (a (x) \frac{\partial^{2} w (x, t)}{\partial x^{2}}) |}_{x = x_{R}} = 0, t \in (0, \hat{T}] . \end{matrix} \end{matrix}

(12)

Where

μ

is a constant,

a (x)

f (x, t)

and

φ (x)

are known smooth functions and assume

0 < a_{0} \leq a (x) \leq a_{1}

4.1. Block-centered finite difference scheme

In this subsection, we give the block-centered difference scheme for the fourth-order parabolic equation.

In order to discuss the errors of fourth-order parabolic partial differential equations, the following lemma is introduced (see reference [2]).

Lemma 2

([2]). Let

ϑ_{i + 1 / 2}, q_{i}

and

ω_{i + 1 / 2}

be any values such that

ϑ_{1 / 2} = ϑ_{N + 1 / 2} = ω_{1 / 2} = ω_{N + 1 / 2} = 0,

ϑ_{i + 1 / 2} = - {[d_{x} q]}_{i + 1 / 2},

then

\begin{matrix} {(ϑ, ω)}_{T} = {(q, D_{x} ω)}_{M} . \end{matrix}

(13)

The block-centered finite difference method of (12) defines

{U_{i + 1 / 2}^{n}}_{n = 1}^{M}

{V_{i + 1 / 2}^{n}}_{n = 1}^{M}

{W_{i}^{n}}_{n = 1}^{M}

and

{P_{i}^{n}}_{n = 1}^{M} (i = 1, 2, \dots, N)

, satisfying

\begin{matrix} {[d_{t} W]}_{i}^{n} + μ {[D_{x} U]}_{i}^{n} = f_{i}^{n}, & U_{1 / 2}^{n} = 0, U_{N + 1 / 2}^{n} = 0, \end{matrix}

(14)

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

(15)

\begin{matrix} {[D_{x} V]}_{i}^{n} = \frac{P_{i}^{n}}{a_{i}}, & V_{1 / 2}^{n} = 0, V_{N + 1 / 2}^{n} = 0, \end{matrix}

(16)

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

(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{matrix} {(d_{t} W^{n}, 1)}_{Ω_{i}} + μ {({(U^{n})}^{'}, 1)}_{Ω_{i}} = {(f^{n}, 1)}_{M, Ω_{i}}, \end{matrix}

(18)

\begin{matrix} {(U^{n}, χ)}_{T} - (P^{n}, χ^{'}) = 0, χ \in \tilde{S}, \end{matrix}

(19)

\begin{matrix} {({(V^{n})}^{'}, 1)}_{Ω_{i}} = {(\frac{P^{n}}{a}, 1)}_{M, Ω_{i}}, \end{matrix}

(20)

\begin{matrix} {(V^{n}, χ)}_{T} - (W^{n}, χ^{'}) = 0, χ \in \tilde{S} . \end{matrix}

(21)

4.2. Stability Analysis

In this subsection, we give the stability analysis of the scheme

(14) - (17)

Theorem 1

For the block-centered difference scheme

(14) - (17)

, the following stable inequality holds unconditionally, for sufficiently small τ

∥ W^{n} ∥_{M} \leq C_{1} ∥ W^{0} ∥_{M} + C_{2} max_{1 \leq k \leq J} {∥ f^{k} ∥}_{M} .

Proof .

For

n = 1

, multiplying (14) by

h_{i} W_{i}^{1}

and summing on i from 1 to N, we have

{(d_{t} W^{1}, W^{1})}_{M} + μ {(W_{x x}^{1}, W_{x x}^{1})}_{M} = {(f^{1}, W^{1})}_{M}

. So,

\frac{1}{τ} {(W^{1}, W^{1})}_{M} \leq \frac{1}{τ} {(W^{0}, W^{1})}_{M} + {(f^{1}, W^{1})}_{M} .

Using the Cauchy-Schwarz inequality, we can get

{∥ W^{n} ∥}_{M}^{2} \leq \frac{1}{1 - τ} {∥ W^{0} ∥}_{M}^{2} + \frac{τ}{1 - τ} {∥ f^{1} ∥}_{M}^{2} .

So, we have

∥ W^{1} ∥_{M} \leq C_{1} ∥ W^{0} ∥_{M} + C_{2} {∥ f^{1} ∥}_{M} .

For

n \geq 2

, we suppose that the stability conclusion of the difference scheme is valid when

k \leq n - 1

Multiplying (14) by

h_{i} W_{i}^{n}

and summing on i from 1 to N, we can get

{(d_{t} W^{n}, W^{n})}_{M} \leq {(f^{n}, W^{n})}_{M}

, so

\frac{1}{τ} {(W^{n}, W^{n})}_{M} \leq \frac{1}{τ} {(W^{n - 1}, W^{n})}_{M} + {(f^{n}, W^{n})}_{M} .

By the Cauchy-Schwarz inequality, we deduce that

∥ W^{n} ∥_{M}^{2} \leq \frac{1}{1 - τ} ∥ W^{n - 1} ∥_{M}^{2} + \frac{τ}{1 - τ} {∥ f^{n} ∥}_{M}^{2} .

Using mathematical induction, we can get

\begin{matrix} ∥ W^{n} ∥_{M}^{2} & \leq \frac{1}{1 - τ} [(C_{1}^{'} ∥ W^{0} ∥_{M}^{2} + C_{2}^{'} max_{1 \leq k \leq J} ∥ f^{k} ∥_{M}^{2}) + τ ∥ f^{n} ∥_{M}^{2}] \\ \leq \frac{1}{1 - τ} [C_{1}^{'} ∥ W^{0} ∥_{M}^{2} + (C_{2}^{'} + τ) max_{1 \leq k \leq J} ∥ f^{k} ∥_{M}^{2}] . \end{matrix}

There are constants

C_{3}

and

C_{4}

, such that

∥ W^{n} ∥_{M}^{2} \leq C_{3} ∥ W^{0} ∥_{M}^{2} + C_{4} max_{1 \leq k \leq J} {∥ f^{k} ∥}_{M}^{2} .

Thus, we have

{∥ W^{n} ∥}_{M}^{2} \leq C_{1} ∥ W^{0} ∥_{M} + C_{2} max_{1 \leq k \leq J} {∥ f^{k} ∥}_{M} .

We complete the proof of stability analysis. □

4.3. Error Analysis

In this subsection, we give the error analysis of the block-centered finite difference scheme of the fourth-order parabolic equation.

Error estimates for this finite difference scheme of

(14) - (17)

are derived using a technique of mixed finite element methods for parabolic partial differential equations. We just have to compare

U^{n}

V^{n}

W^{n}

and

P^{n}

to the solution

E_{h} u^{n}, E_{h} v^{n}, E_{h} w^{n}

and

E_{h} p^{n}

which are block-centered finite difference approximations to the fourth-order ordinary differential equations with Neumann boundary conditions.

\begin{matrix} μ \frac{\partial^{2}}{\partial x^{2}} (a (x) \frac{\partial^{2} w^{n}}{\partial x^{2}}) = ϕ = f^{n} - w_{t}^{n}, & x \in (x_{L}, x_{R}) \end{matrix}

(22)

\begin{matrix} \frac{\partial w^{n}}{\partial x} {|_{x = x_{L}} = 0, \frac{\partial w^{n}}{\partial x} |}_{x = x_{R}} = 0, \end{matrix}

(23)

\begin{matrix} \frac{\partial}{\partial x} (a (x) \frac{\partial^{2} w^{n}}{\partial x^{2}}) {|_{x = x_{L}} = 0, \frac{\partial}{\partial x} (a (x) \frac{\partial^{2} w^{n}}{\partial x^{2}}) |}_{x = x_{R}} = 0 . \end{matrix}

(24)

For fixed n,

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}

are defined as

\begin{matrix} μ {[D_{x} E_{h} u]}_{i}^{n} = f_{i}^{n} - {(w_{t})}_{i}^{n}, & E_{h} u_{1 / 2}^{n} = 0, E_{h} u_{N + 1 / 2}^{n} = 0, \end{matrix}

(25)

\begin{matrix} E_{h} u_{i + 1 / 2}^{n} = {[- d_{x} E_{h} p]}_{i + 1 / 2}^{n}, \end{matrix}

(26)

\begin{matrix} {[D_{x} E_{h} v]}_{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, \end{matrix}

(27)

\begin{matrix} E_{h} v_{i + 1 / 2}^{n} = {[- d_{x} E_{h} w]}_{i + 1 / 2}^{n}, \end{matrix}

(28)

where

E_{h} u^{0} = {(a (x) φ^{''} (x))}^{'}

E_{h} p^{0} = - a (x) φ^{''} (x)

E_{h} v^{0} = - φ^{'} (x)

and

E_{h} w^{0} = φ (x)

Equations

(25) - (28)

can be written as a mixed finite element method with approximate integration

\begin{matrix} μ {({(E_{h} u^{n})}^{'}, 1)}_{Ω_{i}} = {(f^{n}, 1)}_{M, Ω_{i}} - {(w_{t}^{n}, 1)}_{M, Ω_{i}} \end{matrix}

(29)

\begin{matrix} {(E_{h} u^{n}, χ)}_{T} - (E_{h} p^{n}, χ^{'}) = 0, χ \in \tilde{S}, \end{matrix}

(30)

\begin{matrix} {({(E_{h} v^{n})}^{'}, 1)}_{Ω_{i}} = {(\frac{E_{h} p^{n}}{a}, 1)}_{M, Ω_{i}}, \end{matrix}

(31)

\begin{matrix} {(E_{h} v^{n}, χ)}_{T} - (E_{h} w^{n}, χ^{'}) = 0, χ \in \tilde{S} . \end{matrix}

(32)

By the error of the ellipse projection, we have

\begin{matrix} ∥ E_{h} u^{n} - u^{n} ∥_{T} = O (h^{2}), {∥ E_{h} p^{n} - p^{n} ∥}_{M} = O (h^{2}), \\ ∥ E_{h} v^{n} - v^{n} ∥_{T} = O (h^{2}), {∥ E_{h} w^{n} - w^{n} ∥}_{M} = O (h^{2}), \end{matrix}

(33)

hold for w sufficiently smooth.

Differentiating t in equations

(25) - (28)

, we can get the following estimation

\begin{matrix} ∥ E_{h} w_{t}^{n} - w_{t}^{n} ∥_{M} = O (h^{2}) . \end{matrix}

(34)

Set

ξ_{i + 1 / 2}^{n} = U_{i + 1 / 2}^{n} - E_{h} u_{i + 1 / 2}^{n}

η_{i}^{n} = P_{i}^{n} - E_{h} p_{i}^{n}

θ_{i + 1 / 2}^{n} = V_{i + 1 / 2}^{n} - E_{h} v_{i + 1 / 2}^{n}

ζ_{i}^{n} = W_{i}^{n} - E_{h} w_{i}^{n}

σ_{i}^{n} = E_{h} w_{i}^{n} - w_{i}^{n}

and

ϵ_{i}^{n} = d_{t} w_{i}^{n} - {(w_{t})}_{i}^{n}

Subtracting (25) from (14), we get

\begin{matrix} d_{t} ζ_{i}^{n} + μ {[D_{x} ξ]}_{i}^{n} = - d_{t} σ_{i}^{n} - ϵ_{i}^{n}, & ξ_{1 / 2}^{n} = 0, ξ_{N + 1 / 2}^{n} = 0 . \end{matrix}

(35)

Subtracting (26), (27) and (28) from (15), (16) and (17), correspondingly, we obtain

\begin{matrix} ξ_{i + 1 / 2}^{n} = {[- d_{x} η]}_{i + 1 / 2}^{n}, \end{matrix}

(36)

\begin{matrix} η_{i}^{n} = a_{i} {[D_{x} θ]}_{i}^{n}, & θ_{1 / 2}^{n} = 0, θ_{N + 1 / 2}^{n} = 0, \end{matrix}

(37)

\begin{matrix} θ_{i + 1 / 2}^{n} = {[- d_{x} ζ]}_{i + 1 / 2}^{n} . \end{matrix}

(38)

Multiplying (35) by

h_{i} ζ_{i}^{n}

and summing on i from 1 to N, we deduce that

\begin{matrix} {(d_{t} ζ^{n}, ζ^{n})}_{M} + μ {(D_{x} ξ^{n}, ζ^{n})}_{M} = {(- d_{t} σ_{i}^{n} - ϵ_{i}^{n}, ζ^{n})}_{M} . \end{matrix}

(39)

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

\begin{matrix} {(D_{x} ξ^{n}, ζ^{n})}_{M} & = {(D_{x} (- d_{x} η^{n}), ζ^{n})}_{M} \\ = {(d_{x} η^{n}, d_{x} ζ^{n})}_{T} \\ = - {(d_{x} η^{n}, θ^{n})}_{T} \\ = {(η^{n}, D_{x} θ^{n})}_{M} \\ = \frac{1}{a} {(η^{n}, η^{n})}_{M} . \end{matrix}

(40)

Now

\begin{matrix} {(d_{t} ζ^{n}, ζ^{n})}_{M} \geq \frac{∥ ζ^{n} ∥_{M}^{2} - {∥ ζ^{n - 1} ∥}_{M}^{2}}{2 τ} . \end{matrix}

(41)

Then (39) can be written as

\begin{matrix} (∥ ζ^{n} ∥_{M}^{2} - ∥ ζ^{n - 1} ∥_{M}^{2}) + \frac{2}{a_{1}} τ μ {(η^{n}, η^{n})}_{M} \leq 2 τ | {(- d_{t} σ^{n} - ϵ^{n}, ζ^{n})}_{M} | . \end{matrix}

(42)

By Young inequality

\begin{matrix} | {(- d_{t} σ^{n} - ϵ^{n}, ζ^{n})}_{M} | \leq ∥ d_{t} σ^{n} ∥_{M}^{2} + ∥ ϵ^{n} ∥_{M}^{2} + \frac{1}{2} {∥ ζ^{n} ∥}_{M}^{2} . \end{matrix}

(43)

Then by

(34)

, we can derive

\begin{matrix} ∥ d_{t} σ^{n} ∥_{M}^{2} & = ∥ \frac{1}{τ} (σ^{n} - σ^{n - 1}) ∥_{M}^{2} \\ = ∥ \frac{1}{τ} (\int_{t^{n - 1}}^{t^{n}} σ_{t} d t) ∥_{M}^{2} \\ = \frac{1}{τ^{2}} {∥ \int_{t^{n - 1}}^{t^{n}} σ_{t} d t ∥}_{M}^{2} \\ \leq \frac{1}{τ} (\int_{t^{n - 1}}^{t^{n}} ∥ σ_{t} ∥_{M}^{2} d t) \\ \leq C h^{4} . \end{matrix}

(44)

Using

| ϵ^{n} | = | d_{t} w^{n} - w_{t}^{n} | \leq C τ .

We can obtain

\begin{matrix} ∥ ϵ^{n} ∥_{M}^{2} \leq C τ^{2} . \end{matrix}

(45)

Substituting (44) and (45) into (43), we have

\begin{matrix} | {(- d_{t} σ^{n} - ϵ^{n}, ζ^{n})}_{M} | \leq C (τ^{2} + h^{4}) + \frac{1}{2} {∥ ζ^{n} ∥}_{M}^{2} . \end{matrix}

(46)

Then substituting (46) into (42) and summing the result on n,

n = 1, 2, \dots J

, using

ζ^{0} = 0

, we get

\begin{matrix} ∥ ζ^{J} ∥_{M}^{2} + \frac{2}{a_{1}} τ μ \sum_{n = 1}^{J} ∥ η^{n} ∥_{M}^{2} \leq C (τ^{2} + h^{4}) + τ \sum_{n = 1}^{J} {∥ ζ^{n} ∥}_{M}^{2} . \end{matrix}

(47)

By discrete Gronwall inequality, we can obtain

\begin{matrix} (1 - τ) ∥ ζ^{J} ∥_{M}^{2} + \frac{2}{a_{1}} τ μ \sum_{n = 1}^{J} {∥ η^{n} ∥}_{M}^{2} \leq C (τ^{2} + h^{4}) . \end{matrix}

(48)

When

τ

is sufficiently small such that

τ < 1

, we have

\begin{matrix} ∥ ζ^{n} ∥_{M}^{2} \leq C (τ^{2} + h^{4}), \end{matrix}

(49)

\begin{matrix} ∥ η^{n} ∥_{M}^{2} \leq C (τ^{2} + h^{4}) . \end{matrix}

(50)

Multiplying (35) by

h_{i} η_{i}^{n}

, and summing on i from 1 to N, we deduce that

\begin{matrix} {(d_{t} ζ^{n}, η^{n})}_{M} + μ {(D_{x} ξ^{n}, η^{n})}_{M} = {(- d_{t} σ_{i}^{n} - ϵ_{i}^{n}, η^{n})}_{M} . \end{matrix}

(51)

Using (36), (37), (38) and Lemma 2, we derive that

\begin{matrix} a {(d_{t} θ^{n}, θ^{n})}_{T} + μ {(ξ^{n}, ξ^{n})}_{T} = {(- d_{t} σ^{n} - ϵ^{n}, η^{n})}_{M} . \end{matrix}

(52)

Now

\begin{matrix} {(d_{t} θ^{n}, θ^{n})}_{T} \geq \frac{∥ θ^{n} ∥_{T}^{2} - {∥ θ^{n - 1} ∥}_{T}^{2}}{2 τ} . \end{matrix}

(53)

Then (52) can be deformable to

\begin{matrix} a_{0} (∥ θ^{n} ∥_{T}^{2} - ∥ θ^{n - 1} ∥_{T}^{2}) + 2 τ μ ∥ ξ^{n} ∥_{T}^{2} \leq 2 τ | {(- d_{t} σ^{n} - ϵ^{n}, η^{n})}_{M} | . \end{matrix}

(54)

By Young inequality

\begin{matrix} | {(- d_{t} σ^{n} - ϵ^{n}, η^{n})}_{M} | \leq ∥ d_{t} σ^{n} ∥_{M}^{2} + ∥ ϵ^{n} ∥_{M}^{2} + \frac{1}{2} {∥ η^{n} ∥}_{M}^{2} . \end{matrix}

(55)

Using (44), (45), (50) and (55), we can get

\begin{matrix} | {(- d_{t} σ^{n} - ϵ^{n}, η^{n})}_{M} | \leq C (τ^{2} + h^{4}) . \end{matrix}

(56)

Then substituting (56) into (54), we obtain

\begin{matrix} a_{0} (∥ θ^{n} ∥_{T}^{2} - ∥ θ^{n - 1} ∥_{T}^{2}) + 2 τ μ ∥ ξ^{n} ∥_{T}^{2} \leq τ C (τ^{2} + h^{4}) . \end{matrix}

(57)

Summing (57) on n,

n = 1, 2, \dots J

, using

θ^{0} = 0

, we obtain

\begin{matrix} a_{0} ∥ θ^{J} ∥_{T}^{2} + 2 τ μ \sum_{n = 1}^{J} {∥ ξ^{n} ∥}_{T}^{2} \leq C (τ^{2} + h^{4}) . \end{matrix}

(58)

So, we have

∥ θ^{n} ∥_{T}^{2} \leq C (τ^{2} + h^{4}),

∥ ξ^{n} ∥_{T}^{2} \leq C (τ^{2} + h^{4}) .

By (33) and triangle inequality,

∥ U^{n} - u^{n} ∥_{T}^{2} \leq C (τ^{2} + h^{4}), {∥ P^{n} - p^{n} ∥}_{M}^{2} \leq C (τ^{2} + h^{4}),

∥ V^{n} - v^{n} ∥_{T}^{2} \leq C (τ^{2} + h^{4}), {∥ W^{n} - w^{n} ∥}_{M}^{2} \leq C (τ^{2} + h^{4}),

hold for w sufficiently smooth.

We can draw the following conclusion.

Theorem 2.

Let w be sufficiently smooth and satisfy (12). If U, P, V and W satisfy

(14) - (17)

, for all

n, n = 1, 2, \dots, J

, then

∥ U^{n} - u^{n} ∥_{T} = O (τ + h^{2}), {∥ P^{n} - p^{n} ∥}_{M} = O (τ + h^{2}),

∥ V^{n} - v^{n} ∥_{T} = O (τ + h^{2}), {∥ W^{n} - w^{n} ∥}_{M} = O (τ + h^{2}) .

5. Time-fractional fourth-order parabolic equation

In this section, we consider the block-centered finite difference method for time-fractional fourth-order parabolic equation when

0 < α < 1

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

\begin{matrix} \{\begin{matrix} _{0}^{c} D_{t}^{α} w (x, t) + μ \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, \hat{T}], \\ w (x, 0) = φ (x), \\ \frac{\partial w (x, t)}{\partial x} {|_{x = x_{L}} = 0, \frac{\partial w (x, t)}{\partial x} |}_{x = x_{R}} = 0, \\ \frac{\partial}{\partial x} (a (x) \frac{\partial^{2} w (x, t)}{\partial x^{2}}) {|_{x = x_{L}} = 0, \frac{\partial}{\partial x} (a (x) \frac{\partial^{2} w (x, t)}{\partial x^{2}}) |}_{x = x_{R}} = 0 . \end{matrix} \end{matrix}

(59)

Where

μ

is a constant,

a (x)

f (x, t)

and

φ (x)

are known smooth functions and assume

0 < a_{0} \leq a (x) \leq a_{1}

We consider the case

0 < α < 1

_{0}^{c} D_{t}^{α} w

in (59) is defined as the following Caputo fractional derivative of

α

, given by

\begin{matrix} _{0}^{c} D_{t}^{α} w (x, t) = \frac{1}{Γ (1 - α)} \int_{0}^{t} {(t - τ)}^{- α} \frac{\partial w (x, τ)}{\partial τ} d τ . \end{matrix}

5.1. block-centered finite difference scheme

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

Lemma 3.

([36]). Suppose

f (t) \in C^{2} [0, \hat{T}]

0 < α < 1

\begin{matrix} _{0}^{c} D_{t}^{α} f (t_{n}) = D_{τ}^{α} f (t_{n}) + R (f (t_{n})) = \frac{τ^{- α}}{Γ (2 - α)} [a_{0}^{(α)} f (t_{n}) - \sum_{k = 1}^{n - 1} (a_{n - k - 1}^{(α)} - a_{n - k}^{(α)}) f (t_{k}) - a_{n - 1}^{(α)} f (t_{0})] \\ + R (f (t_{n})), \end{matrix}

where

\begin{matrix} a_{k}^{(α)} = {(k + 1)}^{1 - α} - k^{1 - α}, \\ | R (f (t_{n})) | \leq \frac{1}{2 Γ (1 - α)} [\frac{1}{4} + \frac{α}{(1 - α) (2 - α)}] max_{t_{0} \leq t \leq t_{n}} | f^{''} (t) | τ^{2 - α} . \end{matrix}

Lemma 4.

([36]). Given

0 < α < 1

, we have

\underset{l \to + \infty}{l i m} a_{l} = 0

\begin{matrix} 1 = a_{0}^{(α)} > a_{1}^{(α)} > a_{2}^{(α)} > \dots > a_{l}^{(α)} > 0 . \end{matrix}

The block-centered finite difference method (59) defines

{U_{i + 1 / 2}^{n}}_{n = 1}^{M}

{V_{i + 1 / 2}^{n}}_{n = 1}^{M}

{W_{i}^{n}}_{n = 1}^{M}

and

{P_{i}^{n}}_{n = 1}^{M} (i = 1, 2, \dots, N)

, satisfying

\begin{matrix} {[D_{τ}^{α} W]}_{i}^{n} + μ {[D_{x} U]}_{i}^{n} = f_{i}^{n}, & U_{1 / 2}^{n} = 0, U_{N + 1 / 2}^{n} = 0, \end{matrix}

(60)

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

(61)

\begin{matrix} {[D_{x} V]}_{i}^{n} = \frac{P_{i}^{n}}{a_{i}}, & V_{1 / 2}^{n} = 0, V_{N + 1 / 2}^{n} = 0, \end{matrix}

(62)

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

(63)

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{matrix} {(D_{τ}^{α} W^{n}, 1)}_{Ω_{i}} + μ {({(U^{n})}^{'}, 1)}_{Ω_{i}} = {(f^{n}, 1)}_{M, Ω_{i}}, \end{matrix}

(64)

\begin{matrix} {(U^{n}, χ)}_{T} - (P^{n}, χ^{'}) = 0, χ \in \tilde{S}, \end{matrix}

(65)

\begin{matrix} {({(V^{n})}^{'}, 1)}_{Ω_{i}} = {(\frac{P^{n}}{a}, 1)}_{M, Ω_{i}}, \end{matrix}

(66)

\begin{matrix} {(V^{n}, χ)}_{T} - (W^{n}, χ^{'}) = 0, χ \in \tilde{S} . \end{matrix}

(67)

5.2. Stability Analysis

In this subsection, we prove the stability of the scheme

(60) - (63)

when

0 < α < 1

Theorem 3.

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

∥ W^{n} ∥_{M} \leq C_{1} ∥ W^{0} ∥_{M} + C_{2} max_{1 \leq k \leq J} {∥ f^{k} ∥}_{M} .

Proof.

For

n = 1

, multiplying (60) by

h_{i} W_{i}^{1}

and summing on i from 1 to N, we have

{(D_{τ}^{α} W^{1}, W^{1})}_{M} + μ {(W_{x x}^{1}, W_{x x}^{1})}_{M} = {(f^{1}, W^{1})}_{M}

. So,

{(a_{0}^{(α)} W^{1}, W^{1})}_{M} \leq {(a_{0}^{(α)} W^{0}, W^{1})}_{M} + ρ {(f^{1}, W^{1})}_{M},

where

ρ = τ^{α} Γ (2 - α)

. Using the Cauchy-Schwarz inequality and Young inequality, we can get

(1 - \frac{ρ}{2}) ∥ W^{1} ∥_{M}^{2} \leq {∥ W^{0} ∥}_{M}^{2} + 2 ρ {∥ f^{1} ∥}_{M}^{2} .

So, we have

∥ W^{1} ∥_{M}^{2} \leq \frac{2}{2 - ρ} ∥ W^{0} ∥_{M}^{2} + \frac{4 ρ}{2 - ρ} {∥ f^{1} ∥}_{M}^{2} .

Thus, we can get

∥ W^{1} ∥_{M} \leq C_{1} ∥ W^{0} ∥_{M} + C_{2} {∥ f^{1} ∥}_{M} .

For

n \geq 2

, we suppose that the stability conclusion of the difference scheme is valid when

k \leq n - 1

Then multiplying (60) by

h_{i} W_{i}^{n}

and summing on i from 1 to N, we can get

{(D_{τ}^{α} W^{n}, W^{n})}_{M} \leq {(f^{n}, W^{n})}_{M}

, so

{(a_{0}^{(α)} W^{n}, W^{n})}_{M} \leq {(\sum_{k = 1}^{n - 1} (a_{n - k - 1}^{(α)} - a_{n - k}^{(α)}) W^{k}, W^{n})}_{M} + {(a_{n - 1}^{(α)} W^{0}, W^{n})}_{M} + ρ {(f^{n}, W^{n})}_{M} .

By the Cauchy-Schwarz inequality and Young inequality, we obtain

∥ W^{n} ∥_{M}^{2} \leq \frac{2}{2 - ρ} \sum_{k = 1}^{n - 1} (a_{n - k - 1}^{(α)} - a_{n - k}^{(α)}) ∥ W^{k} ∥_{M}^{2} + \frac{2}{2 - ρ} a_{n - 1}^{(α)} ∥ W^{0} ∥_{M}^{2} + \frac{4 ρ}{2 - ρ} {∥ f^{n} ∥}_{M}^{2} .

By mathematical induction and using the relation of coefficient

a_{k}

, we have

\begin{matrix} ∥ W^{n} ∥_{M}^{2} & \leq \frac{2}{2 - ρ} \sum_{k = 1}^{n - 1} (a_{n - k - 1}^{(α)} - a_{n - k}^{(α)}) (C_{1}^{'} ∥ W^{0} ∥_{M}^{2} + C_{2}^{'} max_{1 \leq k \leq J} ∥ f^{k} ∥_{M}^{2}) + \frac{2}{2 - ρ} a_{n - 1}^{(α)} {∥ W^{0} ∥}_{M}^{2} \\ + \frac{4 ρ}{2 - ρ} {∥ f^{n} ∥}_{M}^{2} \\ \leq \frac{2}{2 - ρ} [(C_{1}^{'} (a_{0}^{(α)} - a_{n - 1}^{(α)}) + a_{n - 1}^{(α)}) ∥ W^{0} ∥_{M}^{2} + (C_{2}^{'} (a_{0}^{(α)} - a_{n - 1}^{(α)}) + 2 ρ) max_{1 \leq k \leq J} ∥ f^{k} ∥_{M}^{2}] . \end{matrix}

There are constants

C_{3}

and

C_{4}

that make

∥ W^{n} ∥_{M}^{2} \leq C_{3} ∥ W^{0} ∥_{M}^{2} + C_{4} max_{1 \leq k \leq J} {∥ f^{k} ∥}_{M}^{2} .

So, we can get

∥ W^{n} ∥_{M} \leq C_{1} ∥ W^{0} ∥_{M} + C_{2} max_{1 \leq k \leq J} {∥ f^{k} ∥}_{M} .

We complete the proof. □

5.3. Error Analysis

Then, the error analysis of the block-centered difference scheme for time-fractional parabolic equation is carried out.

Error estimates for this finite difference scheme of

(60) - (63)

are derived using a technique of mixed finite element methods for parabolic partial differential equations.

\begin{matrix} μ \frac{\partial^{2}}{\partial x^{2}} (a (x) \frac{\partial^{2} w^{n}}{\partial x^{2}}) = ϕ = f^{n} -_{0}^{c} D_{t}^{α} w^{n}, & x \in (x_{L}, x_{R}) \end{matrix}

(68)

\begin{matrix} \frac{\partial w^{n}}{\partial x} {|_{x = x_{L}} = 0, \frac{\partial w^{n}}{\partial x} |}_{x = x_{R}} = 0, \end{matrix}

(69)

\begin{matrix} \frac{\partial}{\partial x} (a (x) \frac{\partial^{2} w^{n}}{\partial x^{2}}) {|_{x = x_{L}} = 0, \frac{\partial}{\partial x} (a (x) \frac{\partial^{2} w^{n}}{\partial x^{2}}) |}_{x = x_{R}} = 0 . \end{matrix}

(70)

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{matrix} μ {[D_{x} E_{h} u]}_{i}^{n} = f_{i}^{n} - {(_{0}^{c} D_{t}^{α} w)}_{i}^{n}, & E_{h} u_{1 / 2}^{n} = 0, E_{h} u_{N + 1 / 2}^{n} = 0, \end{matrix}

(71)

\begin{matrix} E_{h} u_{i + 1 / 2}^{n} = {[- d_{x} E_{h} p]}_{i + 1 / 2}^{n}, \end{matrix}

(72)

\begin{matrix} {[D_{x} E_{h} v]}_{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, \end{matrix}

(73)

\begin{matrix} E_{h} v_{i + 1 / 2}^{n} = {[- d_{x} E_{h} w]}_{i + 1 / 2}^{n}, \end{matrix}

(74)

where

E_{h} u^{0} = {(a (x) φ^{''} (x))}^{'}

E_{h} p^{0} = - a (x) φ^{''} (x)

E_{h} v^{0} = - φ^{'} (x)

and

E_{h} w^{0} = φ (x)

Equations

(71 - (74)

can be written as a mixed finite element method with approximate integration

\begin{matrix} μ {({(E_{h} u^{n})}^{'}, 1)}_{Ω_{i}} = {(f^{n}, 1)}_{M, Ω_{i}} - {(_{0}^{c} D_{t}^{α} w^{n}, 1)}_{M, Ω_{i}} \end{matrix}

(75)

\begin{matrix} {(E_{h} u^{n}, χ)}_{T} - (E_{h} p^{n}, χ^{'}) = 0, χ \in \tilde{S}, \end{matrix}

(76)

\begin{matrix} {({(E_{h} v^{n})}^{'}, 1)}_{Ω_{i}} = {(\frac{E_{h} p^{n}}{a}, 1)}_{M, Ω_{i}}, \end{matrix}

(77)

\begin{matrix} {(E_{h} v^{n}, χ)}_{T} - (E_{h} w^{n}, χ^{'}) = 0, χ \in \tilde{S} . \end{matrix}

(78)

Set

ξ_{i + 1 / 2}^{n} = U_{i + 1 / 2}^{n} - E_{h} u_{i + 1 / 2}^{n}

η_{i}^{n} = P_{i}^{n} - E_{h} p_{i}^{n}

θ_{i + 1 / 2}^{n} = V_{i + 1 / 2}^{n} - E_{h} v_{i + 1 / 2}^{n}

ζ_{i}^{n} = W_{i}^{n} - E_{h} w_{i}^{n}

, and

σ_{i}^{n} = E_{h} w_{i}^{n} - w_{i}^{n} .

Subtracting (71) from (60), we get

\begin{matrix} D_{τ}^{α} ζ_{i}^{n} + μ {[D_{x} ξ]}_{i}^{n} = R (w_{i}^{n}) - D_{τ}^{α} σ_{i}^{n}, & ξ_{1 / 2}^{n} = 0, ξ_{N + 1 / 2}^{n} = 0 . \end{matrix}

(79)

Subtracting (72), (73) and (74) from (61), (62) and (63), correspondingly, we obtain

\begin{matrix} ξ_{i + 1 / 2}^{n} = {[- d_{x} η]}_{i + 1 / 2}^{n}, \end{matrix}

(80)

\begin{matrix} η_{i}^{n} = a_{i} {[D_{x} θ]}_{i}^{n}, & θ_{1 / 2}^{n} = 0, θ_{N + 1 / 2}^{n} = 0, \end{matrix}

(81)

\begin{matrix} θ_{i + 1 / 2}^{n} = {[- d_{x} ζ]}_{i + 1 / 2}^{n} . \end{matrix}

(82)

Multiplying (79) by

h_{i} ζ_{i}^{n}

and summing on i from 1 to N, we deduce that

\begin{matrix} {(D_{τ}^{α} ζ^{n}, ζ^{n})}_{M} + μ {(D_{x} ξ^{n}, ζ^{n})}_{M} = {(R (w^{n}), ζ^{n})}_{M} - {(D_{τ}^{α} σ^{n}, ζ^{n})}_{M} . \end{matrix}

(83)

By (80), (81), (82) and Lemma 2, we have

\begin{matrix} {(D_{x} ξ^{n}, ζ^{n})}_{M} & = \frac{1}{a} {(η^{n}, η^{n})}_{M} . \end{matrix}

(84)

Now

\begin{matrix} {(D_{τ}^{α} ζ^{n}, ζ^{n})}_{M} = \frac{τ^{- α}}{Γ (2 - α)} {[a_{0}^{(α)} (ζ^{n}, ζ^{n}) - \sum_{k = 1}^{n - 1} (a_{n - k - 1}^{(α)} - a_{n - k}^{(α)}) (ζ^{k}, ζ^{n}) - a_{n - 1}^{(α)} (ζ^{0}, ζ^{n})]}_{M} . \end{matrix}

(85)

Let

ρ = τ^{α} Γ (2 - α)

, (83) can be written as

\begin{matrix} ∥ ζ^{n} ∥_{M}^{2} + \frac{μ ρ}{a_{1}} {∥ η^{n} ∥}_{M}^{2} \leq \sum_{k = 1}^{n - 1} (a_{n - k - 1}^{(α)} - a_{n - k}^{(α)}) {(ζ^{k}, ζ^{n})}_{M} + a_{n - 1}^{(α)} {(ζ^{0}, ζ^{n})}_{M} \\ + ρ | {(R (w^{n}), ζ^{n})}_{M} | + ρ | {(D_{τ}^{α} σ^{n}, ζ^{n})}_{M} | . \end{matrix}

(86)

By Cauchy-Schwarz inequality and Young inequality, we have

\begin{matrix} (1 - \frac{ρ}{2}) ∥ ζ^{n} ∥_{M}^{2} + \frac{2}{a_{1}} μ ρ ∥ η^{n} ∥_{M}^{2} \leq \sum_{k = 1}^{n - 1} (a_{n - k - 1}^{(α)} - a_{n - k}^{(α)}) ∥ ζ^{k} ∥_{M}^{2} + a_{n - 1}^{(α)} {∥ ζ^{0} ∥}_{M}^{2} \\ + 4 ρ ∥ D_{τ}^{α} σ^{n} ∥_{M}^{2} + 4 ρ {∥ R (w^{n}) ∥}_{M}^{2} . \end{matrix}

(87)

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

\begin{matrix} ∥ D_{τ}^{α} σ^{n} ∥_{M} & = \frac{1}{Γ (1 - α)} {∥ \sum_{k = 1}^{n} \frac{σ^{k} - σ^{k - 1}}{τ} \int_{t_{k - 1}}^{t_{k}} {(t_{n} - t)}^{- α} d t ∥}_{M} \\ = \frac{1}{Γ (1 - α)} {∥ \sum_{k = 1}^{n} \frac{1}{τ} \int_{t_{k - 1}}^{t_{k}} σ_{t} d t \int_{t_{k - 1}}^{t_{k}} {(t_{n} - t)}^{- α} d t ∥}_{M} \\ \leq \frac{1}{Γ (1 - α)} max_{0 \leq t \leq t_{n}} {∥ σ_{t} ∥}_{M} \sum_{k = 1}^{n} \int_{t_{k - 1}}^{t_{k}} {(t_{n} - t)}^{- α} d t \\ \leq \frac{T^{1 - α}}{Γ (2 - α)} max_{0 \leq t \leq t_{n}} {∥ σ_{t} ∥}_{M} \\ \leq C h^{2} . \end{matrix}

(88)

Notice that

ζ^{0} = 0

, by Theorem 3, using inductive hypothesis, we deduce that

\begin{matrix} (1 - \frac{ρ}{2}) ∥ ζ^{n} ∥_{M}^{2} + \frac{2}{a_{1}} μ ρ {∥ η^{n} ∥}_{M}^{2} \leq C (τ^{4 - 2 α} + h^{4}) . \end{matrix}

(89)

Thus, we have

\begin{matrix} ∥ ζ^{n} ∥_{M}^{2} \leq C (τ^{4 - 2 α} + h^{4}), \end{matrix}

(90)

\begin{matrix} ∥ η^{n} ∥_{M}^{2} \leq C (τ^{4 - 2 α} + h^{4}) . \end{matrix}

(91)

Multiplying (79) by

h_{i} η_{i}^{n}

, and summing on i from 1 to N, we can get

\begin{matrix} {(D_{τ}^{α} ζ^{n}, η^{n})}_{M} + μ {(D_{x} ξ^{n}, η^{n})}_{M} = {(R (w^{n}), η^{n})}_{M} - {(D_{τ}^{α} σ^{n}, η^{n})}_{M} . \end{matrix}

(92)

Now

\begin{matrix} {(D_{τ}^{α} ζ^{n}, η^{n})}_{M} = \frac{τ^{- α}}{Γ (2 - α)} {[a_{0}^{(α)} (ζ^{n}, η^{n}) - \sum_{k = 1}^{n - 1} (a_{n - k - 1}^{(α)} - a_{n - k}^{(α)}) (ζ^{k}, η^{n}) - a_{n - 1}^{(α)} (ζ^{0}, η^{n})]}_{M} . \end{matrix}

(93)

Using (80), (81), (82) and Lemma 2, we derive

\begin{matrix} {(ζ^{n}, η^{n})}_{M} = a {(θ^{n}, θ^{n})}_{T}, \end{matrix}

(94)

\begin{matrix} {(D_{x} ξ^{n}, η^{n})}_{M} = {(ξ^{n}, ξ^{n})}_{T} . \end{matrix}

(95)

Substituting (93), (94) and (95) into (92), (92) can be deformable to

\begin{matrix} a_{0} ∥ θ^{n} ∥_{T}^{2} + μ ρ {∥ ξ^{n} ∥}_{T}^{2} \leq a_{1} \sum_{k = 1}^{n - 1} (a_{n - k - 1}^{(α)} - a_{n - k}^{(α)}) {(θ^{k}, θ^{n})}_{T} + a_{1} a_{n - 1}^{(α)} {(θ^{0}, θ^{n})}_{T} \\ + ρ | {(R (w^{n}), η^{n})}_{M} | + ρ | {(D_{τ}^{α} σ^{n}, η^{n})}_{M} | . \end{matrix}

(96)

By Cauchy-Schwarz inequality, we derive that

\begin{matrix} (2 a_{0} - a_{1}) ∥ θ^{n} ∥_{T}^{2} + 2 μ ρ ∥ ξ^{n} ∥_{T}^{2} \leq a_{1} \sum_{k = 1}^{n - 1} (a_{n - k - 1}^{(α)} - a_{n - k}^{(α)}) ∥ θ^{k} ∥_{T}^{2} + a_{1} a_{n - 1}^{(α)} ∥ θ^{0} ∥_{T}^{2} + 2 ρ {∥ η^{n} ∥}_{M}^{2} + \\ ρ ∥ D_{τ}^{α} σ^{n} ∥_{M}^{2} + ρ {∥ R (w^{n}) ∥}_{M}^{2} . \end{matrix}

(97)

Similarly, when

2 a_{0} > a_{1}

, using the above mathematical induction and substituting (88)and (91) into (97), we can get

\begin{matrix} (2 a_{0} - a_{1}) ∥ θ^{n} ∥_{T}^{2} + 2 μ ρ {∥ ξ^{n} ∥}_{T}^{2} \leq C (τ^{4 - 2 α} + h^{4}) . \end{matrix}

(98)

So, we can get

∥ θ^{n} ∥_{T}^{2} \leq C (τ^{4 - 2 α} + h^{4}),

∥ ξ^{n} ∥_{T}^{2} \leq C (τ^{4 - 2 α} + h^{4}) .

By (33) and triangle inequality,

∥ U^{n} - u^{n} ∥_{T}^{2} \leq C (τ^{4 - 2 α} + h^{4}), {∥ P^{n} - p^{n} ∥}_{M}^{2} \leq C (τ^{4 - 2 α} + h^{4}),

∥ V^{n} - v^{n} ∥_{T}^{2} \leq C (τ^{4 - 2 α} + h^{4}), {∥ W^{n} - w^{n} ∥}_{M}^{2} \leq C (τ^{4 - 2 α} + h^{4}),

hold for w sufficiently smooth.

We can get the following conclusion.

Theorem 4.

Let w be sufficiently smooth and satisfy (59). If U, P, V and W satisfy

(60) - (63)

, for all

n, n = 1, 2, \dots, J

, then

∥ U^{n} - u^{n} ∥_{T} = O (τ^{2 - α} + h^{2}), {∥ P^{n} - p^{n} ∥}_{M} = O (τ^{2 - α} + h^{2}),

∥ V^{n} - v^{n} ∥_{T} = O (τ^{2 - α} + h^{2}), {∥ W^{n} - w^{n} ∥}_{M} = O (τ^{2 - α} + h^{2}) .

6. Numerical Experiments

In this section, we give 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{matrix} \{\begin{matrix} w^{''''} = f (x), x \in (0, 1), \\ w^{'} {|_{x = 0} = 0, w^{'} |}_{x = 1} = 0, \\ w^{'''} {|_{x = 0} = 0, w^{'''} |}_{x = 1} = 0 . \end{matrix} \end{matrix}

(99)

We give the exact solution to the problem

w (x) = cos (π x)

and a source term

f (x) = π^{4} cos (π x) .

Take the spatial step is

h = 1 / N .

The space errors and convergence orders of W, V, P and U are shown in Table 1, Table 2, Table 3 and Table 4, 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, while the third and fifth columns show their corresponding spatial convergence order respectively. It can be seen from the table that the order of spatial convergence reaches the second order under the maximum norm and discrete L2 error, which verifies the convergence of the block-centered finite difference scheme.

Figure 1 shows the image of the numerical solution and the exact solution when L=3. It can be seen from the figure that the numerical solution can better fit the exact solution.

Example 2.

In equation (12), let μ,

a (x)

and

\hat{T}

all be equal to 1. Consider the following initial boundary value problem

\begin{matrix} \{\begin{matrix} w_{t} (x, t) + \frac{\partial^{4} w (x, t)}{\partial x^{4}} = f (x, t), x \in (0, 1), t \in (0, 1], \\ w (x, 0) = 0, x \in (0, 1), \\ \frac{\partial w (x, t)}{\partial x} {|_{x = 0} = 0, \frac{\partial w (x, t)}{\partial x} |}_{x = 1} = 0, t \in (0, 1], \\ \frac{\partial^{3} w (x, t)}{\partial x^{3}} {|_{x = 0} = 0, \frac{\partial^{3} w (x, t)}{\partial x^{3}} |}_{x = 1} = 0, t \in (0, 1] . \end{matrix} \end{matrix}

(100)

The exact solution is

w (x, t) = t^{2} cos (π x)

, the initial value is

w (x, 0) = 0

, and the interval is (0,1). The spatial step is

h = 1 / N

and the time step is

τ = 1 / J

. The time and space errors and convergence orders of W, V, P and U are shown in Table 5, Table 6, Table 7, Table 8, Table 9, Table 10, Table 11 and Table 12 , respectively. The second and fourth columns of Table 5, Table 6, Table 7 and Table 8 show the maximum norms error and the discrete L2-norms error of W, V, P and U as the mesh size h is reduced with fixing

τ = 1 / 2^{16}

. The third and fifth columns of Table 5, Table 6, Table 7 and Table 8 show the spatial convergence order of their corresponding the maximum norms error and the discrete L2-norms error. The second and fourth columns of Table 9, Table 10, Table 11 and Table 12 show the maximum norms error and the discrete L2-norms error of W, V, P and U as the mesh size τ is reduced with fixing

h = 1 / 2^{10}

. The third and fifth columns of Table 9, Table 10, Table 11 and Table 12 show the time convergence order of their corresponding the maximum norms error and the discrete L2-norms error. It can be seen from the tables the space convergence order has reached the second order, and the time convergence order has reached the first order, the validity of the block-centered finite difference scheme is verified.

Figure 2 and Figure 3 show the images of the numerical solution W and the analytic solution w when the mesh is divided into

N = 2^{4}

J = N^{2} = 2^{8}

. As can be seen from the figures, the fitting effect between them is good.

Example 3.

In equation (59), let

\hat{T}

and

a (x)

be equal to 1 and μ equal to 10. Consider the following initial boundary value problem

\begin{matrix} \{\begin{matrix} _{0}^{c} D_{t}^{α} w (x, t) + 10 \frac{\partial^{4} w (x, t)}{\partial x^{4}} = f (x, t), x \in (0, 1), t \in (0, 1], \\ w (x, 0) = 0, \\ \frac{\partial w (x, t)}{\partial x} {|_{x = 0} = 0, \frac{\partial w (x, t)}{\partial x} |}_{x = 1} = 0, \\ \frac{\partial^{3} w (x, t)}{\partial x^{3}} {|_{x = 0} = 0, \frac{\partial^{3} w (x, t)}{\partial x^{3}} |}_{x = 1} = 0 . \end{matrix} \end{matrix}

(101)

The exact solution is

w (x, t) = t^{2} cos (π x)

. The spatial step is

h = 1 / N

and the time step is

τ = 1 / J

. When

α = 0.4, 0.6, 0.8

, the time and space errors and convergence orders of W, V, P and U are shown in Table 13, Table 14, Table 15, Table 16, Table 17, Table 18, Table 19 and Table 20. Table 13, Table 14, Table 15 and Table 16 show the maximum norms error, the discrete L2-norms error and convergence order of W, V, P and U as the mesh size h is reduced with fixing

τ = 1 / 2^{10}

. Table 17, Table 18, Table 19 and Table 20 show the maximum norms error, the discrete L2-norms error and convergence order of W, V, P and U as the mesh size τ is reduced with fixing

h = 1 / 2^{10}

. It can be seen from the tables that for the maximum norm and the discrete L2-norm, the space convergence order has reached the second order, and the time convergence order has reached the

2 - α

order, the validity of the block-centered finite difference scheme is verified.

Figure 4 and Figure 5 show the images of the numerical solution W and the analytic solution w when the mesh is divided into

N = J^{2} = 2^{8}

. It can be seen from the two figures that the numerical solution W fits well with the analytic solution w.

7. Conclusions

This paper discusses the block-centered finite difference method for solving fourth-order parabolic equations of integer and 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 proven. The effectiveness of the block-centered finite difference scheme is verified by numerical examples. In future work, we will further study to solve other types of fourth-order partial differential equations using the block-centered finite difference method and numerical examples on non-uniform grids.

Author Contributions

Conceptualization, Z.Y. and A.Z.; writing-original draft, T.Z.; writing-review and editing, T.Z. All authors have read and agreed to the published version of the manuscript.

Funding

The work is supported by the National Natural Science Foundation of China (contract grant number: 12171287) and the Natural Science Foundation of Shandong Province (contract grant number: ZR2021MA063).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to express their sincere thanks to the referees for the valuable comments and suggestions that helped to improve the original paper.

Conflicts of Interest

The authors declare no conflict of interest regarding the publication of this paper.

References

Russell, T.F.; Wheeler, M.F. Finite Element and Finite Difference Methods for Continuous Flows in Porous Media. In The mathematics of reservoir simulation; Society for Industrial Mathematics: Philadelphia, USA, 1987. [Google Scholar]
Weiser, A.; Wheeler, M.F. On Convergence of Block-Centered Finite Differences for Elliptic Problems. SIAM J. Numer. Anal. 1988, 25, 351–375. [Google Scholar] [CrossRef]
Rui, H.; Pan, H. A Block-Centered Finite Difference Method for the Darcy-Forchheimer Model. SIAM J. Numer. Anal. 2012, 50, 2612–2631. [Google Scholar] [CrossRef]
Rui, H.; Pan, H. Block-Centered Finite Difference Methods for Parabolic Equation with Time-Dependent Coefficient. Japan J. Indust. Appl. Math. 2013, 30, 681–699. [Google Scholar] [CrossRef]
Rui, H.; Zhao, D.; Pan, H. A Block-Centered Finite Difference Method for Darcy-Forchheimer Model with Variable Forchheimer Number. Numer. Methods Partial. Differ. Equ. 2015, 31, 1603–1622. [Google Scholar] [CrossRef]
Zhai, S.; Qian, L.; Gui, D. A Block-Centered Characteristic Finite Difference Method for Convection-Dominated Diffusion Equation. Int. Commun. Heat Mass Transf. 2015, 61, 1–7. [Google Scholar] [CrossRef]
Liu, Y. Block-Centered Finite Difference Method for the Numerical Simulation of a Semiconductor Device. Numer. Math. J. Chinese Univ 1995, 6, 108–118. [Google Scholar]
Yuan, Y.; Liu, Y.; Li, C. Analysis on Block-Centered Finite Differences of Numerical Simulation of Semiconductor Device Detector. Appl. Math. Comput. 2016, 279, 1–15. [Google Scholar] [CrossRef]
Li, X.; Rui, H. Characteristic Block-Centred Finite Difference Methods for Nonlinear Convection-Dominated Diffusion Equation. Int. J Comput. Math. 2017, 94, 386–404. [Google Scholar] [CrossRef]
Zhai, S.; Weng, Z.; Feng, X. An Adaptive Local Grid Refinement Method for 2D Diffusion Equation with Variable Coefficients Based on Block-Centered Finite Difference. Appl. Math. Comput. 2015, 268, 284–294. [Google Scholar] [CrossRef]
Li, X.; Rui, H. Block-centered Finite Difference Method for Simulating Compressible Wormhole Propagation. J. Sci. Comput. 2018, 74, 1115–1145. [Google Scholar] [CrossRef]
Ren, Z.; Zhang, X.; Yin, Z. Crank-Nicolson Block-centered Difference Scheme for Parabolic Problems. J. Henan Norm Univ: Nat SCi Ed. 2011, 39, 15–19. [Google Scholar]
Xu, W. Cell-Centered Finite Difference Method for Parabolic Equation. Appl. Math. Comput. 2014, 235, 66–79. [Google Scholar] [CrossRef]
Liu, Y. Error Estimates of Block-Centered Finite Difference Methods For a Class of Parabolic Integro-Differential Equations. J. Shandong Univ. 1994, 29, 257–265. [Google Scholar]
Dawson, C.N.; Dupont, T.F. Explicit/Implicit, Conservative Domain Decomposition Procedures for Parabolic Problems Based on Block-Centered Finite Differences. SIAM J. Numer. Anal. 1994, 31, 1045–1061. [Google Scholar] [CrossRef]
Li, X.; Rui, H. A Block-Centered Finite Difference Method for the Distributed-Order Time-Fractional Diffusion-Wave Equation. Appl. Numer. Math. 2018, 131, 123–139. [Google Scholar] [CrossRef]
Li, X.; Rui, H. A Two-Grid Block-Centered Finite Difference Method for the Nonlinear Time-Fractional Parabolic Equation. J. Sci. Comput. 2017, 72, 863–891. [Google Scholar] [CrossRef]
Liu, Z.; Li, X. A Parallel CGS Block-Centered Finite Difference Method for a Nonlinear Time-Fractional Parabolic Equation. Comput. Methods Appl. Mech. Eng. 2016, 308, 330–348. [Google Scholar] [CrossRef]
Li, X.; Chen,Y. ; Chen, C. An Improved Two-Grid Technique for the Nonlinear Time-Fractional Parabolic Equation Based on the Block-Centered Finite Difference Method. J. Comput. Math. 2022, 40, 455–473. [Google Scholar] [CrossRef]
Zhai, S.; Feng, X. A Block-Centered Finite-Difference Method for the Time-Fractional Diffusion Equation on Nonuniform Grids. Numer. Heat Transf. Part B Fundam. 2016, 69, 217–233. [Google Scholar] [CrossRef]
Guo, Q.; Rui, H. Block-Centered Local Refinement Methods for the Time-Fractional Equations. Chaos Solitons Fractals 2021, 152, 111–314. [Google Scholar] [CrossRef]
Shi, Y.; Xie, S.; Liang, D.; Fu, K. High Order Compact Block-Centered Finite Difference Schemes for Elliptic and Parabolic Problems. J. Sci. Comput. 2021, 87, 1–26. [Google Scholar] [CrossRef]
Dee, G.T.; Saarloos, Van. W. Bistable Systems with Propagating Fronts Leading to Pattern Formation. Phys. Rev. Lett. 1988, 60, 2641–2644. [Google Scholar] [CrossRef] [PubMed]
Tayler, A.B. Mathematical Model in Applied Mechanics; Oxford Clarendon Press: New York, USA, 1986. [Google Scholar]
Jüngel, A. A Positivity-Preserving Numerical Scheme for a Nonlinear Fourth Order Parabolic System. SIAM J. Numer. Anal. 2001, 39, 385–406. [Google Scholar] [CrossRef]
Kaur, D.; Mohanty, R.K. Highly accurate compact difference scheme for fourth order parabolic equation with Dirichlet and Neumann boundary conditions: Application to good Boussinesq equation. Appl. Math. Comput. 2020, 378, 125–202. [Google Scholar] [CrossRef]
Liang, B.; Peng, X.; Shen, H. Study of Solutions to a Fourth Order Parabolic Equation. Math. Modell. Anal. 2016, 21, 1–15. [Google Scholar] [CrossRef]
Acharya, S.; Porwal, K. Primal hybrid finite element method for fourth order parabolic problems. Appl. Numer. Math. 2020, 152, 12–28. [Google Scholar] [CrossRef]
Keita, S.; Beljadid, A.; Bourgault, Y. Efficient second-order semi-implicit finite element method for fourth-order nonlinear diffusion equations. Comput. Phys. Commun. 2021, 258, 107588. [Google Scholar] [CrossRef]
Aziz, S.; Malik, S. Identification of an unknown source term for a time fractional fourth-order parabolic equation. Electron. J. Differ. Equ. 2016, 293, 1–28. [Google Scholar]
Li, D.; Liao, H.; Sun, W. Analysis of L1-Galerkin FEMs for time-fractional nonlinear parabolic problems. Commun. Comput. Phys. 2018, 24, 86–103. [Google Scholar] [CrossRef]
Zhu, A.; Wang, Y.; Xu, Q. A weak Galerkin finite element approximation of two-dimensional sub-diffusion equation with time-fractional derivative. AIMS Math. 2020, 5, 4297–4310. [Google Scholar] [CrossRef]
Liu, Y.; Du, Y.; Li, H. Finite difference/finite element method for a nonlinear time-fractional fourth-order reaction-diffusion problem. Comput. Math. Appl. 2015, 70, 573–591. [Google Scholar] [CrossRef]
Hu, X.; Zhang, L. A new implicit compact difference scheme for the fourth-order fractional diffusion-wave system. Int. J Comput. Math. 2014, 91, 2215–2231. [Google Scholar] [CrossRef]
Ji, C.; Sun, Z.; Hao, Z. Numerical algorithms with high spatial accuracy for the fourth-order fractional sub-diffusion equations with the first Dirichlet boundary conditions. J. Sci. Comput. 2016, 66, 1148–1174. [Google Scholar] [CrossRef]
Sun, Z.; Gao, G. Finite Difference Method for Fractional Differential Equations; Science Press: Beijing, China, 2015. [Google Scholar]

Figure 1. Comparison of the numerical solution and the exact solution when

L = 3, h = \frac{1}{32}

Figure 1. Comparison of the numerical solution and the exact solution when

L = 3, h = \frac{1}{32}

Figure 2. The numerical solution W when

h = \frac{1}{16}

Figure 2. The numerical solution W when

h = \frac{1}{16}

Figure 3. The exact solution w when

h = \frac{1}{16}

Figure 3. The exact solution w when

h = \frac{1}{16}

Figure 4. The numerical solution W when

α = 0.4, h = \frac{1}{16}

Figure 4. The numerical solution W when

α = 0.4, h = \frac{1}{16}

Figure 5. The exact solution w when

α = 0.4, h = \frac{1}{16}

Figure 5. The exact solution w when

α = 0.4, h = \frac{1}{16}

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

h	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
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	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
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	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
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

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

h	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
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

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

h	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
1/16	$6.2860 \times 10^{- 3}$	-	$4.4664 \times 10^{- 3}$	-
1/32	$1.5733 \times 10^{- 3}$	1.9984	$1.1138 \times 10^{- 3}$	2.0036
1/64	$3.9353 \times 10^{- 4}$	1.9992	$2.7836 \times 10^{- 4}$	2.0005
1/128	$9.8512 \times 10^{- 5}$	1.9982	$6.9662 \times 10^{- 5}$	1.9985

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

h	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
1/16	$1.4768 \times 10^{- 2}$	-	$1.0442 \times 10^{- 2}$	-
1/32	$3.6851 \times 10^{- 3}$	2.0026	$2.6058 \times 10^{- 3}$	2.0026
1/64	$9.2121 \times 10^{- 4}$	2.0001	$6.5139 \times 10^{- 4}$	2.0001
1/128	$2.3064 \times 10^{- 4}$	1.9979	$1.6309 \times 10^{- 4}$	1.9979

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

h	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
1/16	$3.0326 \times 10^{- 2}$	-	$2.1548 \times 10^{- 2}$	-
1/32	$7.6004 \times 10^{- 3}$	1.9964	$5.3808 \times 10^{- 3}$	2.0016
1/64	$1.9024 \times 10^{- 3}$	1.9983	$1.3456 \times 10^{- 3}$	1.9996
1/128	$4.7682 \times 10^{- 4}$	1.9963	$3.3719 \times 10^{- 4}$	1.9966

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

h	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
1/16	$4.5798 \times 10^{- 2}$	-	$3.2384 \times 10^{- 2}$	-
1/32	$1.1447 \times 10^{- 2}$	2.0004	$8.0940 \times 10^{- 3}$	2.0004
1/64	$2.8649 \times 10^{- 3}$	1.9984	$2.0258 \times 10^{- 3}$	1.9984
1/128	$7.1983 \times 10^{- 4}$	1.9927	$5.0900 e \times 10^{- 4}$	1.9927

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

$τ$	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
1/16	$6.4326 \times 10^{- 4}$	-	$4.5462 \times 10^{- 4}$	-
1/32	$3.2405 \times 10^{- 4}$	0.9892	$2.2779 \times 10^{- 4}$	0.9970
1/64	$1.6350 \times 10^{- 4}$	0.9869	$1.1446 \times 10^{- 4}$	0.9928
1/128	$8.2199 \times 10^{- 5}$	0.9921	$5.7822 \times 10^{- 5}$	0.9852

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

$τ$	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
1/16	$2.0187 \times 10^{- 3}$	-	$1.4275 \times 10^{- 3}$	-
1/32	$1.0108 \times 10^{- 3}$	0.9979	$7.1477 \times 10^{- 4}$	0.9979
1/64	$5.0729 \times 10^{- 4}$	0.9947	$3.5870 \times 10^{- 4}$	0.9947
1/128	$2.5567 \times 10^{- 4}$	0.9885	$1.8079 \times 10^{- 4}$	0.9885

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

$τ$	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
1/16	$6.3888 \times 10^{- 3}$	-	$4.4821 \times 10^{- 3}$	-
1/32	$3.1721 \times 10^{- 3}$	0.9988	$2.2429 \times 10^{- 3}$	0.9988
1/64	$1.5901 \times 10^{- 3}$	0.9963	$1.1242 \times 10^{- 3}$	0.9964
1/128	$7.9948 \times 10^{- 4}$	0.9920	$5.6527 \times 10^{- 4}$	0.9919

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

$τ$	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
1/16	$1.9902 \times 10^{- 2}$	-	$1.4073 \times 10^{- 2}$	-
1/32	$9.9541 \times 10^{- 3}$	0.9996	$7.0384 \times 10^{- 3}$	0.9997
1/64	$4.9829 \times 10^{- 3}$	0.9983	$3.5236 \times 10^{- 3}$	0.9982
1/128	$2.4992 \times 10^{- 3}$	0.9955	$1.7674 \times 10^{- 3}$	0.9954

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

$α$	h	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
	1/4	$8.7432 \times 10^{- 2}$	-	$4.3716 \times 10^{- 2}$	-
0.4	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
	1/4	$8.4477 \times 10^{- 2}$	-	$4.2238 \times 10^{- 2}$	-
0.6	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
	1/4	$8.2380 \times 10^{- 2}$	-	$4.1190 \times 10^{- 2}$	-
0.8	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 14. The computational errors and convergence orders of V in space.

$α$	h	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
	1/4	$2.0975 \times 10^{- 1}$	-	$1.4832 \times 10^{- 1}$	-
0.4	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
	1/4	$1.9996 \times 10^{- 1}$	-	$1.4139 \times 10^{- 1}$	-
0.6	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
	1/4	$1.9302 \times 10^{- 1}$	-	$1.3648 \times 10^{- 1}$	-
0.8	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 15. The computational errors and convergence orders of P in space.

$α$	h	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
	1/4	$3.6117 \times 10^{- 1}$	-	$2.7643 \times 10^{- 1}$	-
0.4	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
	1/4	$3.3347 \times 10^{- 1}$	-	$2.5523 \times 10^{- 1}$	-
0.6	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
	1/4	$3.1382 \times 10^{- 1}$	-	$2.4019 \times 10^{- 1}$	-
0.8	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 16. The computational errors and convergence orders of U in space.

$α$	h	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
	1/4	$4.0755 \times 10^{- 1}$	-	$2.8818 \times 10^{- 1}$	-
0.4	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
	1/4	$3.1577 \times 10^{- 1}$	-	$2.2328 \times 10^{- 1}$	-
0.6	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
	1/4	$2.5066 \times 10^{- 1}$	-	$1.7724 \times 10^{- 1}$	-
0.8	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 17. The computational errors and convergence orders of W in time.

$α$	$τ$	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
	1/8	$1.1526 \times 10^{- 3}$	-	$3.6020 \times 10^{- 5}$	-
0.4	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
	1/8	$2.7928 \times 10^{- 3}$	-	$8.7275 \times 10^{- 5}$	-
0.6	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
	1/8	$6.1370 \times 10^{- 3}$	-	$1.9178 \times 10^{- 4}$	-
0.8	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 18. The computational errors and convergence orders of V in time.

$α$	$τ$	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
	1/8	$3.6206 \times 10^{- 3}$	-	$2.5601 \times 10^{- 3}$	-
0.4	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
	1/8	$8.7728 \times 10^{- 3}$	-	$6.2033 \times 10^{- 3}$	-
0.6	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
	1/8	$1.9279 \times 10^{- 2}$	-	$1.3632 \times 10^{- 2}$	-
0.8	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 19. The computational errors and convergence orders of P in time.

$α$	$τ$	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
	1/8	$1.1372 \times 10^{- 2}$	-	$8.0408 \times 10^{- 3}$	-
0.4	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
	1/8	$2.7558 \times 10^{- 2}$	-	$1.9486 \times 10^{- 2}$	-
0.6	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
	1/8	$6.0563 \times 10^{- 2}$	-	$4.2825 \times 10^{- 2}$	-
0.8	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

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

$α$	$τ$	$l^{\infty} - E r r o r$	$O r d e r$	$l^{2} - E r r o r$	$O r d e r$
	1/8	$3.5715 \times 10^{- 2}$	-	$2.5254 \times 10^{- 2}$	-
0.4	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
	1/8	$8.6565 \times 10^{- 2}$	-	$6.1211 \times 10^{- 2}$	-
0.6	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
	1/8	$1.9026 \times 10^{- 1}$	-	$1.3453 \times 10^{- 1}$	-
0.8	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

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

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

Abstract

1. Introduction

2. Notation

3. Fourth-order ordinary differential equation

4. Time integer fourth-order parabolic equation

4.1. Block-centered finite difference scheme

4.2. Stability Analysis

4.3. Error Analysis

5. Time-fractional fourth-order parabolic equation

5.1. block-centered finite difference scheme

5.2. Stability Analysis

5.3. Error Analysis

6. Numerical Experiments

7. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Acknowledgments

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe