version 1.4.0.0 (5.28 KB) by
Roberto Garrappa

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.

Roberto Garrappa (2021). Predictor-corrector PECE method for fractional differential equations (https://www.mathworks.com/matlabcentral/fileexchange/32918-predictor-corrector-pece-method-for-fractional-differential-equations), MATLAB Central File Exchange. Retrieved .

Created with
R2009b

Compatible with any release

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Jadon ClugstonFahimeh AkhavanThank you very much for your response, dear professor.

Roberto Garrappa@Fahimeh Akhavan; please, just read the long documentation provided with the code. The fractional order must be a single value, thus run the code with just alpha = 0.8

Roberto Garrappa@Fahimeh Akhavan: this code does not allow to solve multi-order systems, i.e. systems in which each equation has a different order as it is your system. You can find some codes for multi-order systems in my web-page https://www.dm.uniba.it/members/garrappa/software

Fahimeh Akhavan@Roberto Garrappa: I'm trying to implement the fde12.m on a classical fractional-order dynamical system consisting of the nonlinear Brusselator system but I'm getting an error, saying ">>Matrix dimensions must agree.

>>Error in fde12 (line 149)

>>nalpha = nvett.^alpha ;".

The inputs are as follow:

alpha = [0.8,0.7];

A=1;

B=3;

param=[A,B];

f_fun = @(t,y,par) [...

par(1)-(par(2)+1)*y(1)+ y(1)^2*y(2); ...

par(2)*y(1)-y(1)^2*y(2)] ;

t0=0;

T=100;

y0=[1.2;2.8];

h=2^(-5);

[t, y] = fde12(alpha,f_fun,t0,T,y0,h,param);

Could you please help me how can I fix it?

Thanks a lot

Roberto Garrappa@Fatemeh Norouzi: I think the message is clear: the vector field you are using (the function F_FUN) has 5 rows, hence you have 5 equations. Consequently, also the the initial condition must have 5 values. Maybe, something is wrong in your definition of F_FUN

Fatemeh NorouziI'm trying to implement the fde12.m on a fractional differential system of 4 variables but I'm getting an error, saying "Size 4 of the problem as obtained from initial conditions (i.e. the number of rows of Y0) not compatible with the size 5 of the output of the vector field F_FUN".

Could you please help me how can I fix it?

Also, by which function I can plot the phase portrait and time histories (series) of each output?

Thanks

Brad HaddinI used the fde12 solver to plot bifurcation diagram for the following problem:

clear; clc;

A=1; eta=0.1; mu=0.1; u=1; d=0.01; rho=0.5; delta=0.6;

global beta;

betarange = 0.5:0.1:1; % Range for parameter a

alpha = 0.99; k = 0; t0=0; tfinal=200; h = 2^(-6) ; % Time interval for solving Rossler system

xmax = []; % A matrix for storing the sorted value of x1

for beta = betarange

f_fun = @(t,x) [A-eta*x(1)*x(2)+(mu+beta*u/(1+rho*u*x(2)))*x(2)-d*x(1); eta*x(1)*x(2)-(mu+beta*u/(1+rho*u*x(2)))*x(2)-(d+delta)*x(2)]; % Rossler system

x0 = [1.75; 1.5]; % initial condition for Rossler system

k = k + 1;

[t,x] = fde12(alpha,f_fun,t0,tfinal,x0,h); % call ode() to solve Rossler system

count = find(t>100); % find all the t_values which is >100

x = x(count,:);

j = 1;

n = length(x(:,1)); % find the length of vector x1(x in our problem)

for i = 2:n-1

% check for the min value in 1st column of sol matrix

if (x(i-1,1)+eps) < x(i,1) && x(i,1) > (x(i+1,1)+eps)

xmax(k,j)=x(i,1); % Sorting the values of x1 in increasing order

j=j+1;

end

end

% generating bifurcation map by plotting j-1 element of kth row each time

if j>1

plot(a,xmax(k,1:j-1),'b.');

end

hold on;

index(k)=j-1;

end

xlabel('Bifurcation parameter \beta');

ylabel('y max');

Still I got the following error:

Index exceeds matrix

dimensions.

Error in bifurcation_beta

(line 13)

x = x(count,:);

can anyone help me in this regard? Thanks in advance.

Ramashis Banerjee@Roberto Garrappa,Sir how to solve FODE backward in time using this code? Please help me.

Thank You

Ramashis BanerjeeSir, how to solve FODE using this code when value of alpha is greater than 1? How to do error analysis and find EOC?

Thank you in advance

Ramashis BanerjeeThank you very much sir

Roberto Garrappa@Ramashis Banerjee: if the code show a "not enough input arguments" message, it clearly means that you are not providing all the parameters which are necessary to run the code. These are not errors in the code but errors in using the code. Please, carefully read the description provided with the code or the final part of the paper https://www.mdpi.com/2227-7390/6/2/16 where some examples are reported

Ramashis BanerjeeIt is showing not enough input arguments

Roberto Garrappa@Ramashis Banerjee: it would be useful, to the author and to the community, that you report in a detailed way what kind of errors you found in the code. The FDE12.m code has been intensively tested by myself and several other colleagues and so far no errors were detected; obviously, I can not exclude the possible presence of some residual bug. Your comment however does not help to identify if the code actually has some error or if your issues are due to an improper use of the code.

Ramashis BanerjeeThere are errors in the code

Zulqr MRoberto Garrappa@Tareq, you can find Matlab routines for multi-order systems (systems in which each equation has a different fractional order) and multi-term fractional equations (equations with several derivatives of different order) in the software section of my web page https://www.dm.uniba.it/Members/garrappa/Software

Tareq AbuaishaWell done Roberto!

Do you have an idea how can the code be upgraded in order to solve an FODE with alpha and beta?

fuchang wangasma brwe cannot solve Backward Integration using fde12

is there any suggestion how we could calculate Backward Integration for fractional order nonlinear systems

fuchang wanga = 1 ; mu = 4 ;

fdefun = @(t,y) [ a-(mu+1)*y(1)+y(1)^2*y(2) ; mu*y(1)-y(1)^2*y(2) ] ;

Jfdefun = @(t,y) [ -(mu+1)+2*y(1)*y(2) , y(1)^2 ; mu-2*y(1)*y(2) , -y(1)^2 ] ;

alpha = 0.8 ;

t0 = 0 ; tfinal = 100 ; y0 = [ 0.2 ; 0.03] ;

h = 2^(-6) ;

[t, y_flmm2] = flmm2(alpha,fdefun,Jfdefun,t0,tfinal,y0,h) ;

[t, y_fde12] = fde12(alpha,fdefun,t0,tfinal,y0,h) ;

figure(1)

subplot(1,2,1),plot(t,y_flmm2(1,:),t,y_flmm2(2,:)) ;

xlabel('t') ; ylabel('y(t)') ;

legend('y_1(t)','y_2(t)') ;

title('FDE solved by the FLMM2.m code');

subplot(1,2,2),plot(y_flmm2(1,:),y_flmm2(2,:)) ;

xlabel('y1(t)') ; ylabel('y2(t)') ;

title('FDE solved by the FLMM2.m code');

figure(2)

subplot(1,2,1),plot(t,y_fde12(1,:),t,y_fde12(2,:)) ;

xlabel('t') ; ylabel('y(t)') ;

legend('y_1(t)','y_2(t)') ;

title('FDE solved by the FDE12.m code');

subplot(1,2,2),plot(y_fde12(1,:),y_fde12(2,:)) ;

xlabel('y1(t)') ; ylabel('y2(t)') ;

title('FDE solved by the FDE12.m code');

Mohamed AburakhisWhen alpha <1, it is easy to use but when alpha >1 and i have more that one state equations, how can i use it

MVOGO AlainCould you please provide an example how to use this code in the case of the fractional Lorentz system for example

James506Hi. I was interested how/if you could use this to solve a system of mixed-order FDEs? For example, if I had an FDE which also had standard first and second order derivatives in it. Only a single value can be prescribed for the equation order (alpha). Thanks

Roberto GarrappaDear Sasan Attar, the usage of this code is very similar to those of other Matlab codes for ODEs (for instance, ode15s or ode45). Anyway, on the software page of my website http://www.dm.uniba.it/Members/garrappa/Software you can find a short tutorial to the use of this code and to the use of other codes for FDEs. Thank you for your interest.

sasan attarCould you please provide an example how to use this code

Manashita Borahvery helpful

Roberto GarrappaOn the software section of my home page http://www.dm.uniba.it/Members/garrappa/ I have uploaded a short tutorial on the use of this code for solving FDEs

qiqhow to solve fde system with 2 or more equations?

amaalhow i can write the input fdefun?

MirzaCould you please provide an example how to use this code.