An Explicit Numerical Method for the Fractional Cable Equation

An explicit numerical method to solve a fractional cable equation which involves two temporal Riemann-Liouville derivatives is studied. The numerical difference scheme is obtained by approximating the first-order derivative by a forward difference formula, the Riemann-Liouville derivatives by the Grünwald-Letnikov formula, and the spatial derivative by a three-point centered formula. The accuracy, stability, and convergence of the method are considered. The stability analysis is carried out by means of a kind of von Neumann method adapted to fractional equations. The convergence analysis is accomplished with a similar procedure. The von-Neumann stability analysis predicted very accurately the conditions under which the present explicit method is stable. This was thoroughly checked by means of extensive numerical integrations.


Introduction
Fractional calculus is a key tool for solving some relevant scientific problems in physics, engineering, biology, chemistry, hydrology, and so on 1-6 .A field of research in which the fractional formalism has been particularly useful is that related to anomalous diffusion processes 1, 7-13 .This kind of process is singularly abundant and important in biological media 14-16 .In this context, the electrodiffusion of ions in neurons is an anomalous diffusion problem to which the fractional calculus has recently been applied.The precise origin of the anomalous character of this diffusion process is not clear see 17 and references therein , but in any case the consideration of anomalous diffusion in the modeling of electrodiffusion of ions in neurons seems pertinent.This problem has been addressed recently by Langlands et al. 17,18 .An equation that plays a key role in their analysis is the following fractional cable or telegrapher's or Cattaneo equation model II : where with n − 1 < γ < n and n integer, is the Riemann-Liouville fractional derivative.Here u is the difference between the membrane potential and the resting membrane potential, γ 1 is the exponent characterizing the anomalous flux of ions along the nerve cell, and γ 2 is the exponent characterizing the anomalous flux across the membrane 17, 18 .Some earlier fractional cable equations were discussed in 19, 20 .A variety of analytical and numerical methods to solve many classes of fractional equations have been proposed and studied over the last few years 10, 21-30 .Of the numerical methods, finite difference methods have been particularly fruitful 31-38 .These methods can be broadly classified as explicit or implicit 39 .An implicit method for dealing with 1.1 has recently been considered by Liu et al. 38 .Although implicit methods are more cumbersome than explicit methods, they usually remain stable over a larger range of parameters, especially for large timesteps, which makes them particularly suitable for fractional diffusion problems.Nevertheless, explicit methods have some features that make them widely appreciated 32, 39 : flexibility, simplicity, small computational demand, and easy generalization to spatial dimensions higher than one.Unfortunately, they can become unstable in some cases, so that it is of great importance to determine the conditions under which these methods are stable.In this paper we will discuss an explicit finite difference scheme for solving the fractional cable equation, which is close to the methods studied in 32, 33 .We shall address two main questions: i whether this kind of method can cope with fractional equations involving different fractional derivatives, such as the fractional cable equation; ii whether the von Neumann stability analysis put forward in 32, 34 is suitable for this kind of equation.

The Numerical Method
Henceforth, we will use the notation x j jΔx, t m mΔt, and u x j , t m u m j U m j , where U m j is the numerical estimate of the exact solution u x, t at x x j and t t m .
In order to get the numerical difference algorithm, we discretize the continuous differential and integro-differential operators as follows.For the discretization of the fractional Riemann-Liouville derivative we use the Gr ünwald-Letnikov formula and ω α 0 1.These coefficients come from the generating function 40 To discretize the integer derivatives we use standard formulas: for the second-order spatial derivative we employ the three-point centered formula with and for the first-order time derivative we use the forward derivative where Inserting 2.1 , 2.5 , and 2.7 into 1.1 , one gets where, as can easily be proved, the truncating error T x, t is

2.10
Neglecting the truncating error we get the finite difference scheme we are seeking: where

2.13
To test this algorithm, we solved 1.1 in the interval −L/2 ≤ x ≤ L/2, with absorbing boundary conditions, u x −L/2, t u x L/2, t 0, and initial condition given by a Dirac's delta function centered at x 0: u x, 0 δ x .The exact solution of this problem for L → ∞ is 17 where H denotes the Fox H function 10, 41 .In our numerical procedure, the exact initial condition u x, 0 δ x is approximated by

2.15
The explicit difference scheme 2.12 is tested by comparing the analytical solution with the numerical solution for several cases of the problem described following 2.13 with different values of γ 1 and γ 2 .We have computed the analytical solution by means of 2.14 truncating the series at k 20.The corresponding Fox H function was evaluated by means of the series expansion described in 10, 41 truncating the infinite series after the first 50 terms.In Figures 1 and 2 we show the analytical and numerical solution for two values of γ 1 and γ 2 at x 0 and x 0.5.The differences between the exact and the numerical solution are shown in Figures 3  and 4. One sees that, except for very short times, the agreement is quite good.The large value of the error for small times is due in part to the approximation of the Dirac's delta function at x 0 by 2.15 .This is clearly appreciated when noticing the quite different scales of Figures 3 and 4: the error is much smaller for x 0.5 than for x 0. For the cases with γ 1 1/2 we used a smaller value of Δt and, simultaneously, a larger value of Δx than for the cases with γ 1 1 in order to keep the numerical scheme stable.This issue will be discussed in Section 3.

Stability
As usual for explicit methods, the present explicit difference scheme 2.12 is not unconditionally stable, that is, for any given set of values of γ 1 , γ 2 , μ, and K there are choices of Δx and Δt for which the method is unstable.Therefore, it is important to determine the conditions under which the method is stable.To this end, here we shall employ the fractional von Neumann stability analysis or Fourier analysis put forward in 32 see also [33][34][35] .
The question we address is to what extent this procedure is valid for fractional diffusion equations that involve fractional derivatives of different order.Proceeding as 32 , we start by recognizing that the solution of our problem can be written as the linear combination of subdiffusive modes, u m j q ζ m q e iqjΔx , where the sum is over all the wave numbers q supported by the lattice.Therefore, following the von Neumann ideas, we reduce the problem of analyzing the stability of the solution to the problem of analyzing the stability of a single generic subdiffusion mode, ζ m e iqjΔx .Inserting this expression into 2.12 one gets The stability of the mode is determined by the behavior of ζ m .Writing of the numerical method for the problem considered in Figures 1 and 2 at x 0.5.Squares: γ 2 1; circles: γ 2 1/2; hollow symbols: γ 1 1; filled symbols: γ 1 1/2.and assuming that the amplification factor ξ of the subdiffusive mode is independent of time, we get If |ξ| > 1 for some q, the temporal factor of the solution grows to infinity c.f., 3.2 , and the mode is unstable.Considering the extreme value ξ −1, one gets from 3.3 that the numerical method is stable if this inequality holds: International Journal of Differential Equations 7 where S S sin 2 qΔx 2 .

3.5
If one defines S × lim m → ∞ S m × , one gets Therefore, because S ≤ S, we find that a sufficient condition for the present method to be stable is that S ≤ S × .In Figures 5 and 6 we show two representative examples of the problem considered in Figure 2 but for two values of S, respectively, larger and smaller than the stability bound provided by 3.7 .One sees that the value of S that one chooses is crucial: when S is smaller than S × one is inside the stable region and gets a sensible numerical solution Figure 5 ; otherwise one gets an evidently unstable and nonsensical solution Figure 6 .

Numerical Check of the Stability Analysis
In this section we describe a comprehensive check of the validity of our stability analysis by using many different values of the parameters γ 1 , γ 2 , Δt, and Δx and testing whether the stability of the numerical method is as predicted by 3.7 .Without loss of generality, we assume μ K 1 in all cases.We proceed in the following way.First, we choose a set of values of γ 1 , γ 2 , Δx, and S and integrate the corresponding fractional cable equation.If for λ 10 within the first 1000 integrations, then we say the method is unstable; otherwise, we label the method as stable.We generated Figure 7 by starting the integration for values of S well below the theoretical stability limit given by 3.7 and kept increasing its value by 0.001 until condition 4.1 was first reached.The last value for which the method was stable is recorded and plotted in Figure 7.The limit value λ 10 is arbitrary, but choosing any other reasonable value does not significantly change these plots.

Convergence Analysis
In this section we show that the present numerical method is convergent, that is, that the numerical solution converges towards the exact solution when the size of the spatiotemporal where |ζ {m} | is the maximum value of |ζ k | for k 0, . . ., m. Taking into account 2.4 , using the value z 1, and because

Conclusions
An explicit method for solving a kind of fractional diffusion equation that involves several fractional Riemann-Liouville derivatives, which are approximated by means of the Gr ünwald-Letnikov formula, has been considered.The method was used to solve a class of equations of this type fractional cable equations with free boundary conditions, Dirac's delta initial condition, and different fractional exponents.The error of the numerical method is compatible with the truncating error, which is of order O Δt O Δx 2 .It was also proved that the method is convergent.Besides, it was also found that a fractional von-Neumann stability analysis, which provides very precise stability conditions for standard fractional diffusion equations, leads also to a very accurate estimate of the stability conditions for cable equations.

3 . Therefore m k 0 |ω 1 −γ k | is bounded in 10 International
Journal of Differential Equations fact, it is smaller than 2 .Using this result in 5.3 , together with |ζ k | ≤ C Δt Δx 2 and |χ k | ≤ C Δt Δx 2 , we find thatζ m 1 ≤ C Δt Δx 2 .5.4Therefore the amplitude of the subdiffusive modes goes to zero when the spatiotemporal mesh goes to zero.Employing the Parseval relation, this means that the norm of the error e k 2 ≡ j |e when Δt and Δx go to zero.This is what we aimed to prove.