Research

# A generalized groundwater flow equation using the concept of variable-order derivative

Abdon Atangana* and Joseph Francois Botha

Author Affiliations

Institute for Groundwater Studies, Faculty of Natural and Agricultural Sciences, University of the Free State, P.O. Box 9301, Bloemfontein, South Africa

For all author emails, please log on.

Boundary Value Problems 2013, 2013:53  doi:10.1186/1687-2770-2013-53

 Received: 30 January 2013 Accepted: 23 February 2013 Published: 14 March 2013

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### Abstract

In this paper, the groundwater flow equation is generalized using the concept of the variational order derivative. We present a numerical solution of the modified groundwater flow equation with the variational order derivative. We solve the generalized equation with the Crank-Nicholson technique. Numerical methods typically yield approximate solutions to the governing equation through the discretization of space and time and can relax the rigid idealized conditions of analytical models or lumped-parameter models. They can therefore be more realistic and flexible for simulating field conditions. Within the discredited problem domain, the variable internal properties, boundaries, and stresses of the system are approximated. We perform the stability and convergence analysis of the Crank-Nicholson method and complete the paper with some illustrative computational examples and their simulations.

##### Keywords:
groundwater flow equation; variable order derivative; Crank-Nicholson scheme; stability; convergence

### 1 Introduction

A problem that arises naturally in groundwater investigations is to choose an appropriate geometry for the geological system in which the flow occurs. This geological formation, through which the groundwater flows, changes in time and space. Figure 1 shows the groundwater flow direction toward the pumping well during the extraction of water from the aquifer. The main effort is to establish a suitable mathematical relation between different observables, or what [1] calls a mathematical model for the observables. An ideal mathematical model should not merely provide links between different observables, but also lead to a better understanding of the phenomenon. The attractiveness of this approach is that the mathematical model could, in principle, be used to investigate the future behavior of a given phenomenon under various conditions. This leads us to the groundwater flow equation. This model is based on the conventional, saturated groundwater flow equation for density-independent flow:

S 0 ( x ̲ , t ) t Φ ( x ̲ , t ) = [ K ̲ ̲ Φ ] + f ( x , t ) . (1.1)

Figure 1. Groundwater flow in the borehole during pumping test.

Here S 0 is the specific storativity, K ̲ ̲ is the hydraulic conductivity tensor of the aquifer, Φ is the piezometric head, f ( x , t ) is the strength of any sources or sinks. As a review of the derivation of Eq. (1.1), we will show the Darcy law

q ( x , t ) = K ̲ ̲ [ Φ ( x , t ) ] (1.2)

is used as a keystone in the derivation of Eq. (1.1). This law, proposed by Darcy early in the nineteenth century, relies on experimental results obtained from the flow of water through a one-dimensional sand column. Alternatively, Darcy’s law states that the rate of flow through a porous medium is proportional to the loss of head and inversely proportional to the length of the flow path. Note that the specific discharge q ( x , t ) has the dimensions of velocity. Recent investigations [2] suggest that the flow is also influenced by the geometry of the bedding parallel fractures, the feature that equation (1.1) cannot account for. It is therefore possible that equation (1.1) may not be applicable to the flow in these fractured aquifers. In an attempt to circumvent this problem, Barker [3] introduced the model in which the geometry of the aquifer is regarded as a fractal. Although this model has been applied with reasonable success in the analysis of hydraulic tests from boreholes in Karoo aquifers [4], it introduces parameters for which no sound definition exists in the case of non-integer flow dimensions. Recently [5,6], the concept of a fractional-order derivative was used to generalize the groundwater flow equation. However, it has been found that the constant-order fractional diffusion equations are not capable of characterizing some complex diffusion processes, for instance, diffusion process in aninhomogeneous or heterogeneous medium [7].

In addition, when we consider the diffusion process in a porous medium, if the medium structure or external field changes with time, in this situation, the constant-order fractional diffusion equation model cannot be used to well characterize such a phenomenon [8-19]. This is the case of the groundwater flow problem, the medium through which the flow occurs is heterogeneous and changes with time. Still in some biology diffusion processes, the concentration of particles will determine the diffusion pattern [10,11]. To solve the above problems, the variable-order (VO) fractional diffusion equation models have been suggested for use [12]. This present work is therefore devoted to the discussion underpinning the description of the groundwater flow equation with the variable-order derivative.

### 2 Modified groundwater flow equation

For the readers that are not acquainted with the concept of the variational order derivative, we start this section by presenting the basic definition of this derivative.

#### 2.1 Variational order differential operator

Let f : R R , x f ( x ) denote a continuous and necessary differentiable, let α ( x ) be a continuous function in ( 0 , 1 ] . Then its variational order differential in [ a , ) is defined as

D 0 α ( x ) ( f ( x ) ) = 1 Γ ( 1 α ( x ) ) 0 x ( x t ) α ( t ) d f ( t ) d t d t . (2.1)

The above derivative is called the Caputo variational order differential operator; in addition, the derivative of the constant is zero.

#### 2.2 Problem formulation

Groundwater models describe the groundwater flow and transport processes using mathematical equations based on certain simplifying assumptions. These assumptions typically involve the direction of flow, geometry of the aquifer, the heterogeneity or anisotropy of sediments or bedrock within the aquifer. This geological formation, through which the groundwater flows, changes in time and space.

The simplest generalization of the groundwater flow equation, which incidentally is also in accord with true physics of the phenomenon, is to assume that water level is not in a steady but transient state. Theis (1935) [13] was the first to develop a formula for unsteady-state flow that introduces the time factor and the storativity. He noted that when a well penetrating an extensive confined aquifer is pumped at a constant rate, the influence of the discharge extends outward with time. The rate of decline of head multiplied by the storativity and summed over the area of influence equals the discharge. The unsteady-state (or Theis) equation, which was derived from the analogy between the flow of groundwater and the conduction of heat, is perhaps the most widely used partial differential equation in groundwater investigations

S D t Φ ( r , t ) = T D r r Φ ( r , t ) + 1 r D r Φ ( r , t ) . (2.2)

The above equation is classified under a parabolic equation. To include explicitly the variability of the medium through which the flow takes place, the standard version of the partial derivative with respect to time is replaced here with variable-order (VO) fractional to obtain

{ S D t α ( x , t ) Φ ( r , t ) = T D r r Φ ( r , t ) + 1 r D r Φ ( r , t ) , 0 < α ( x , t ) 1 . (2.3)

### 3 Numerical solution

Environmental phenomena such as groundwater flow described by variational order derivative are highly complex phenomena, which do not lend themselves readily to the analysis of analytical models. The discussion presented in this section will therefore be devoted to the derivation of a numerical solution to groundwater flow equation (2.3).

Numerical methods yield approximate solutions to the governing equation through the discretization of space and time. Within the discredited problem domain, the variable internal properties, boundaries, and stresses of the system are approximated. Deterministic, distributed-parameter, numerical models can relax the rigid idealized conditions of analytical models or lumped-parameter models, and they can therefore be more realistic and flexible for simulating field conditions. The finite difference schemes for constant-order time or space fractional diffusion equations have been widely studied [14-19]. For constant-order time fractional diffusion equations, Chen et al. proposed an implicit difference approximation scheme [20]. Yuste et al. introduced weighted average finite difference methods [21]. Podlubny et al. proposed the matrix approach for fractional diffusion equations [22], and Hanert proposed a flexible numerical scheme for the discretization of the space-time fractional diffusion equation [23]. Recently, Zhuang et al. considered the numerical schemes for VO space fractional advection-dispersion equation [16], Lin et al. investigated the explicit scheme for VO nonlinear space fractional diffusion equation [24]. Before applying the numerical methods, we assume Eq. (2.3) has a unique and sufficiently smooth solution. To establish the numerical schemes for the above equation, we let x l = l h , 0 l M , M h = L , t k = k τ , 0 k N , N τ = T , h is the step and τ is the time size, M and N are grid points.

#### 3.1 Crank-Nicholson scheme [24]

We introduce the Crank-Nicholson scheme as follows. Firstly, the discretization of first- and second-order space derivative is stated as follows:

(3.1)

(3.2)

(3.3)

The Crank-Nicholson scheme for the VO time fractional groundwater flow model can be stated as follows:

α l k + 1 Φ ( r l , t k + 1 ) t α l k + 1 = τ α l k + 1 Γ ( 2 α l k + 1 ) ( Φ ( r l , t k + 1 ) Φ ( r l , t k ) + j = 1 k [ Φ ( r l , t k + 1 j ) Φ ( r l , t k j ) ] [ ( j + 1 ) 1 α l k + 1 ( j ) 1 α l k + 1 ] ) . (3.4)

Now, replacing equations (3.1), (3.2), (3.3) and (3.4) in (2.3), we obtain the following:

[ S τ α l k + 1 Γ ( 2 α l k + 1 ) ( Φ ( r l , t k + 1 ) Φ ( r l , t k ) + j = 1 k [ Φ ( r l , t k + 1 j ) Φ ( r l , t k j ) ] [ ( j + 1 ) 1 α l k + 1 ( j ) 1 α l k + 1 ] ) ] = T [ 1 2 ( ( Φ ( r l + 1 , t k + 1 ) 2 Φ ( r l , t k + 1 ) + Φ ( r l 1 , t k + 1 ) ( h ) 2 ) + ( Φ ( r l + 1 , t k ) 2 Φ ( r l , t k ) + Φ ( r l 1 , t k ) ( h ) 2 ) ) ] + 1 r l [ 1 2 ( ( Φ ( r l + 1 , t k + 1 ) Φ ( r l 1 , t k + 1 ) 2 ( h ) ) + ( Φ ( r l + 1 , t k ) Φ ( r l 1 , t k ) 2 ( h ) ) ) ] .

For simplicity, let us put:

b j l , k + 1 = ( j + 1 ) 1 α l k + 1 ( j ) 1 α l k + 1 ; T l k + 1 = Γ ( 2 α l k + 1 ) τ α l k + 1 S h 2 T ; G l k + 1 = Γ ( 2 α l k + 1 ) τ α l k + 1 S h and λ j l , k + 1 = b j 1 l , k + 1 b j l , k + 1 . (3.5)

Equation (3.5) becomes

Φ l k + 1 ( 1 + 2 T l k + 1 ) = Φ l + 1 k + 1 ( T l k + 1 + G l k + 1 r l ) + Φ l 1 k + 1 ( T l k + 1 G l k + 1 r l ) + Φ l + 1 k ( T l k + 1 G l k + 1 r l ) + Φ l k + 1 ( 1 + 2 T l k + 1 ) + j = 1 k ( Φ l k + 1 j Φ l k j ) λ j l , k + 1 G l k + 1 . (3.6)

### 4 Stability analysis of the Crank-Nicholson scheme

In this section, we will analyze the stability conditions of the Crank-Nicholson scheme for the generalized groundwater flow equation.

Let ζ l k = Φ l k Θ l k . Here Θ l k is the approximate solution at the point ( x l , t k ) ( k = 1 , 2 , , N , l = 1 , 2 , , M 1 ) and, in addition, ζ k = [ ζ 1 k , ζ 2 k , , ζ M 1 k ] T and the function ζ k ( x ) is chosen to be

ζ k ( x ) = { ζ l k if  x l h 2 < x x l + h 2 , l = 1 , 2 , , M 1 , 0 if  L h 2 < x L . (4.1)

Then the function ζ k ( x ) can be expressed in Fourier series as follows:

ζ k ( x ) = m = m = δ m ( m ) exp [ 2 i π m k / L ] , δ k ( x ) = 1 L 0 L ρ k ( x ) exp [ 2 i π m x L ] d x . (4.2)

It was established by Chen et al.[19] that

ρ 2 2 2 = m = m = δ k ( m ) 2 . (4.3)

Observe that for all k , l 1 , 0 1 α l k + 1 < 1 , in addition, according to the problem in point, the storativity S, and transmissivity T are positive constants. Then the following properties of the coefficients T l k + 1 , G l k + 1 , λ j l , k + 1 , and b l k + 1 can be established:

1 . G l k + 1 , T l k + 1 are positive for all  l = 1 , 2 , , M 1 ; 2 . 0 < λ j l , k λ j 1 l , k 1 for all  l = 1 , 2 , , M 1 ; 3 . 0 b j l , k 1 , j = 0 k 1 b j + 1 l , k + 1 = 1 λ k l , k + 1 for all  l = 1 , 2 , , M 1 . (4.4)

It is customary in groundwater investigations to choose a point on the centerline of the pumped borehole as a reference for the observations, and therefore neither the drawdown nor its derivatives will vanish at the origin, as required. In such situations where the distribution of the piezometric head in the aquifer is a decreasing function of the distance from the borehole, the expression 1 r l 0 . Under this situation, the error committed while approximating the solution of the generalized groundwater flow equation with the Crank-Nicholson scheme can be presented as follows:

ζ l k + 1 ( 1 + 2 T l k + 1 ) = ζ l + 1 k + 1 ( T l k + 1 ) + ζ l 1 k + 1 ( T l k + 1 ) + ζ l + 1 k ( T l k + 1 ) + ζ l k ( 1 + 2 T l k + 1 ) + j = 1 k ( ζ l k + 1 j ζ l k j ) λ j l , k + 1 G l k + 1 . (4.5a)

If we assume that ζ l k in equation (4.1) can be put in the delta-exponential form as follows:

ζ l k = δ k exp [ i ϕ l k ] , (4.5b)

where ϕ is a real spatial wave number, new replacing the above equation (4.5b) into (4.5a), we obtain

(4.6)

Equation (4.6) can be written in the following form:

δ 1 = [ 1 4 T l 1 sin 2 ( ϕ h 2 ) ] δ 0 [ 1 + 4 T l 1 sin 2 ( ϕ h 2 ) ] , δ k + 1 = [ 1 4 T l 1 + k sin 2 ( ϕ h 2 ) e 1 l , k + 1 ] δ k + j = 0 k 1 λ j + 1 l , k + 1 δ k j + λ k l , k + 1 δ 0 [ 1 + 4 T l 1 + k sin 2 ( ϕ h 2 ) ] . (4.7)

Our next concern here is to show that for all k = 1 , 2 , , N 1 , the solution of equation (4.7) satisfies the following condition:

| δ k | < | δ 0 | .

To achieve this, we make use of the recurrence technique on the natural number k.

For k = 1 and remembering that d l k + 1 , b l k + 1 are positive for all l = 1 , 2 , , M 1 , we obtain

| δ 1 | | δ 0 | = | [ 1 4 T l 1 sin 2 ( ϕ h 2 ) ] [ 1 + 4 T l 1 sin 2 ( ϕ h 2 ) ] | < 1 . (4.8)

Assuming that for m = 2 , 3 , , k the property is verified, we get

| δ k + 1 | = | [ 1 4 T l 1 + k sin 2 ( ϕ h 2 ) e 1 l , k + 1 ] δ k + j = 0 k 1 λ j + 1 l , k + 1 δ k j + λ k l , k + 1 δ 0 [ 1 + 4 T l 1 + k sin 2 ( ϕ h 2 ) ] | . (4.9)

Making use of the triangular inequality we obtain

| δ k + 1 | | 1 4 T l 1 + k sin 2 ( ϕ h 2 ) e 1 l , k + 1 | | δ k | + | j = 0 k 1 p j + 1 l , k + 1 δ k j | + | e k l , k + 1 δ 0 | | 1 + 4 T l 1 + k sin 2 ( ϕ h 2 ) | . (4.10)

Using the recurrence hypothesis, we have

(4.11)

which this completes the proof.

### 5 Convergence analysis of the Crank-Nicholson scheme

If we assume that Φ ( r l , t k ) ( l = 1 , 2 , , M , k = 1 , 2 , , N 1 ) is the exact solution of our problem at the point ( r l , t k ) , by letting Ω l k = Φ ( r l , t k ) Φ l k and Ω k = ( 0 , Ω 1 k , Ω 2 k , , Ω M 1 k ) , substituting this in equation (4.7), we obtain

ζ l 1 ( 1 + 2 T l 1 ) ζ l + 1 1 ( T l 1 ) ζ l 1 1 ( T l 1 ) = R l 1 for  k = 0 , ζ l 1 + k ( 1 + 2 T l 1 + k ) ζ l + 1 k + 1 ( T l 1 + k ) ζ l 1 1 + k ( T l 1 + l ) = R l 1 + k + j = 0 k 1 Ω l k j λ j + 1 l , k + 1 for  k 1 . (5.1)

Here,

R l 1 + k = Φ ( r l , t k + 1 ) j = 0 k 1 Φ ( r l , t k j ) λ j + 1 l , k + 1 + b 1 l , k + 1 Φ ( r l , t 0 ) T l 1 + k [ Φ ( r l + 1 , t k + 1 ) 2 Φ ( r l , t k + 1 ) + Φ ( r l 1 , t k + 1 ) ] . (5.2)

From equations (3.1) and (3.4), we have

2 Φ ( r l , t k + 1 ) r 2 + h 2 V 1 = 1 2 ( ( Φ ( r l + 1 , t k + 1 ) 2 Φ ( r l , t k + 1 ) + Φ ( r l 1 , t k + 1 ) ) h 2 2 Φ ( r l , t k + 1 ) r 2 + h 2 V 1 = + ( Φ ( r l + 1 , t k ) 2 Φ ( r l , t k ) + Φ ( r l 1 , t k ) ) h 2 ) , α l k + 1 Φ ( r l , t k + 1 ) t α l k + 1 + τ V 2 = τ α l k + 1 Γ ( 2 α l k + 1 ) ( Φ ( r l , t k + 1 ) Φ ( r l , t k ) + j = 1 k [ Φ ( r l , t k + 1 j ) Φ ( r l , t k j ) ] λ j l , k ) .

From the above, we have that

R l k + 1 K ( τ 1 + α l k + 1 + h 2 τ α l k ) , (5.3)

where K 1 , K 2 , and K are constants. Taking into account the Caputo-type fractional derivative, the detailed error analysis on the above schemes can refer to the work by Diethelm et al.[25] and further work by Li and Tao [26].

Lemma 1 Ω k + 1 K ( τ 1 + α l k + 1 + h 2 τ α l k ) ( Ω j l , k + 1 ) 1 is true for ( k = 0 , 1 , 2 , , N 1 ), where w k = max 1 l M 1 ( Ω k ) , Kis a constant. In addition,

α k + 1 = { min 1 l M 1 α l k + 1 , if   τ < 1 , max 1 l M 1 α l k + 1 , if   τ > 1 .

This can be achieved via the recurrence technique on the natural numberk.

When k = 0 , we have the following:

| Ω l 1 | ( c l 1 + b l 1 ) | w l + 1 1 | + ( c l 1 2 b l 1 ) | w l 1 1 | = | F l 1 | V ( τ 1 + α l k + 1 + h 2 τ α l k ) ( λ j l , k + 1 ) 1 . (5.4)

Now suppose that Ω i + 1 K ( τ 1 + α l i + 1 + h 2 τ α l i ) ( λ j l , i + 1 ) 1 , i = 1 , , N 2 . Then

(5.5)

which completes the proof.

Theorem 1The Crank-Nicholson scheme is convergent, and there exists a positive constantVsuch that

| Φ l k Φ ( x l , t k ) | K ( τ + h 2 ) , l = 1 , 2 , , M 1 , k = 1 , 2 , , N . (5.6)

An interested reader can find the solvability of the Crank-Nicholson scheme in the work done by [24]. Therefore, the details of the proof will not be presented in this paper.

### 6 Numerical results

An image is worth ten thousand words; therefore, we devote this section to the numerical simulations of the solution of the generalized groundwater flow equation. The parameters used in the simulation are given as S = 0.0086 , T = 400 , r [ 0.5 , 35 ] , and t [ 0 , 10 ] .

Example

{ S D t α ( x , t ) Φ ( r , t ) = T D r r Φ ( r , t ) + 1 r D r Φ ( r , t ) , α ( x , t ) = 1 sin ( r t ) , Φ ( r , 0 ) = 2 . (6.1)

The numerical simulations of the approximate solution are presented in figures. First, Figure 2 shows the solution as a function of time and space. Figure 3 shows the solution as a function of time for a fixed distance. Figure 4 shows the solution as a function of space for fixed time.

Figure 2. Drawdown as a function of space and time.

Figure 3. Drawdown as a function of time for a fixed distance.

Figure 4. Drawdown as a function of space for fixed time.

### 7 Conclusion

In this paper, the groundwater flow equation was generalized using the concept of a variational order derivative. The new equation was solved numerically via the Crank-Nicholson technique. We presented in detail the stability and the convergence of this problem. We presented numerical simulations.

### Competing interests

The authors declare that they have no competing interests.

### Authors’ contributions

AA wrote the first draft and FJB corrected the final version. Both authors read and approved the final draft.

### Acknowledgements

Dedicated to Professor Hari M Srivastava.

The authors would like to thank the referee for some valuable comments and helpful suggestions. This study was supported by the National Research Fund of South Africa.

### References

1. Botha, JF, Buys, J, Verwey, JP, Tredoux, G, Moodie, JW, Hodgkiss, M: Modelling groundwater contamination in the Atlantis aquifer. WRC Report No 175/1/90. Water Research Commission, P.O. Box 824, Pretoria 0001 (2004)

2. Van Tonder, GJ, Botha, JF, Chiang, WH, Kunstmann, H, Xu, Y: Estimation of sustainable yields of boreholes in fractured rock formations. J. Hydrol.. 241, 70–90 (2001). Publisher Full Text

3. Barker, JA: A generalised radial flow model for hydraulic tests in fractured rock. Water Resour. Res.. 24(10), 1796–1804 (1988). Publisher Full Text

4. Van Der Voort, I: Risk-based decision tools for managing and protecting groundwater resources. PhD dissertation, University of the Free State, P.O. Box 339, Bloemfontein (2001)

5. Atangana, A: Numerical solution of space-time fractional order derivative of groundwater flow equation. Istanbul, June 20-24. (2012)

6. Botha, JF, Cloot, AH: A generalized groundwater flow equation using the concept of non-integer order. Water SA. 32(1), 1–7 (2006)

7. Solomon, TH, Weeks, ER, Swinney, HL: Observation of anomalous diffusion and Lévy flights in a two-dimensional rotating flow. Phys. Rev. Lett.. 71, 3975–3978 (1993). PubMed Abstract | Publisher Full Text

8. Magin, RL, Abdullah, O, Baleanu, D, Zhou, XJ: Anomalous diffusion expressed through fractional order differential operators in the Bloch-Torrey equation. J. Magn. Reson.. 190, 255–270 (2008). PubMed Abstract | Publisher Full Text

9. Sun, HG, Chen, W, Chen, YQ: Variable order fractional differential operators in anomalous diffusion modeling. Phys. A. 388, 4586–4592 (2009). Publisher Full Text

10. Umarov, S, Steinberg, S: Variable order differential equations and diffusion with changing modes. Z. Anal. Anwend.. 28, 431–450 (2009)

11. Ross, B, Samko, SG: Fractional integration operator of variable order in the Hölder space H λ ( x ) . Int. J. Math. Math. Sci.. 18, 777–788 (1995). Publisher Full Text

12. Chen, YQ, Moore, KL: Discretization schemes for fractional-order differentiators and integrators. IEEE Trans. Circ. Syst. I: Fund. Theoret. Appl.. 49, 363–367 (2002). Publisher Full Text

13. Theis, CV: The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage. Trans. Am. Geophys. Union. 16, 519–524 (1935). Publisher Full Text

14. Zhang, Y: A finite difference method for fractional partial differential equation. Appl. Math. Comput.. 215, 524–529 (2009). Publisher Full Text

15. Tadjeran, C, Meerschaert, MM, Scheffler, HP: A second order accurate numerical approximation for the fractional diffusion equation. J. Comput. Phys.. 213, 205–213 (2006). Publisher Full Text

16. Meerschaert, MM, Tadjeran, C: Finite difference approximations for fractional advection dispersion equations. J. Comput. Appl. Math.. 172, 65–77 (2004). PubMed Abstract | Publisher Full Text | PubMed Central Full Text

17. Deng, WH: Numerical algorithm for the time fractional Fokker-Planck equation. J. Comput. Phys.. 227, 1510–1522 (2007). Publisher Full Text

18. Li, CP, Chen, A, Ye, JJ: Numerical approaches to fractional calculus and fractional ordinary differential equation. J. Comput. Phys.. 230, 3352–3368 (2011). PubMed Abstract | Publisher Full Text | PubMed Central Full Text

19. Chen, CM, Liu, F, Turner, I, Anh, V: A Fourier method for the fractional diffusion equation describing sub-diffusion. J. Comput. Phys.. 227, 886–897 (2007). PubMed Abstract | Publisher Full Text | PubMed Central Full Text

20. Yuste, SB, Acedo, L: An explicit finite difference method and a new Von Neumann-type stability analysis for fractional diffusion equations. SIAM J. Numer. Anal.. 42, 1862–1874 (2005). Publisher Full Text

21. Podlubny, I, Chechkin, A, Skovranek, T, Chen, YQ, Vinagre Jara, BM: Matrix approach to discrete fractional calculus II: partial fractional differential equations. J. Comput. Phys.. 228, 3137–3153 (2009). Publisher Full Text

22. Hanert, E: On the numerical solution of space time fractional diffusion models. Comput. Fluids. 46, 33–39 (2011). Publisher Full Text

23. Lin, R, Liu, F, Anh, V, Turner, I: Stability and convergence of a new explicit finite-difference approximation for the variable-order nonlinear fractional diffusion equation. Appl. Math. Comput.. 212, 435–445 (2009). Publisher Full Text

24. Crank, J, Nicolson, P: A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type. Proc. Camb. Philol. Soc.. 43(1), 50–67 (1947). Publisher Full Text

25. Diethelm, K, Ford, NJ, Freed, AD: Detailed error analysis for a fractional Adams method. Numer. Algorithms. 36, 31–52 (2004)

26. Li, CP, Tao, CX: On the fractional Adams method. Comput. Math. Appl.. 58, 1573–1588 (2009). Publisher Full Text