Finite Volume Element Method for Fractional Order Neutral Time-Delay Differential Equations ()
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
formula to approximate the fractional order derivative of time based on the
formula, which had a convergence accuracy of the order
.
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
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
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
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
, 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
(1)
The initial and boundary conditions are:
(2)
where
is the delay quantity,
is the neutral time-delay term,
is nonlinear source term with time delay,
,
,
and
are all sufficiently smooth functions.
The Caputo time fractional order derivative is defined as follows
(3)
Assume that the partial derivatives
are continuous in a
neighborhood of
, where
is a positive constant. Denote
(4)
(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
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
be a quasi-uniformly regular triangulation of
, and the triangular element in
is denoted by
.
denotes the set of all nodes in
, and
is the set of all inner points. Let
be the diameter of the element
,
.
For any
, take the adjacent nodes of
to be
, and the midpoint of
to be
. Let the barycenter of the triangular element
be
and connect
to get the polygonal region
surrounding
. Similarly, for any
, we can get the corresponding semi-enclosed polygonal region
. The polygonal region
is called a dual element, and all the dual elements form the barycenter dual subdivision of
, denoted by
.
2.1. Trial Function Space and Test Function Space
The trial function space
is chose as the standard linear element space on
, then
, where
is the set of all linear polynomials on the triangle element
.
The test function space
is chosed as the piecewise constant space on
, then
. Thus,
can be viewed as a space spanned by the following basis functions: for
,
(6)
2.2. Finite Volume Element Format
Firstly, the intermediate function
is introduced, then the original equation can be rewritten as
(7)
Multiplying
and
at both sides of the first equation and the second equation of (7), respectively and utilizing Green’s formula to get:
(8)
where the bilinear form
.
Choosing the trial function space
as described above, and we can obtain the semi-discrete finite volume element format: find
such that
(9)
where
and
are the elliptic projections of
and
on
, respectively.
Next, we discrete time fractional order derivatives and deal with nonlinear terms. Let
be a positive integer, and make
,
,
. For convenience, assume that
is a positive integer.
Lemma 1 [27] Let
, then
(I)
(10)
(II)
(11)
Lemma 2 [27] For
order Caputo fractional order derivatives, there are
(12)
where
.
When
, for the right nonlinear term, Taylor’s formula yields:
(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
such that
(14)
where
.
For the purpose of subsequent research, the following definitions and lemmas are presented.
Definition 1. [28] The interpolation operator
is defined as follows:
(15)
Definition 2. [28] Introduce the elliptic projection operator
, defined by the following format:
(16)
Lemma 3. [28] Let
, and
is the elliptic projection of
, then
(17)
Lemma 4. [18] [28] Let
, we have
(I)
(18)
(II)
(19)
(III)
(20)
Lemma 5. [18] [28]
and
are equivalent on
, i.e., there exist constants
such that
(21)
Lemma 6. [28] There exist positive constants
and
such that the following equations hold
(22)
(23)
(24)
Lemma 7. [29] Let
be two non-negative sequences. If
(25)
then
(26)
where
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
(27)
where
.
To prove that (14) is uniquely solvable, it is sufficient to prove that its corresponding homogeneous equations (27) has only zero solutions. Let
in (27), then
(28)
(29)
According to Lemma 4
(30)
Substituting (30) into (29) gives
(31)
Substituting (31) into (28) gives
(32)
According to Lemma 6
(33)
We get
(34)
So
. 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
and
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
such that
(35)
(36)
Proof. Let
(37)
(38)
then
(39)
(40)
Let
in (8) and subtract from (14). According to the elliptic projection in Definition 2, we get
(41)
(42)
where
and
is a positive constant that does not depend on
.
Let
and
in (41) and (42), respectively, and use Lemma 6 to get
(43)
(44)
Let
in (42), and use Lemma 6 to get
(45)
Substitute (44) and (45) into (43) to obtain
(46)
From the definition of
, multiplying
on both sides of the (46) to get
(47)
In order to estimate the order of convergence of the solution
, we analyze the terms at the right side of the equation (47). According to the properties of the function
, it follows that
(48)
where
.
According to the Cauchy-Schwarz inequality and (48), we get
(49)
For the remaining terms at the right side, using the Cauchy-Schwartz inequality and substituting (49) into (47) to obtain
(50)
Then, we can get
(51)
According to the definition of
, and take a sufficiently small
to get
(52)
(52) will be used in the subsequent proof to estimate the convergence order of the solution
.
In order to estimate the order of convergence of the solution
, the terms at the right side of (47) are analyzed. Similar to the derivation of (52), we get
(53)
(53) will be used in the subsequent proof to estimate the convergence order of the solution
.
Let
. If
is an integer, then the time period
is divided into
segments, where the length of each segment is
. If
is not an integer, then the time period
is divided into
segments, where the length of the first
segments are
and the length of the last segment is
. Since
, then
on the first segment,
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
. The function is known at
and
, so
,
on the first segment. Substitute
,
into (52) to get
(54)
By Lemma 1 and Lemma 7, it follows that
(55)
According to the equivalence of norms in Lemma 5 and the property of elliptic projection in Lemma 3,
(56)
(57)
For
,
(58)
According to (58) and Lemma 3, it follows that
(59)
Substituting (56), (57) and (59) into (55), we get
(60)
Let
, and then
(61)
Using the equivalence of norms in Lemma 5 and the definition of
to get
(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
. Substituting
,
into (53) yields
(63)
The right side terms of (63) are analyzed below. From (56), (57) and (59), we can obtain
(64)
From (61), we can obtain
(65)
(66)
Substituting (64)-(66) into (63) and dividing both sides by
to get
(67)
Using the equivalence of the norms in Lemma 5, we can obtain
(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
(69)
According to lemma 7,
(70)
From (67),
(71)
Let
. According to the properties of elliptic projection and (71), we can get
(72)
Using the equivalence of the norms in Lemma 5 and the definition of
(73)
Finally, by using the triangle inequality, it can be proved that (35) holds on the second segment.
Next, from (53) we can get
(74)
The right side terms of (74) are analyzed below. From (61) and (67), we get
(75)
(76)
According to the properties of elliptical projection and (76), it follows that
(77)
From (72) and (75), we have
(78)
From (72), we get
(79)
Substituting (77), (78), (79) into (74) and dividing both sides by
, we get
(80)
Using the equivalence of the norms in Lemma 5, we can obtain
(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
be the solution of the following equations
(82)
where
are the perturbation terms of the initial function.
Let
.
Theorem 3 There exists a constant C such that
(83)
(84)
Proof. Subtracting (82) from (14), we get
(85)
(86)
According to the derivation of (52) and (53) in the convergence analysis, and similarly for (85) and (86), we can get
(87)
(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
,
, and define
(89)
The convergence order is calculated as follows:
(90)
Example 1. Consider the exact solution of problem (1) as:
(91)
The nonlinear source term with time delay is
(92)
where
(93)
Taking
, the above problem can be solved by using the finite volume element format constructed in this paper. When
, we fix the time step
, and change the spatial step to get the errors and the spatial convergence order of
. We fix the spatial step
, and change the time step to get the errors and the time convergence order of
. 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
-normand
-norm are approximately equal to 2, and the temporal convergence order in
-normand
-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
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
2.087e−01 |
− |
9.456e−02 |
− |
2.483 |
− |
1.102 |
− |
|
9.305e−02 |
1.995 |
4.204e−02 |
1.994 |
1.102 |
2.001 |
4.882e−01 |
2.004 |
|
5.224e−02 |
1.998 |
2.357e−02 |
2.004 |
6.180e−01 |
2.006 |
2.728e−01 |
2.014 |
|
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
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
1.338e−01 |
− |
6.550e−02 |
− |
2.759 |
− |
1.271 |
− |
|
4.157e−02 |
1.686 |
2.011e−02 |
1.703 |
8.244e−01 |
1.746 |
3.743e−01 |
1.764 |
|
1.363e−02 |
1.608 |
6.393e−03 |
1.653 |
2.356 e−01 |
1.803 |
1.032e−01 |
1.858 |
When
, we fix the time step
, and change the spatial step to get the errors and spatial convergence order of
. We fix the spatial step
, and change the time step to get the errors and temporal convergence order of
. The results are given by Table 3 and Table 4, respectively.
Table 3. Errors and order of spatial convergence of finite volume element solutions
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
2.449e−01 |
− |
1.110e−01 |
− |
2.964 |
− |
1.309 |
− |
|
1.103e−01 |
1.995 |
4.992e−02 |
1.994 |
1.341 |
2.001 |
5.902e−01 |
2.004 |
|
6.286e−02 |
1.974 |
2.842e−02 |
1.972 |
7.705e−01 |
1.965 |
3.378e−01 |
1.970 |
|
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
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
1.985e−01 |
− |
9.699e−02 |
− |
4.063 |
− |
1.870 |
− |
|
6.604e−02 |
1.587 |
3.182e−02 |
1.607 |
1.280 |
1.665 |
5.815e−01 |
1.685 |
|
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
-normand
-norm are approximately equal to 2, and the temporal convergence order in
-normand
-norm are close to 1.5. This is consistent with the theoretical analysis.
When
, we fix the time step
, and change the spatial step to get the errors and spatial convergence order of
. We fix the spatial step
, and change the time step to get the errors and temporal convergence order of
. The results are given by Table 5 and Table 6, respectively.
Table 5. Errors and spatial convergence order of finite volume element solutions
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
2.825e−01 |
− |
1.280e−01 |
− |
3.431 |
− |
1.514 |
− |
|
1.266e−01 |
1.984 |
5.728e−02 |
1.982 |
1.538 |
1.983 |
6.767e−01 |
1.986 |
|
7.162e−02 |
1.979 |
3.237e−02 |
1.984 |
8.728e−01 |
1.974 |
3.830e−01 |
1.982 |
|
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
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
2.912e−01 |
− |
1.420e−01 |
− |
5.917 |
− |
2.718 |
− |
|
1.028e−01 |
1.501 |
4.938e−02 |
1.524 |
1.958 |
1.595 |
8.858e−01 |
1.617 |
|
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
-normand
-norm are approximately equal to 2, and the temporal convergence order in
-normand
-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
. (a) The finite volume element solution
; (b) The exact solution
.
Figure 2. Comparison of
. (a) The finite volume element solution
; (b) The exact solution
.
(94)
The nonlinear source term with time delay is
(95)
where
(96)
Taking
, and we choose the numerical solution when
,
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
when
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
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
3.207e−01 |
− |
1.598e−02 |
− |
1.102 |
− |
5.393e−01 |
− |
|
8.030e−02 |
1.998 |
4.029e−02 |
1.9878 |
2.753e−01 |
2.001 |
1.356e−01 |
1.991 |
|
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
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
3.756e−02 |
− |
1.847e−02 |
− |
7.886e−01 |
− |
3.653e−01 |
− |
|
9.360e−03 |
2.004 |
4.604e−03 |
2.003 |
1.965e−01 |
2.004 |
9.105e−02 |
2.004 |
|
2.354e−03 |
1.991 |
1.158e−03 |
1.990 |
4.946 e−02 |
1.990 |
2.290e−02 |
1.990 |
order in
-normand
-norm are approximately equal to 2, and the time convergence order of
-normand
-norm are close to 2.
When
, the errors, temporal and spatial convergence order of
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
-normand
-norm are approximately equal to 2, and the time convergence order in
-normand
-norm are close to 2.
When
, the errors, temporal and spatial convergence order of
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
-normand
-norm are approximately equal to 2, and the time convergence order in
-normand
-norm are close to 2.
Table 9. Errors and order of spatial convergence of finite volume element solutions
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
3.206e−01 |
− |
1.597e−01 |
− |
1.099 |
− |
5.377e−01 |
− |
|
8.025e−02 |
1.998 |
4.026e−02 |
1.987 |
2.745e−01 |
2.001 |
1.352e−01 |
1.991 |
|
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
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
3.702e−02 |
− |
1.820e−02 |
− |
7.779e−01 |
− |
3.601e−01 |
− |
|
9.149e−03 |
2.016 |
4.499e−03 |
2.017 |
1.923e−01 |
2.015 |
8.897e−02 |
2.016 |
|
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
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
3.204e−01 |
− |
1.596e−01 |
− |
1.095 |
− |
5.362e−01 |
− |
|
8.021e−02 |
1.998 |
4.024e−02 |
1.987 |
2.737e−01 |
2.001 |
1.348e−01 |
1.991 |
|
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
,
when
.
|
|
order |
|
order |
|
order |
|
order |
|
3.601e−02 |
− |
1.770e−02 |
− |
7.579e−01 |
− |
3.501e−01 |
− |
|
8.822e−03 |
2.029 |
4.335e−03 |
2.030 |
1.859e−01 |
2.026 |
8.574e−02 |
2.030 |
|
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.