Numerical treatment of the generalized time-fractional Huxley-Burgers’ equation and its stability examination

1 Introduction

Fractional calculus is a generalization of classical calculus that is concerned with operations of differentiation and integration to fractional order. The concept of the fractional operator was first raised by Marquis de L'Hôpital and G. W. Leibniz in the year 1695. From now on there have been several fundamental works on the fractional derivative and fractional differential equations, and many books were performed in this field, in which we can refer to the books of Ross and Miller [1], Samko et al. [2], and Podlubny [3]. Several numerical methods have been developed to obtain an approximate solution of fractional differential equations. These methods include finite difference method [4,5], homotopy perturbation method [6,7,8], variational iteration method [9,10], Adomian decomposition method [11,12], and spectral methods [13,14,15,16]. The generalized Huxley-Burgers’ equation is in the following form [17]:

where β , ν , η , γ , and δ are parameters, δ > 0 , η ≥ 0 , and γ ∈ ( 0 , 1 ) . Equation (1) is an extended form of the famous Huxley and Burgers’ equations. If we take β = 0 and δ = 1 , equation (1) is reduced to the Huxley equation which describes nerve pulse propagation in nerve fibers and wall motion in liquid crystals [17]:

If we take η = 0 , β = 1 , and δ = 1 , equation (1) is reduced to Burgers’ equation which describes the far field of wave propagation in nonlinear dissipative systems [18]:

It is known that the nonlinear diffusion equations (2) and (3) play important roles in physics. They are of special significance for studying nonlinear phenomena. If we take η ≠ 0 , β ≠ 0 , and δ = 1 , equation (1) becomes the following Huxley and Burgers’ equation:

Equation (1) shows a prototype model for describing the interaction among reaction mechanisms, convection effects, and diffusion transport. This equation was investigated by Satsuma [19] in 1987. The generalized time-fractional Huxley-Burgers’ equation whose solutions will be approximated is of the form [20]:

and is subject to the conditions:

and

We can obtain approximate solutions of the generalized time-fractional Huxley-Burgers’ equation that have applications in various fields of science and engineering. The time-fractional Huxley-Burgers’ equation illustrates intercommunications between reaction mechanisms and diffusion transports and an evolutionary model depicting nerve pulse propagation [21]. Many scholars have recently worked on the fractional Burger-Huxley equation due to the necessity of finding a solution to this equation. In [22], the space-time-fractional Burgers-Huxley problem was solved using a Legendre spectral finite difference method. The author in [23] looked into solitary wave solutions to this problem. In [24], the problem is solved using the finite difference method. In [20], the residual power series method is used to obtain an analytical solution to this problem.

In this work, we have proposed a collocation method for solving the time-fractional Burgers-Huxley equation using the mean value theorem for integrals and cubic B-spline basis functions. In this method, we use cubic B-splines for spatial variables and their derivatives which produce a system of fractional ordinary differential equations. The fractional derivative is simplified by using the mean value theorem for integrals and the finite difference method. B-splines have two useful qualities. One distinguishing characteristic is that the continuity conditions are inherent. As a result, when compared to other piecewise polynomial interpolating functions, the B-spline is the smoothest. B-splines also have small local support, which means that each B-spline function is only non-zero over a few mesh subintervals, resulting in a tightly banded resultant matrix for the discretization equation. B-splines have special advantages due to their smoothness and capacity to handle local occurrences. In combination with the collocation method, this significantly simplifies the solution procedure of differential equations. There is a great reduction of the numerical effort because there is no need to calculate integrals (as in variational methods) to build the final set of algebraic equations that replaces the provided nonlinear differential equation. The current method does not require extra effort to deal with the nonlinear parts, unlike some previous techniques that used various transformations to reduce the problem to a simpler equation. As a result, the equations are readily and elegantly solved utilizing the current method. The Caputo fractional derivative of order α which is used in this paper is defined by,

where Г ( ⋅ ) is the Gamma function.

This paper is organized as follows. In Section 2, a proposed method depends on the cubic B-spline collocation method and the mean value theorem is derived. In Section 3, we obtain the initial vector that is required to start the iterative process of the proposed method. In Section 4, the stability is theoretically discussed by using the von Neumann technique. In Section 5, we illustrate three numerical examples that are introduced to illustrate the efficiency and accuracy of the proposed method. Finally, the conclusion is given in Section 6.

2 Derivation of the method

To approximate u ( x , t ) by cubic B-spline collocation method over finite elements, let the region R = [ a , b ] × [ 0 , T ] be discretized by a set of points R i j which are the vertices of a grid of points ( x i , t j ) where x i = a + i h , h = Δ x for i = 0 , 1 , … , N , and t j = j k , k = Δ t for j = 0 , 1 , … , M . Let ∅ i ( x ) be cubic B-splines with nodes at x N − 2 < x N − 1 < … < x N + 1 < x N + 2 . The cubic B-spline functions ∅ i ( x ) at these knots are given as follows:

The collocation method for approximately solving equation (5) consists in seeking approximations U ( x ,  t ) of u ( x , t ) from the finite dimensional subspace of C 2 [ a ,  b ] which is spanned by the linearly independent set of cubic B-spline functions { ∅ − 1 ( x ) , ∅ 0 ( x ) , … , ∅ N ( x ) , ∅ N + 1 ( x ) } having the form [25]:

where w i ( t ) are the unknown time-dependent quantities. The quantities w i ( t ) are chosen such that U ( x ,  t ) satisfies equation (5) and its conditions (6) and (7) at the points ( x i ,  t j ) for i = 0 , 1 , … ,  N , j = 0 , 1 , … , M . Then U ( x ,  t ) satisfies the following equations:

Substitute equation (9) and the values in Table 1 into equation (10) to obtain

where μ i j = ( w i − 1 j + 4 w i j + w i + 1 j ) δ and ρ i j = ( 1 − ( w i − 1 j + 4 w i j + w i + 1 j ) δ ) ( ( w i − 1 j + 4 w i j + w i + 1 j ) δ − γ ) .

Using the Caputo fractional derivative definition, the time-fractional derivative can be discretized as

By using a piecewise technique, equation (14) becomes

Since ( t j − τ ) − α does not change the sign on the interval [ ( m − 1 ) k ,  m k ] , the mean value theorem for integrals [26] will be applied as follows:

which implies that

where τ j , m α = ( j − m + 1 ) 1 − α − ( j − m ) 1 − α . This enables us to write the Caputo fractional derivative in the following form:

then

with a local truncation error

for 0 < θ < 1 , t m − 1 < ξ m − < τ m ∗ , and τ m ∗ < ξ m + < t m .

Substituting equation (16) in equation (13), to obtain the following system:

From equation (18) we obtain the following:

For simplicity, equation (19) can be rewritten as a recurrence relationship:

where A i = σ − 3 β h μ i j − 6 ν h 2 − η ρ i j , B i = 4 σ + 12 ν h 2 − 4 η ρ i j , C i = σ + 3 β h μ i j − 6 ν h 2 − η ρ i j , and D i j − 1 = σ ∑ m = 1 j τ j , m α ( w i − 1 m − 1 + 4 w i m − 1 + w i + 1 m − 1 ) − σ ∑ m = 1 j − 1 τ j , m α ( w i − 1 m + 4 w i m + w i + 1 m ) + f ( x i , t j ) .

System (20) has to be complemented by the boundary conditions

For each i = 0 , 1 , … ,  N and j = 0 , 1 , … ,  M . The cubic B-spline collocation method produces the ( N + 3 ) × ( N + 3 ) matrix system, when we add the boundary conditions. From equations (20) and (21), we obtain the following system of linear equation:

where

Q = 1 4 1 0 ⋯ 0 0 0 A 0 B 0 C 0 0 ⋯ 0 0 0 0 A 1 B 1 C 1 ⋯ 0 0 0 ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ 0 0 0 0 ⋯ A N B N C N 0 0 0 0 ⋯ 1 4 1 , w j = w − 1 j w 0 j w 1 j ⋮ w N j w N + 1 j , and Q ∗ = φ 1 ( t j ) D 0 j − 1 D 1 j − 1 ⋮ D N j − 1 φ 2 ( t j ) .

3 Initial vector

The initial condition helps us to obtain the initial vector ( w − 1 0 , w 0 0 , w 1 0 , … ,  w N 0 ,  w N + 1 0 ) , which is required to start the iterative process of the proposed method

Using the values in Table 1 into equation (23) gives

The last system (24) consists of ( N + 1 ) equations and ( N + 3 ) unknowns. To solve the system of linear equation (24) we need two additional equations. These equations are obtained as follows:

By using Table 1, in equation (25) we obtain

Combining equations (24) and (26) we obtain the following system of linear equation:

where

which enables us to obtain the initial vector ( w − 1 0 , w 0 0 , w 1 0 , … , w N 0 ,  w N + 1 0 ) .

4 Stability analysis

The von Neumann technique is suitable for linear problems with constant coefficients. To study the stability of our numerical scheme by using the von Neumann technique, we must linearize the nonlinear terms u δ ( x ,  t ) u ′ ( x ,  t ) and ( 1 − u δ ( x ,  t ) ) ( u δ ( x ,  t ) − γ ) u ( x ,  t ) of the time-fractional Huxley and Burgers equation (5), by assuming that the corresponding quantities ρ i − 1 j ,   ρ i j , and ρ i + 1 j are equal to a local constant χ ∗ and μ i − 1 j , μ i j , and μ i + 1 j to λ ∗ in equation (20). According to the von Neumann technique, we assume a solution of the form [27,28]:

where q 2 = − 1 , h is the element size, φ is the mode number, and ξ j is the amplification factor at time level j . Inserting the latter expression for w i j into the system (20) we have the following characteristic equation:

where ϑ i j − 1 = K ( ξ j − 1 e φ q ( i − 1 ) h + 4 ξ j − 1 e φ q i h + ξ j − 1 e φ q ( i + 1 ) h ) + K ∑ m = 1 j − 1 τ j , m α [ ( ξ m − 1 e φ q ( i − 1 ) h + 4 ξ m − 1 e φ q i h + ξ m − 1 e φ q ( i + 1 ) h ) − ( ξ m e φ q ( i − 1 ) h + 4 ξ m e φ q i h + ξ m e φ q ( i + 1 ) h ) ] and K = h 2 6 σ .

After dividing by e φ q i h and using Euler’s formula e φ h q = cos φ h + q sin φ h , the last equation can be simplified as

Equation (30) has been rewritten in a simple form as

where Ω = E H + q Z , E = 2 K ( cos ϕ + 2 ) , H = E + 2 ν ( 1 − cos ϕ ) − h 2 3 η χ ∗ ( cos ϕ + 2 ) , Z = h β λ ∗ sin ϕ , and ϕ = φ h .

For j = 1 in equation (31), we obtain

thus we obtain the following:

where ∣ Ω ∣ = E 2 H 2 + Z 2 . The quantity ( cos ϕ + 2 ) is positive and ( 1 − cos ϕ ) is equal to zero or positive. If we choose h small enough to make h 2 3 η χ ⁎ ( cos ϕ + 2 ) → 0 , then we have

If we put j = 2 in equation (31), we obtain

where ( 1 − τ 2 , 1 α ) , and τ 2 , 1 α are positive quantities. Using inequality (34) implies

then

For j = 3 in equation (31), we obtain

then