Finite Volume Element Method for Fractional Order Neutral Time-Delay Differential Equations

Abstract

Fractional-order time-delay differential equations can describe many complex physical phenomena with memory or delay effects, which are widely used in the fields of cell biology, control systems, signal processing, etc. Therefore, it is of great significance to study fractional-order time-delay differential equations. In this paper, we discuss a finite volume element method for a class of fractional-order neutral time-delay differential equations. By introducing an intermediate variable, the fourth-order problem is transformed into a system of equations consisting of two second-order partial differential equations. The L1 formula is used to approximate the time fractional order derivative terms, and the finite volume element method is used in space. A fully discrete format of the equations is established, and we prove the existence, uniqueness, convergence and stability of the solution. Finally, the validity of the format is verified by numerical examples.

Share and Cite:

Wei, Z. and Yang, Q. (2025) Finite Volume Element Method for Fractional Order Neutral Time-Delay Differential Equations. Engineering, 17, 30-52. doi: 10.4236/eng.2025.171003.

1. Introduction

The nonlinear time-delay diffusion equation is of great importance as a mathematical model for studying diffusion phenomena affected by time delay in many fields, such as the evolution of species populations, chemical reaction processes, and engineering mechanics, etc. [1]-[4].

In recent years, fractional-order differential equations have played a very important role in many disciplines, such as complex mechanics, porous media flow, hydrology and so on. The research on the numerical solution of fractional order equations has been developing rapidly. Cui et al. [5] investigated the higher-order compact finite difference format for solving one-dimensional fractional order diffusion equations, approximated the spatial second-order derivative by using compact finite difference, and discretized the fractional order derivatives of time by using the Grünwald-Letnikov formula. Gao et al. [6] established a compact finite difference format for a class of fractional-order diffusion equations, discretized the fractional-order derivatives of time using the L1 formulation, and approximated the spatial derivatives using a fourth-order accuracy compact difference. In [7], the variable order time fractional order derivatives in anomalous diffusion and wave propagation were investigated, and two second-order approximation formulas were given. Liu et al. [8] established a second-order time discretization format containing the parameter θ for nonlinear cable equations with time fractional order derivatives, which was approximated in time by the WSGD method and discretized in space by the Galerkin finite element method. In [9], Alikhanov et al. proposed the L2 1 σ formula to approximate the fractional order derivative of time based on the L12 formula, which had a convergence accuracy of the order 3α .

In some practical problems, it is necessary to use equations containing space fourth-order derivatives for mathematical modeling, such as fluid solidification process [10] [11], diffusion growth model of material surfaces [12], propagation of strong laser beams through waves [13] [14] and so on. There has been a lot of work in the study of numerical methods for fourth-order partial differential equations. In [15], Du et al. studied the finite volume element solution of fourth-order partial differential equations on a class of smooth surfaces. In [16], Liu et al. studied a mixed finite element method for time fractional order fourth-order partial differential equations. By introducing an auxiliary variable, the fourth-order equations are decomposed into a coupled system of two second-order equations. In [17], Liu and Yang proposed a finite-difference/finite-element method for a class of time-fractional order fourth-order reaction-diffusion equations with nonlinear reaction terms. The finite-difference approximation is used in the time direction, and the finite-element method is used for approximation in the spatial direction. Zhang et al. [18] proposed a fully discrete numerical format based on the linear finite-volume method for a class of time-fractional order diffusion equations, which is spatially discretized by the finite-volume method and temporally approximated by L1 interpolation. In [19], Huang et al. studied a mixed finite element method for a class of time fractional order biharmonic equations.

With the in-depth study of the anomalous diffusion phenomenon of matter in nature, it is found that the anomalous diffusion phenomenon in nature will have the effect of time delay, so there are many studies on the fractional-order nonlinear time-delay diffusion equations, such as fractional-order time-delay infectious disease model, fractional-order system dynamics [20] [21], etc. In recent years, more and more attention has been paid to the research on numerical methods of fractional-order nonlinear time-delay differential equations. In [22], Nandal et al. studied the numerical solution of a class of one-dimensional fourth-order time fractional-order nonlinear neutral time-delay differential equations, and constructed a second-order convergent differential format in the time fractional-order by approximating the fractional-order derivative term with the L2 1 σ formula. In [23], Nandal et al. studied a class of one-dimensional fourth-order fractional-order nonlinear subdiffusion equations with time delay and variable coefficients. Constructing a linearized compact difference format for the equations, approximated in time by the L2 1 σ formula and discretized in space by a compact linear operator. Hendy et al. [24] constructed a compact difference format for the multinomial time fractional-order convection-diffusion wave equation by using an interpolation method for the nonlinear source term of the time delay. Wang et al. [25] investigated a compact difference method for a class of two-dimensional time fractional order neutral time-delay diffusion equation. Xie et al. [26] studied a class of time-fractional nonlinear time-delay fourth-order diffusion equations. The time-fractional derivative is approximated by L2 1 σ , the space fourth derivative is approximated by compact difference method, and the nonlinear source term with time-delay is treated by extrapolation method.

At present, the numerical methods for solving fractional order time-delay differential equations mainly focus on the finite difference method, and the finite volume element method is rarely used. The finite volume element method is more flexible and can maintain the conservation of local physical quantities, and can achieve higher convergence accuracy than the finite difference method. Therefore, in this paper, we consider using the finite volume element method to numerically analyze the following two-dimensional fractional-order nonlinear neutral time-delay fourth-order diffusion equation

0 c D t α u( x,y,t )+ Δ 2 u( x,y,t )+ Δ 2 u( x,y,ts ) =f( x,y,t,u( x,y,t ),u( x,y,ts ) ),( x,y,t )Ω×( 0,T ]. (1)

The initial and boundary conditions are:

{ u( x,y,t )=0, ( x,y,t )Ω×[ 0,T ], Δu( x,y,t )=0, ( x,y,t )Ω×[ 0,T ], u( x,y,t )=ω( x,y,t ), ( x,y,t )Ω×[ s,0 ], (2)

where s>0 is the delay quantity, Δ 2 u( x,y,ts ) is the neutral time-delay term, f( x,y,t,u( x,y,t ),u( x,y,ts ) ) is nonlinear source term with time delay, Ω=( 0, L 1 )×( 0, L 2 ) , α( 0,1 ) , f and ω are all sufficiently smooth functions.

The Caputo time fractional order derivative is defined as follows

0 c D t α u( x,y,t )= 1 Γ( 1α ) 0 t ( tξ ) α u( x,y,ξ ) ξ dξ,0<α<1. (3)

Assume that the partial derivatives f μ ( x,y,t,μ,ν ), f ν ( x,y,t,μ,ν ) are continuous in a ϵ 0 neighborhood of μ,ν , where ϵ 0 is a positive constant. Denote

c 1 = sup 0<x< L 1 ,0<y< L 2 ,0<t<T,| ϵ 1 |,| ϵ 2 |< ϵ 0 | f μ ( x,y,t,u( x,y,t )+ ϵ 1 ,u( x,y,ts )+ ϵ 2 ) |, (4)

c 2 = sup 0<x< L 1 ,0<y< L 2 ,0<t<T,| ϵ 1 |,| ϵ 2 |< ϵ 0 | f ν ( x,y,t,u( x,y,t )+ ϵ 1 ,u( x,y,ts )+ ϵ 2 ) |. (5)

This paper focuses on the finite volume element method for the problem (1)-(2). In section 2, the fourth-order problem is transformed into a second-order system of equations by introducing an intermediate variable. The fully discrete finite volume element format of the problem is constructed by using the L1 formula to approximate the time-fractional order derivative term and using the finite volume element method to discretize the spatial derivative term. In section 3, the existence, uniqueness, convergence and stability of the finite volume element solution are proved. In section 4, specific numerical examples are given to verify the accuracy and validity of the format.

2. Fully Discrete Finite Volume Element Numerical Format

Let T h be a quasi-uniformly regular triangulation of Ω , and the triangular element in T h is denoted by K . Ω h denotes the set of all nodes in T h , and Ω h * = Ω h \Ω is the set of all inner points. Let h K be the diameter of the element K , h= max K T h { h K } .

For any P 0 Ω h * , take the adjacent nodes of P 0 to be P i ( i=1,2,,n ) , and the midpoint of P 0 P i to be M i . Let the barycenter of the triangular element Δ P 0 P i P i+1 ( P 1 = P n+1 ) be Q i and connect M 1 , Q 1 ,, M n , Q n , M 1 to get the polygonal region K P 0 * surrounding P 0 . Similarly, for any P 0 Ω , we can get the corresponding semi-enclosed polygonal region K P 0 * . The polygonal region K P 0 * is called a dual element, and all the dual elements form the barycenter dual subdivision of Ω , denoted by T h * .

2.1. Trial Function Space and Test Function Space

The trial function space U h is chose as the standard linear element space on T h , then U h ={ u h C( Ω ): u h | K P 1 ; u h | Ω =0 } H 0 1 ( Ω ) , where P 1 is the set of all linear polynomials on the triangle element K .

The test function space V h is chosed as the piecewise constant space on T h * , then V h ={ v h L 2 ( Ω ): P 0 Ω h * , v h | K P 0 * =C; P 0 Ω, v h | K P 0 * =0 } . Thus, V h can be viewed as a space spanned by the following basis functions: for P 0 Ω h * ,

ϕ P 0 ( P )={ 1, P K P 0 * , 0, elsewhere. (6)

2.2. Finite Volume Element Format

Firstly, the intermediate function v( x,y,t )=Δu( x,y,t ) is introduced, then the original equation can be rewritten as

{ 0 c D t α u( x,y,t )Δv( x,y,t )Δv( x,y,ts ) =f( x,y,t,u( x,y,t ),u( x,y,ts ) ), ( x,y,t )Ω×( 0,T ], v( x,y,t )=Δu( x,y,t ), ( x,y,t )Ω×( 0,T ], u( x,y,t )=0, ( x,y,t )Ω×[ 0,T ], v( x,y,t )=0, ( x,y,t )Ω×[ 0,T ], u( x,y,t )=ω( x,y,t ), ( x,y,t )Ω×[ s,0 ], v( x,y,t )=Δω( x,y,t ), ( x,y,t )Ω×[ s,0 ]. (7)

Multiplying ϕ h V h and ψ h V h at both sides of the first equation and the second equation of (7), respectively and utilizing Green’s formula to get:

{ ( 0 c D t α u, ϕ h )+a( v, ϕ h )+a( v( x,y,ts ), ϕ h ) =( f( x,y,t,u,u( x,y,ts ) ), ϕ h ), ϕ h V h , ( v, ψ h )=a( u, ψ h ), ψ h V h , (8)

where the bilinear form a( v, ϕ h )= P 0 Ω h ϕ h ( P 0 ) K P 0 * v x dy v y dx .

Choosing the trial function space U h as described above, and we can obtain the semi-discrete finite volume element format: find ( u h ( ,t ), v h ( ,t ) ) U h × U h ( 0tT ) such that

{ ( 0 c D t α u h , ϕ h )+a( v h , ϕ h )+a( v h ( ,ts ), ϕ h ) =( f( x,y,t, u h , u h ( x,y,ts ) ), ϕ h ), ϕ h V h , ( v h , ψ h )=a( u h , ψ h ), ψ h V h . u h ( x,y,t )=0, ( x,y,t )Ω×[ 0,T ], v h ( x,y,t )=0, ( x,y,t )Ω×[ 0,T ], u h ( x,y,t )= ω h ( x,y,t ), ( x,y,t )Ω×[ s,0 ], v h ( x,y,t )= ϖ h ( x,y,t ), ( x,y,t )Ω×[ s,0 ], (9)

where ω h ( x,y,t ) and ϖ h ( x,y,t ) are the elliptic projections of ω( x,y,t ) and Δω( x,y,t ) on U h , respectively.

Next, we discrete time fractional order derivatives and deal with nonlinear terms. Let m be a positive integer, and make τ= s m , N= T τ , t n =nτ( n=0,1,2,,N ) . For convenience, assume that N is a positive integer.

Lemma 1 [27] Let 0<α<1, a l α = ( l+1 ) 1α l 1α ,l=0,1,2, , then

(I) 1= a 0 α > a 1 α > a 2 α >> a l α >0; a l α 0,ifl, (10)

(II) ( 1α ) l α < a l1 α <( 1α ) ( l1 ) α ,l1. (11)

Lemma 2 [27] For α( 0<α<1 ) order Caputo fractional order derivatives, there are

0 c D t α f( t n )= τ α Γ( 2α ) [ a 0 α f( t n ) k=1 n1 ( a nk1 α a nk α )f( t k ) a n1 α f( t 0 ) ]+ R n , (12)

where | R n | 1 2Γ( 1α ) [ 1 4 + α ( 1α )( 2α ) ] max t 0 t t n | f tt | τ 2α .

When t= t n , for the right nonlinear term, Taylor’s formula yields:

f( x,y,t, u h ( x,y, t n ), u h ( x,y, t nm ) ) =f( x,y,t,2 u h ( x,y, t n1 ) u h ( x,y, t n2 ), u h ( x,y, t nm ) )+O( τ 2 ). (13)

Substitute (13) into (9). The applying the Lemma 2 to the fractional derivative term and omitting small quantity terms, we get the fully discrete finite volume element format: find ( u h n , v h n ) U h × U h ( n=1,2, ) such that

{ ( D τ α u h n , ϕ h )+a( v h n , ϕ h )+a( v h nm , ϕ h ) =( f( x,y, t n ,2 u h n1 u h n2 , u h nm ), ϕ h ), ϕ h V h , ( v h n , ψ h )=a( u h n , ψ h ), ψ h V h , u h n =0, ( x,y )Ω, v h n =0, ( x,y )Ω, u h l = ω h l , l=m,,1,0, v h l = ϖ h l , l=m,,1,0, (14)

where D τ α u h n = τ α Γ( 2α ) [ a 0 α u h n k=1 n1 ( a nk1 α a nk α ) u h k a n1 α u h 0 ] .

For the purpose of subsequent research, the following definitions and lemmas are presented.

Definition 1. [28] The interpolation operator Π * is defined as follows:

Π * : H 0 1 V h , Π * u= P 0 Ω h * u( P 0 ) ϕ P 0 ,u H 0 1 . (15)

Definition 2. [28] Introduce the elliptic projection operator P h : H 2 ( Ω ) H 0 1 ( Ω ) U h , defined by the following format:

a( P h u, v h )=a( u, v h ), v h V h . (16)

Lemma 3. [28] Let u W 3,p ( Ω ) H 0 1 ( Ω ) , and P h u is the elliptic projection of u , then

u P h u C h 2 u 3,p . (17)

Lemma 4. [18] [28] Let | u h |= ( u h , Π * u h ) 1 2 , we have

(I) ( u h , Π * w h )=( w h , Π * u h ), u h , w h U h , (18)

(II) ( u h , Π * u h )0, u h U h , (19)

(III) | ( u h , Π * w h ) || u h || w h |, u h , w h U h . (20)

Lemma 5. [18] [28] | | and are equivalent on U h , i.e., there exist constants C 1 , C 2 >0 such that

C 1 u h | u h | C 2 u h , u h U h . (21)

Lemma 6. [28] There exist positive constants α and M such that the following equations hold

a( u h , Π * w h )=a( w h , Π * u h ), u h , w h U h , (22)

a( u h , Π * u h )α u h 1 2 , u h U h , (23)

| a( u h , Π * w h ) |M u h 1 w h 1 , u h , w h U h . (24)

Lemma 7. [29] Let { z k },{ g k } be two non-negative sequences. If

z k K+ i=0 k g i z i ,k0, (25)

then

z k Kexp( i=0 k g i ),k0, (26)

where K is a positive constant.

3. Theoretical Analysis of the Format

3.1. Uniqueness of the Solution

Theorem 1. The fully discrete finite volume element format (14) has a unique solution.

Proof. The homogeneous equations corresponding to (14) are

{ ( u h n , ϕ h )+μa( v h n , ϕ h )=0, ϕ h V h , ( v h n , ψ h )=a( u h n , ψ h ), ψ h V h , (27)

where μ= τ α Γ( 2α ) .

To prove that (14) is uniquely solvable, it is sufficient to prove that its corresponding homogeneous equations (27) has only zero solutions. Let ϕ h = Π * v h n , ψ h = Π * u h n in (27), then

( u h n , Π * v h n )+μa( v h n , Π * v h n )=0, (28)

( v h n , Π * u h n )=a( u h n , Π * u h n ). (29)

According to Lemma 4

( v h n , Π * u h n )=( u h n , Π * v h n ). (30)

Substituting (30) into (29) gives

( u h n , Π * v h n )=a( u h n , Π * u h n ). (31)

Substituting (31) into (28) gives

a( u h n , Π * u h n )+μa( v h n , Π * v h n )=0. (32)

According to Lemma 6

a( u h n , Π * u h n )α u h n 1 2 0,a( v h n , Π * v h n )α v h n 1 2 0. (33)

We get

u h n 1 = v h n 1 =0. (34)

So u h n =0, v h n =0 . According to the induction principle, the homogeneous equations has only zero solutions, so the format (14) is uniquely solvable.

3.2. Convergence Analysis

Theorem 2. Assuming that { u,v } and { u h n , v h n } are solutions of the original Equation (7) and the finite volume element format (14), respectively, and that the requirement of sufficient regularity is satisfied, then there exists a constant C such that

u n u h n C( h 2 + τ 2α ), (35)

v n v h n C( h 2 + τ 2α ). (36)

Proof. Let

η n = u n P h u n , ε n = P h u n u h n , (37)

θ n = v n P h v n , ρ n = P h v n v h n , (38)

then

u n u h n = η n + ε n , (39)

v n v h n = θ n + ρ n . (40)

Let t= t n in (8) and subtract from (14). According to the elliptic projection in Definition 2, we get

( D τ α η n + D τ α ε n , ϕ h )+a( ρ n , ϕ h )+a( ρ nm , ϕ h ) =( f( x,y, t n ,2 u n1 u n2 , u nm )f( x,y, t n ,2 u h n1 u h n2 , u h nm ), ϕ h ) +( R 1 , ϕ h ), ϕ h V h , (41)

( θ n + ρ n , ψ h )=a( ε n , ψ h ), ψ h V h , (42)

where R 1 c 0 τ 2α and c 0 is a positive constant that does not depend on τ,h .

Let ϕ h = Π * ε n and ψ h = Π * ρ n in (41) and (42), respectively, and use Lemma 6 to get

( D τ α η n + D τ α ε n , Π * ε n )+a( ρ n , Π * ε n )+a( ρ nm , Π * ε n ) =( f( x,y, t n ,2 u n1 u n2 , u nm )f( x,y, t n ,2 u h n1 u h n2 , u h nm ), Π * ε n ) +( R 1 , Π * ε n ), (43)

a( ρ n , Π * ε n )=a( ε n , Π * ρ n )=( θ n , Π * ρ n )+( ρ n , Π * ρ n ). (44)

Let ψ h = Π * ρ nm in (42), and use Lemma 6 to get

a( ρ nm , Π * ε n )=a( ε n , Π * ρ nm )=( θ n , Π * ρ nm )+( ρ n , Π * ρ nm ). (45)

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

( D τ α ε n , Π * ε n )+( ρ n , Π * ρ n ) =( D τ α η n , Π * ε n )( θ n , Π * ρ n )( θ n , Π * ρ nm )( ρ n , Π * ρ nm ) +( f( x,y, t n ,2 u n1 u n2 , u nm )f( x,y, t n ,2 u h n1 u h n2 , u h nm ), Π * ε n ) +( R 1 , Π * ε n ). (46)

From the definition of D τ α ε n , multiplying μ= τ α Γ( 2α ) on both sides of the (46) to get

( ε n , Π * ε n )+μ( ρ n , Π * ρ n ) =μ( D τ α η n , Π * ε n )μ( θ n , Π * ρ n )μ( θ n , Π * ρ nm )μ( ρ n , Π * ρ nm ) +μ( f( x,y, t n ,2 u n1 u n2 , u nm )f( x,y, t n ,2 u h n1 u h n2 , u h nm ), Π * ε n ) +μ( R 1 , Π * ε n )+ k=1 n1 ( a nk1 α a nk α )( ε k , Π * ε n ). (47)

In order to estimate the order of convergence of the solution u h n , we analyze the terms at the right side of the equation (47). According to the properties of the function f , it follows that

| f( x,y, t n ,2 u n1 u n2 , u nm )f( x,y, t n ,2 u h n1 u h n2 , u h nm ) | c 1 | ( 2 u n1 u n2 )( 2 u h n1 u h n2 ) |+ c 2 | u nm u h nm | c 1 | 2( η n1 + ε n1 )( η n2 + ε n2 ) |+ c 2 | η nm + ε nm | 2 c 3 | η n1 |+2 c 3 | ε n1 |+ c 3 | η n2 |+ c 3 | ε n2 |+ c 3 | η nm |+ c 3 | ε nm |, (48)

where c 3 =max{ c 1 , c 2 } .

According to the Cauchy-Schwarz inequality and (48), we get

| μ( f( x,y, t n ,2 u n1 u n2 , u nm )f( x,y, t n ,2 u h n1 u h n2 , u h nm ), Π * ε n ) | 1 48 c 3 2 μ | f( x,y, t n ,2 u n1 u n2 , u nm )f( x,y, t n ,2 u h n1 u h n2 , u h nm ) | 2 +12 c 3 2 μ | ε n | 2 1 48 c 3 2 ( 4 c 3 2 μ | η n1 | 2 +4 c 3 2 μ | ε n1 | 2 + c 3 2 μ | η n2 | 2 + c 3 2 μ | ε n2 | 2 + c 3 2 μ | η nm | 2 + c 3 2 μ | ε nm | 2 )+12 c 3 2 μ | ε n | 2 1 12 μ | η n1 | 2 + 1 12 μ | ε n1 | 2 + 1 48 μ | η n2 | 2 + 1 48 μ | ε n2 | 2 + 1 48 μ | η nm | 2 + 1 48 μ | ε nm | 2 +12 c 3 2 μ | ε n | 2 . (49)

For the remaining terms at the right side, using the Cauchy-Schwartz inequality and substituting (49) into (47) to obtain

( ε n , Π * ε n )+μ( ρ n , Π * ρ n ) 1 2 μ | D τ α η n | 2 + 3 4 μ | ε n | 2 +3μ | θ n | 2 +μ | ρ n | 2 + 11 24 μ | ρ nm | 2 + 1 12 μ | η n1 | 2 + 1 12 μ | ε n1 | 2 + 1 48 μ | η n2 | 2 + 1 48 μ | ε n2 | 2 + 1 48 μ | η nm | 2 + 1 48 μ | ε nm | 2 +12 c 3 2 μ | ε n | 2 +μ | R 1 | 2 + k=1 n1 2( a nk1 α a nk α ) | ε k | 2 + 1 8 | ε n | 2 . (50)

Then, we can get

( 1 3 4 μ12 c 3 2 μ 1 8 ) | ε n | 2 1 2 μ | D τ α η n | 2 +3μ | θ n | 2 + 11 24 μ | ρ nm | 2 + 1 12 μ | η n1 | 2 + 1 12 μ | ε n1 | 2 + 1 48 μ | η n2 | 2 + 1 48 μ | ε n2 | 2 + 1 48 μ | η nm | 2 + 1 48 μ | ε nm | 2 +μ | R 1 | 2 + k=1 n1 2( a nk1 α a nk α ) | ε k | 2 . (51)

According to the definition of μ , and take a sufficiently small τ to get

| ε n | 2 4μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 ) + 11 18 μ | ρ nm | 2 + 1 9 μ | ε n1 | 2 + 1 36 μ | ε n2 | 2 + 1 36 μ | ε nm | 2 + k=1 n1 8 3 ( a nk1 α a nk α ) | ε k | 2 . (52)

(52) will be used in the subsequent proof to estimate the convergence order of the solution u h n .

In order to estimate the order of convergence of the solution v h n , the terms at the right side of (47) are analyzed. Similar to the derivation of (52), we get

μ | ρ n | 2 6μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 ) + 9 4 μ | ρ nm | 2 + 1 6 μ | ε n1 | 2 + 1 24 μ | ε n2 | 2 + 1 24 μ | ε nm | 2 + k=1 n1 ( a nk1 α a nk α ) | ε k | 2 . (53)

(53) will be used in the subsequent proof to estimate the convergence order of the solution v h n .

Let d= T s . If d is an integer, then the time period T is divided into d segments, where the length of each segment is s . If d is not an integer, then the time period T is divided into [ d ]+1 segments, where the length of the first [ d ] segments are s and the length of the last segment is T[ d ]s . Since m= s τ , then n=1,2,,m on the first segment, n=m+1,m+2,,2m on the second segment, and so on. Now let’s first prove that (35)-(36) are true on the first segment.

Starting from (52) we estimate the order of convergence of the solution u h n . The function is known at [ s,0 ] and n=1,2,,m , so ρ nm =0 , ε nm =0 on the first segment. Substitute ρ nm =0 , ε nm =0 into (52) to get

| ε n | 2 4μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 ) + 1 9 μ | ε n1 | 2 + 1 36 μ | ε n2 | 2 + k=1 n1 8 3 ( a nk1 α a nk α ) | ε k | 2 . (54)

By Lemma 1 and Lemma 7, it follows that

| ε n | 2 4μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 )exp( 1 9 μ+ 1 36 μ+ k=1 n1 8 3 ( a nk1 α a nk α ) ) 4exp( 3 )μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 ). (55)

According to the equivalence of norms in Lemma 5 and the property of elliptic projection in Lemma 3,

| θ n |C h 2 max 0tT v 3,p c 4 h 2 , (56)

| η nk |C h 2 max 0tT u 3,p c 5 h 2 ,k=1,2,m. (57)

For D τ α η n ,

| D τ α η n |=| τ α Γ( 2α ) k=1 n a nk α δ t η k | =| τ α Γ( 2α ) k=1 n a nk α t k1 t k ( I P h ) u t dt | τ α Γ( 2α ) k=1 n a nk α τ max t 0 t t n | ( I P h ) u t | (nτ) 1α Γ( 2α ) max t 0 t t n | ( I P h ) u t | T 1α Γ( 2α ) max t 0 t t n | ( I P h ) u t |. (58)

According to (58) and Lemma 3, it follows that

| D τ α η n | T 1α Γ( 2α ) C h 2 max 0tT u t 3,p c 6 h 2 . (59)

Substituting (56), (57) and (59) into (55), we get

| ε n | 2 4exp( 3 )μ( c 4 2 h 4 + c 5 2 h 4 + c 5 2 h 4 + c 5 2 h 4 + c 6 2 h 4 + c 0 2 ( τ 2α ) 2 ). (60)

Let C 1 =max{ c 4 2 +3 c 5 2 + c 6 2 + c 6 2 , c 0 2 } , and then

| ε n | 2 4exp( 3 ) C 1 μ( h 4 + ( τ 2α ) 2 ). (61)

Using the equivalence of norms in Lemma 5 and the definition of μ to get

ε n C( h 2 + τ 2α ),n=1,2,,m. (62)

Finally, it can be proved that (35) holds on the first segment by using the triangle inequality.

Then starting from (53) we estimate the order of convergence of the solution v h n . Substituting ρ nm =0 , ε nm =0 into (53) yields

μ | ρ n | 2 6μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 ) + 1 6 μ | ε n1 | 2 + 1 24 μ | ε n2 | 2 + k=1 n1 ( a nk1 α a nk α ) | ε k | 2 . (63)

The right side terms of (63) are analyzed below. From (56), (57) and (59), we can obtain

6μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 ) 6 C 1 μ( h 4 + ( τ 2α ) 2 ). (64)

From (61), we can obtain

1 6 μ | ε n1 | 2 + 1 24 μ | ε n2 | 2 5 6 exp( 3 ) C 1 μ( h 4 + ( τ 2α ) 2 ), (65)

k=1 n1 ( a nk1 α a nk α ) | ε k | 2 k=1 n1 ( a nk1 α a nk α )4exp( 3 ) C 1 μ( h 4 + ( τ 2α ) 2 ) 4exp( 3 ) C 1 μ( h 4 + ( τ 2α ) 2 ). (66)

Substituting (64)-(66) into (63) and dividing both sides by μ to get

| ρ n | 2 ( 6+ 5 6 exp( 3 )+4exp( 3 ) ) C 1 ( h 4 + ( τ 2α ) 2 ) 6exp( 3 ) C 1 ( h 4 + ( τ 2α ) 2 ). (67)

Using the equivalence of the norms in Lemma 5, we can obtain

ρ n C( h 2 + τ 2α ),n=1,2,,m. (68)

Finally, by using the triangle inequality, it can be proved that (36) holds on the first segment.

Continue to prove that (35) (36) hold on the second segment. Starting from (52), we get

| ε n | 2 4μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 + | ρ nm | 2 )+ 1 9 μ | ε n1 | 2 + 1 36 μ | ε n2 | 2 + 1 36 μ | ε nm | 2 + k=1 n1 8 3 ( a nk1 α a nk α ) | ε k | 2 . (69)

According to lemma 7,

| ε n | 2 4μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 + | ρ nm | 2 )exp( 1 9 μ+ 1 36 μ++ 1 36 μ+ k=1 n1 8 3 ( a nk1 α a nk α ) ) 4exp( 3 )μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 + | ρ nm | 2 ), (70)

From (67),

| ρ nm | 2 6exp( 3 ) C 1 ( h 4 + ( τ 2α ) 2 ). (71)

Let C 2 =max{ c 4 2 +3 c 5 2 + c 6 2 +6exp( 3 ) C 1 , c 0 2 +6exp( 3 ) C 1 } . According to the properties of elliptic projection and (71), we can get

| ε n | 2 4exp( 3 ) C 2 μ( h 4 + ( τ 2α ) 2 ). (72)

Using the equivalence of the norms in Lemma 5 and the definition of μ

ε n C( h 2 + τ 2α ),n=m+1,m+2,,2m. (73)

Finally, by using the triangle inequality, it can be proved that (35) holds on the second segment.

Next, from (53) we can get

μ | ρ n | 2 6μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 + | ρ nm | 2 )+ 1 6 μ | ε n1 | 2 + 1 24 μ | ε n2 | 2 + 1 24 μ | ε nm | 2 + k=1 n1 ( a nk1 α a nk α ) | ε k | 2 . (74)

The right side terms of (74) are analyzed below. From (61) and (67), we get

| ε nm | 2 4exp( 3 ) C 1 μ( h 4 + ( τ 2α ) 2 ), (75)

| ρ nm | 2 6exp( 3 ) C 1 ( h 4 + ( τ 2α ) 2 ). (76)

According to the properties of elliptical projection and (76), it follows that

6μ( | D τ α η n | 2 + | θ n | 2 + | η n1 | 2 + | η n2 | 2 + | η nm | 2 + | R 1 | 2 + | ρ nm | 2 ) 6 C 2 μ( h 4 + ( τ 2α ) 2 ). (77)

From (72) and (75), we have

1 6 μ | ε n1 | 2 + 1 24 μ | ε n2 | 2 + 1 24 μ | ε nm | 2 exp( 3 ) C 2 μ( h 4 + ( τ 2α ) 2 ). (78)

From (72), we get

k=1 n1 ( a nk1 α a nk α ) | ε k | 2 k=1 n1 ( a nk1 α a nk α )4exp( 3 ) C 2 μ( h 4 + ( τ 2α ) 2 ) 4exp( 3 ) C 2 μ( h 4 + ( τ 2α ) 2 ). (79)

Substituting (77), (78), (79) into (74) and dividing both sides by μ , we get

| ρ n | 2 ( 6+exp( 3 )+4exp( 3 ) ) C 2 ( h 4 + ( τ 2α ) 2 ) 6exp( 3 ) C 2 ( h 4 + ( τ 2α ) 2 ). (80)

Using the equivalence of the norms in Lemma 5, we can obtain

ρ n C( h 2 + τ 2α ),n=m+1,m+2,,2m. (81)

Finally, by using the triangle inequality, it can be proved that (36) holds on the second segment.

Since there are a finite number of segments, by analogy, we can prove that Theorem 2 holds at each segment. The proof is completed.

3.3. Stability Analysis

Let { p h n , q h n } be the solution of the following equations

{ ( D τ α p h n , ϕ h )+a( q h n , ϕ h )+a( q h nm , ϕ h ) =( f( x,y,t,2 p h n1 p h n2 , p h nm ), ϕ h ), ϕ h V h , ( q h n , ψ h )=a( p h n , ψ h ), ψ h V h , p h n =0, ( x,y )Ω, q h n =0, ( x,y )Ω, p h l = ω h l + f 1 l , l=m,,1,0, q h l = ϖ h l + f 2 l , l=m,,1,0, (82)

where f 1 , f 2 are the perturbation terms of the initial function.

Let θ n = u h n p h n , λ n = v h n q h n ,n=1,2,,N .

Theorem 3 There exists a constant C such that

θ n 2 C( max ml0 | f 1 l | 2 + max ml0 | f 2 l | 2 ), (83)

λ n 2 C( max ml0 | f 1 l | 2 + max ml0 | f 2 l | 2 ). (84)

Proof. Subtracting (82) from (14), we get

( D τ α θ n , ϕ h )+a( λ n , ϕ h )+a( λ nm , ϕ h ) =( f( x,y,t,2 u h n1 u h n2 , u h nm )f( x,y,t,2 p h n1 p h n2 , p h nm ), ϕ h ), ϕ h V h , (85)

( λ n , ψ h )=a( θ n , ψ h ), ψ h V h . (86)

According to the derivation of (52) and (53) in the convergence analysis, and similarly for (85) and (86), we can get

| θ n | 2 1 3 μ | λ nm | 2 + 1 9 μ | θ n1 | 2 + 1 36 μ | θ n2 | 2 + 1 36 μ | θ nm | 2 + k=1 n1 8 3 ( a nk1 α a nk α ) | θ k | 2 + 8 3 a n1 α | θ 0 | 2 , (87)

μ | λ n | 2 μ | λ nm | 2 + 1 6 μ | θ n1 | 2 + 1 24 μ | θ n2 | 2 + 1 24 μ | θ nm | 2 + k=1 n1 ( a nk1 α a nk α ) | θ k | 2 + a n1 α | θ 0 | 2 . (88)

Similar to the proof procedure of Theorem 2, the stability result can be obtained.

4. Numerical Examples

In this section, two numerical examples will be given to verify the validity of the finite volume element format. Let Ω= ( 0,1 ) 2 , T=1 , and define

E L 2 ( h,τ )= u N u h N , E ( h,τ )= u N u h N . (89)

The convergence order is calculated as follows:

order( h )= log 2 ( E L 2 ( 2h,τ ) E L 2 ( h,τ ) ),order( τ )= log 2 ( E L 2 ( h,2τ ) E L 2 ( h,τ ) ). (90)

Example 1. Consider the exact solution of problem (1) as:

u= ( t+1 ) 4+α sin( πx )sin( πy ),( x,y,t )( x,y,t )Ω×( 0,T ], (91)

The nonlinear source term with time delay is

f=u ( x,y,t ) 2 +u( x,y,ts )+G( x,y,t ), (92)

where

G( x,y,t )= Γ( 5+α ) Γ( 5 ) ( t+1 ) 4 sin( πx )sin( πy )+4 π 4 ( t+1 ) 4+α sin( πx )sin( πy ) +4 π 4 ( t+10.4 ) 4+α sin( πx )sin( πy ) ( ( t+1 ) 4+α sin( πx )sin( πy ) ) 2 ( t+10.4 ) 4+α sin( πx )sin( πy ). (93)

Taking s=0.4 , the above problem can be solved by using the finite volume element format constructed in this paper. When α=0.3 , we fix the time step τ= 1 40 , and change the spatial step to get the errors and the spatial convergence order of u h , v h . We fix the spatial step h= 1 120 , and change the time step to get the errors and the time convergence order of u h , v h . The results are given in Table 1 and Table 2, respectively.

From the data in Table 1 and Table 2, it can be seen that the spatial convergence order in L -normand L 2 -norm are approximately equal to 2, and the temporal convergence order in L -normand L 2 -norm are close to 1.7. This is consistent with the theoretical analysis.

Table 1. Errors and spatial convergence order of finite volume element solutions u h , v h when α=0.3 .

h

u u h

order

u u h

order

v v h

order

v v h

order

1 20

2.087e−01

9.456e−02

2.483

1.102

1 30

9.305e−02

1.995

4.204e−02

1.994

1.102

2.001

4.882e−01

2.004

1 40

5.224e−02

1.998

2.357e−02

2.004

6.180e−01

2.006

2.728e−01

2.014

1 50

3.329e−02

2.011

1.501e−02

2.012

3.932e−01

2.016

1.731e−01

2.024

Table 2. Errors and time convergence order of finite volume element solutions u h , v h when α=0.3 .

τ

u u h

order

u u h

order

v v h

order

v v h

order

1 5

1.338e−01

6.550e−02

2.759

1.271

1 10

4.157e−02

1.686

2.011e−02

1.703

8.244e−01

1.746

3.743e−01

1.764

1 20

1.363e−02

1.608

6.393e−03

1.653

2.356 e−01

1.803

1.032e−01

1.858

When α=0.5 , we fix the time step τ= 1 40 , and change the spatial step to get the errors and spatial convergence order of u h , v h . We fix the spatial step h= 1 90 , and change the time step to get the errors and temporal convergence order of u h , v h . The results are given by Table 3 and Table 4, respectively.

Table 3. Errors and order of spatial convergence of finite volume element solutions u h , v h when α=0.5 .

h

u u h

order

u u h

order

v v h

order

v v h

order

1 20

2.449e−01

1.110e−01

2.964

1.309

1 30

1.103e−01

1.995

4.992e−02

1.994

1.341

2.001

5.902e−01

2.004

1 40

6.286e−02

1.974

2.842e−02

1.972

7.705e−01

1.965

3.378e−01

1.970

1 50

4.081e−02

1.955

1.846e−02

1.953

5.060e−01

1.921

2.209e−01

1.933

Table 4. Errors and time convergence order of finite volume element solutions u h , v h when α=0.5 .

τ

u u h

order

u u h

order

v v h

order

v v h

order

1 5

1.985e−01

9.699e−02

4.063

1.870

1 10

6.604e−02

1.587

3.182e−02

1.607

1.280

1.665

5.815e−01

1.685

1 20

2.502e−02

1.401

1.170e−02

1.443

4.190e−01

1.611

1.841e−01

1.659

From the data in Table 3 and Table 4, it can be seen that the spatial convergence order in L -normand L 2 -norm are approximately equal to 2, and the temporal convergence order in L -normand L 2 -norm are close to 1.5. This is consistent with the theoretical analysis.

When α=0.7 , we fix the time step τ= 1 60 , and change the spatial step to get the errors and spatial convergence order of u h , v h . We fix the spatial step h= 1 70 , and change the time step to get the errors and temporal convergence order of u h , v h . The results are given by Table 5 and Table 6, respectively.

Table 5. Errors and spatial convergence order of finite volume element solutions u h , v h when α=0.7 .

h

u u h

order

u u h

order

v v h

order

v v h

order

1 20

2.825e−01

1.280e−01

3.431

1.514

1 30

1.266e−01

1.984

5.728e−02

1.982

1.538

1.983

6.767e−01

1.986

1 40

7.162e−02

1.979

3.237e−02

1.984

8.728e−01

1.974

3.830e−01

1.982

1 50

4.609e−02

1.982

2.083e−02

1.981

5.642e−01

1.969

2.469e−01

1.976

Table 6. Errors and time convergence order of finite volume element solutions u h , v h when α=0.7 .

τ

u u h

order

u u h

order

v v h

order

v v h

order

1 5

2.912e−01

1.420e−01

5.917

2.718

1 10

1.028e−01

1.501

4.938e−02

1.524

1.958

1.595

8.858e−01

1.617

1 20

4.332e−02

1.247

2.018e−02

1.290

7.057e−01

1.472

3.097e−01

1.515

From the data in Table 5 and Table 6, it can be seen that the spatial convergence order in L -normand L 2 -norm are approximately equal to 2, and the temporal convergence order in L -normand L 2 -norm are close to 1.3. This is consistent with the theoretical analysis.

The finite volume element solutions and exact solutions are shown in Figure 1 and Figure 2.

Example 2. Consider the case where the exact solution of problem (1) is unknown, with the following boundary and initial conditions:

Figure 1. Comparison of u h ,u . (a) The finite volume element solution u h ; (b) The exact solution u .

Figure 2. Comparison of v h ,v . (a) The finite volume element solution v h ; (b) The exact solution v .

{ u( x,y,t )=0, ( x,y,t )Ω×[ 0,T ], Δu( x,y,t )=0, ( x,y,t )Ω×[ 0,T ], u( x,y,t )=2 ( t+1 ) 3 sin( πx )sin( πy ), ( x,y,t )Ω×[ s,0 ], Δu( x,y,t )=8 π 4 ( t+1 ) 3 sin( πx )sin( πy ), ( x,y,t )Ω×[ s,0 ]. (94)

The nonlinear source term with time delay is

f=u ( x,y,t ) 2 +u( x,y,ts )+G( x,y,t ), (95)

where

G( x,y,t )= 2Γ( 4 ) Γ( 4α ) ( t+1 ) 3α sin( πx )sin( πy )+8 π 4 (t+1) 3 sin( πx )sin( πy ) +8 π 4 ( t+10.4 ) 3 sin( πx )sin( πy ) ( 2 ( t+1 ) 3 sin( πx )sin( πy ) ) 2 2 ( t+10.4 ) 3 sin( πx )sin( πy ). (96)

Taking s=0.4 , and we choose the numerical solution when h= 1 120 , τ= 1 60 as the reference solution. Using the finite volume element format constructed in this paper to solve the above problem, the errors, temporal and spatial convergence order of u h , v h when α=0.3 are given in Table 7 and Table 8, respectively.

From the data in Table 7 and Table 8, we can see that the spatial convergence

Table 7. Errors and spatial convergence order of finite volume element solutions u h , v h when α=0.3 .

h

u u h

order

u u h

order

v v h

order

v v h

order

1 10

3.207e−01

1.598e−02

1.102

5.393e−01

1 20

8.030e−02

1.998

4.029e−02

1.9878

2.753e−01

2.001

1.356e−01

1.991

1 40

1.845e−02

2.121

9.278e−03

2.118

6.327e−02

2.121

3.122e−02

2.119

Table 8. Errors and time convergence order of finite volume element solutions u h , v h when α=0.3 .

τ

u u h

order

u u h

order

v v h

order

v v h

order

1 5

3.756e−02

1.847e−02

7.886e−01

3.653e−01

1 10

9.360e−03

2.004

4.604e−03

2.003

1.965e−01

2.004

9.105e−02

2.004

1 20

2.354e−03

1.991

1.158e−03

1.990

4.946 e−02

1.990

2.290e−02

1.990

order in L -normand L 2 -norm are approximately equal to 2, and the time convergence order of L -normand L 2 -norm are close to 2.

When α=0.5 , the errors, temporal and spatial convergence order of u h , v h are given in Table 9 and Table 10, respectively.

From the data in Table 9 and Table 10, we can see that the spatial convergence order in L -normand L 2 -norm are approximately equal to 2, and the time convergence order in L -normand L 2 -norm are close to 2.

When α=0.7 , the errors, temporal and spatial convergence order of u h , v h are given in Table 11 and Table 12, respectively.

From the data in Table 11 and Table 12, we can see that the spatial convergence order in L -normand L 2 -norm are approximately equal to 2, and the time convergence order in L -normand L 2 -norm are close to 2.

Table 9. Errors and order of spatial convergence of finite volume element solutions u h , v h when α=0.5 .

h

u u h

order

u u h

order

v v h

order

v v h

order

1 10

3.206e−01

1.597e−01

1.099

5.377e−01

1 20

8.025e−02

1.998

4.026e−02

1.987

2.745e−01

2.001

1.352e−01

1.991

1 40

1.844e−02

2.121

9.271e−03

2.118

6.308e−02

2.121

3.113e−02

2.119

Table 10. Errors and time convergence order of finite volume element solutions u h , v h when α=0.5 .

τ

u u h

order

u u h

order

v v h

order

v v h

order

1 5

3.702e−02

1.820e−02

7.779e−01

3.601e−01

1 10

9.149e−03

2.016

4.499e−03

2.017

1.923e−01

2.015

8.897e−02

2.016

1 20

2.294e−03

1.995

1.128e−03

1.995

4.828e−02

1.994

2.231e−02

1.995

Table 11. Errors and spatial convergence order of finite volume element solutions u h , v h when α=0.7 .

h

u u h

order

u u h

order

v v h

order

v v h

order

1 10

3.204e−01

1.596e−01

1.095

5.362e−01

1 20

8.021e−02

1.998

4.024e−02

1.987

2.737e−01

2.001

1.348e−01

1.991

1 40

1.843e−02

2.121

9.267e−03

2.118

6.289e−02

2.121

3.104e−02

2.119

Table 12. Errors and time convergence order of finite volume element solutions u h , v h when α=0.7 .

τ

u u h

order

u u h

order

v v h

order

v v h

order

1 5

3.601e−02

1.770e−02

7.579e−01

3.501e−01

1 10

8.822e−03

2.029

4.335e−03

2.030

1.859e−01

2.026

8.574e−02

2.030

1 20

2.136e−03

2.045

1.049e−03

2.046

4.514e−02

2.042

2.075e−02

2.046

5. Conclusion

In this paper, a finite volume element method for solving fractional-order neutral time-delay differential equations is proposed. By introducing an intermediate variable, the fourth-order problem is reduced to a second-order system of equations, and the fully discrete finite volume element format of the problem is constructed. We prove the existence, uniqueness, convergence and stability of the solution. Specific numerical examples are given to verify the validity of the format. In this paper, the finite volume element method is used to solve the fractional order time-delay problem in a novel way, and satisfactory theoretical results are obtained, which provide a new way for the subsequent research on the numerical method of fractional order time-delay differential equations.

Conflicts of Interest

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

References

[1] Rezounenko, A.V. and Wu, J. (2006) A Non-Local PDE Model for Population Dynamics with State-Selective Delay: Local Theory and Global Attractors. Journal of Computational and Applied Mathematics, 190, 99-113.
https://doi.org/10.1016/j.cam.2005.01.047
[2] Beta, C., Bertram, M., Mikhailov, A.S., Rotermund, H.H. and Ertl, G. (2003) Controlling Turbulence in a Surface Chemical Reaction by Time-Delay Autosynchronization. Physical Review E, 67, Article ID: 046224.
https://doi.org/10.1103/physreve.67.046224
[3] Li, D., Zhang, C. and Wang, W. (2012) Long Time Behavior of Non-Fickian Delay Reaction-Diffusion Equations. Nonlinear Analysis: Real World Applications, 13, 1401-1415.
https://doi.org/10.1016/j.nonrwa.2011.11.005
[4] Benchohra, M. and Berhoun, F. (2010) Impulsive Fractional Differential Equations with State-Dependent Delay. Communications in Applied Analysis, 14, 213-224.
[5] Cui, M. (2009) Compact Finite Difference Method for the Fractional Diffusion Equation. Journal of Computational Physics, 228, 7792-7804.
https://doi.org/10.1016/j.jcp.2009.07.021
[6] Gao, G. and Sun, Z. (2011) A Compact Finite Difference Scheme for the Fractional Sub-Diffusion Equations. Journal of Computational Physics, 230, 586-595.
https://doi.org/10.1016/j.jcp.2010.10.007
[7] Zhao, X., Sun, Z. and Karniadakis, G.E. (2015) Second-Order Approximations for Variable Order Fractional Derivatives: Algorithms and Applications. Journal of Computational Physics, 293, 184-200.
https://doi.org/10.1016/j.jcp.2014.08.015
[8] Liu, Y., Du, Y., Li, H., Liu, F. and Wang, Y. (2018) Some Second-Order �� Schemes Combined with Finite Element Method for Nonlinear Fractional Cable Equation. Numerical Algorithms, 80, 533-555.
https://doi.org/10.1007/s11075-018-0496-0
[9] Alikhanov, A.A. (2015) A New Difference Scheme for the Time Fractional Diffusion Equation. Journal of Computational Physics, 280, 424-438.
https://doi.org/10.1016/j.jcp.2014.09.031
[10] Myers, T.G., Charpin, J.P.F. and Chapman, S.J. (2002) The Flow and Solidification of a Thin Fluid Film on an Arbitrary Three-Dimensional Surface. Physics of Fluids, 14, 2788-2803.
https://doi.org/10.1063/1.1488599
[11] Myers, T.G. and Charpin, J.P.F. (2004) A Mathematical Model for Atmospheric Ice Accretion and Water Flow on a Cold Surface. International Journal of Heat and Mass Transfer, 47, 5483-5500.
https://doi.org/10.1016/j.ijheatmasstransfer.2004.06.037
[12] Leith, J.R. (2003) Fractal Scaling of Fractional Diffusion Processes. Signal Processing, 83, 2397-2409.
https://doi.org/10.1016/s0165-1684(03)00192-0
[13] Karpman, V.I. (1996) Stabilization of Soliton Instabilities by Higher-Order Dispersion: Fourth-Order Nonlinear Schrödinger-Type Equations. Physical Review E, 53, R1336-R1339.
https://doi.org/10.1103/physreve.53.r1336
[14] Karpman, V.I. and Shagalov, A.G. (2000) Stability of Solitons Described by Nonlinear Schrödinger-Type Equations with Higher-Order Dispersion. Physica D: Nonlinear Phenomena, 144, 194-210.
https://doi.org/10.1016/s0167-2789(00)00078-6
[15] Du, Q., Ju, L. and Tian, L. (2008) Analysis of a Mixed Finite-Volume Discretization of Fourth-Order Equations on General Surfaces. IMA Journal of Numerical Analysis, 29, 376-403.
https://doi.org/10.1093/imanum/drn021
[16] Liu, Y., Fang, Z., Li, H. and He, S. (2014) A Mixed Finite Element Method for a Time-Fractional Fourth-Order Partial Differential Equation. Applied Mathematics and Computation, 243, 703-717.
https://doi.org/10.1016/j.amc.2014.06.023
[17] Liu, Y., Du, Y., Li, H., He, S. and Gao, W. (2015) Finite Difference/Finite Element Method for a Nonlinear Time-Fractional Fourth-Order Reaction-Diffusion Problem. Computers & Mathematics with Applications, 70, 573-591.
https://doi.org/10.1016/j.camwa.2015.05.015
[18] Zhang, T. and Guo, Q. (2018) The Finite Difference/Finite Volume Method for Solving the Fractional Diffusion Equation. Journal of Computational Physics, 375, 120-134.
https://doi.org/10.1016/j.jcp.2018.08.033
[19] Huang, C. and Stynes, M. (2020) Α-Robust Error Analysis of a Mixed Finite Element Method for a Time-Fractional Biharmonic Equation. Numerical Algorithms, 87, 1749-1766.
https://doi.org/10.1007/s11075-020-01036-y
[20] Arenas, A.J., González-Parra, G. and Chen-Charpentier, B.M. (2016) Construction of Nonstandard Finite Difference Schemes for the SI and SIR Epidemic Models of Fractional Order. Mathematics and Computers in Simulation, 121, 48-63.
https://doi.org/10.1016/j.matcom.2015.09.001
[21] Boukal, Y., Zasadzinski, M., Darouach, M. and Radhy, N.E. (2016). Robust Functional Observer Design for Uncertain Fractional-Order Time-Varying Delay Systems. 2016 American Control Conference (ACC), Boston, 6-8 July 2016, 2741-2746.
https://doi.org/10.1109/acc.2016.7525333
[22] Nandal, S. and Pandey, D. (2019) Numerical Solution of Time Fractional Non-Linear Neutral Delay Differential Equations of Fourth-Order. Malaya Journal of Matematik, 7, 579-589.
https://doi.org/10.26637/mjm0703/0035
[23] Nandal, S. and Narain Pandey, D. (2020) Numerical Solution of Non-Linear Fourth Order Fractional Sub-Diffusion Wave Equation with Time Delay. Applied Mathematics and Computation, 369, Article ID: 124900.
https://doi.org/10.1016/j.amc.2019.124900
[24] Hendy, A.S. and De Staelen, R.H. (2020) Theoretical Analysis (Convergence and Stability) of a Difference Approximation for Multiterm Time Fractional Convection Diffusion-Wave Equations with Delay. Mathematics, 8, Article No. 1696.
https://doi.org/10.3390/math8101696
[25] Wang, H. and Yang, Q. (2022) Compact Difference Method for Time-Fractional Neutral Delay Nonlinear Fourth-Order Equation. Engineering, 14, 544-566.
https://doi.org/10.4236/eng.2022.1412041
[26] Xie, H. and Yang, Q. (2022) Compact Difference Scheme for Time-Fractional Nonlinear Fourth-Order Diffusion Equation with Time Delay. Results in Applied Mathematics, 16, Article ID: 100339.
https://doi.org/10.1016/j.rinam.2022.100339
[27] Sun, Z. and Wu, X. (2006) A Fully Discrete Difference Scheme for a Diffusion-Wave System. Applied Numerical Mathematics, 56, 193-209.
https://doi.org/10.1016/j.apnum.2005.03.003
[28] Li, R., Chen, Z. and Wu, W. (2000) Generalized Difference Methods for Differential Equations: Numerical Analysis of Finite Volume Methods. CRC Press.
[29] Zhang, Q., Ran, M. and Xu, D. (2016) Analysis of the Compact Difference Scheme for the Semilinear Fractional Partial Differential Equation with Time Delay. Applicable Analysis, 96, 1867-1884.
https://doi.org/10.1080/00036811.2016.1197914

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.