Predictor-corrector PECE method for fractional differential equations

Solves initial value problems for fractional differential equations

FDE12 solves an initial value problem for a non-linear differential equation of fractional order (FDE).

This is an implementation of the predictor-corrector method of Adams-Bashforth-Moulton described in [1]. Convergence and accuracy of the method are studied in [2]. The implementation with multiple corrector iterations has been proposed and discussed for multiterm FDEs in [3]. In this implementation the discrete convolutions are evaluated by means of the FFT algorithm described in [4] allowing to keep the computational cost proportional to N*log(N)^2 instead of N^2 as in the classical implementation; N is the number of time-point in which the solution is evaluated, i.e. N = (TFINAL-T)/H. The stability properties of the method implemented by FDE12 have been studied in [5].

USAGE:

[T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,h) integrates the initial value problem for the FDE, or the system of FDEs, of order ALPHA > 0
D^ALPHA Y(t) = FDEFUN(T,Y(T))
Y^(k)(T0) = Y0(:,k+1), k=0,...,m-1
where m is the smallest integer grater than ALPHA and D^ALPHA is the fractional derivative according to the Caputo's definition. FDEFUN is a function handle corresponding to the vector field of the FDE and for a scalar T and a vector Y, FDEFUN(T,Y) must return a column vector. The set of initial conditions Y0 is a matrix with a number of rows equal to the size of the problem (hence equal to the number of rows of the output of FDEFUN) and a number of columns depending on ALPHA and given by m. The step-size H>0 is assumed constant throughout the integration.

[T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM) solves as above with the additional set of parameters for the FDEFUN as FDEFUN(T,Y,PARAM).

[T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU) solves the FDE with the selected number MU of multiple corrector iterations. The following values for MU are admissible:
MU = 0 : the corrector is not evaluated and the solution is provided just by the predictor method (the first order rectangular rule);
MU > 0 : the corrector is evaluated by the selected number MU of times; the classical PECE method is obtained for MU=1;
MU = Inf : the corrector is evaluated for a certain number of times until convergence of the iterations is reached (for convergence the difference between two consecutive iterates is tested).
The defalut value for MU is 1

[T,Y] = FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU,MU_TOL) allows to specify the tolerance for testing convergence when MU = Inf. If not specified, the default value MU_TOL = 1.E-6 is used.

REFERENCES:

[1] K. Diethelm, A.D. Freed, The Frac PECE subroutine for the numerical solution of differential equations of fractional order, in: S. Heinzel, T. Plesser (Eds.), Forschung und Wissenschaftliches Rechnen 1998, Gessellschaft fur Wissenschaftliche Datenverarbeitung, Gottingen, 1999, pp. 57-71.

[2] K. Diethelm, N.J. Ford, A.D. Freed, Detailed error analysis for a fractional Adams method, Numer. Algorithms 36 (1) (2004) 31-52.

[3] K. Diethelm, Efficient solution of multi-term fractional differential equations using P(EC)mE methods, Computing 71 (2003), pp. 305-319.

[4] E. Hairer, C. Lubich, M. Schlichte, Fast numerical solution of nonlinear Volterra convolution equations, SIAM J. Sci. Statist. Comput. 6 (3) (1985) 532-541.

[5] R. Garrappa, On linear stability of predictor-corrector algorithms for fractional differential equations, Internat. J. Comput. Math. 87 (10) (2010) 2281-2290.