function [t, y] = fde12(alpha,fdefun,t0,tfinal,y0,h,param,mu,mu_tol)
%FDE12 Solves an initial value problem for a non-linear differential
% equation of fractional order (FDE). The code implements the
% predictor-corrector PECE method of Adams-Bashforth-Moulton type
% described in [1].
%
% [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.
%
% FDE12 is an implementation of the predictor-corrector method of
% Adams-Bashforth-Moulton studied 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].
%
% [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.
%
% Copyright (c) 2011-2012, Roberto Garrappa, University of Bari, Italy
% garrappa at dm dot uniba dot it
% Revision: 1.2 - Date: July, 6 2012
% Check inputs
if nargin < 9
mu_tol = 1.0e-6 ;
if nargin < 8
mu = 1 ;
if nargin < 7
param = [] ;
end
end
end
% Check order of the FDE
if alpha <= 0
error('MATLAB:fde12:NegativeOrder',...
['The order ALPHA of the FDE must be positive. The value ' ...
'ALPHA = %f can not be accepted. See FDE12.'], alpha);
end
% Check the step--size of the method
if h <= 0
error('MATLAB:fde12:NegativeStepSize',...
['The step-size H for the method must be positive. The value ' ...
'H = %e can not be accepted. See FDE12.'], h);
end
% Structure for storing initial conditions
ic.t0 = t0 ;
ic.y0 = y0 ;
ic.m_alpha = ceil(alpha) ;
ic.m_alpha_factorial = factorial(0:ic.m_alpha-1) ;
% Structure for storing information on the problem
Probl.ic = ic ;
Probl.fdefun = fdefun ;
Probl.problem_size = size(y0,1) ;
Probl.param = param ;
% Check number of initial conditions
if size(y0,2) < ic.m_alpha
error('MATLAB:fde12:NotEnoughInputs', ...
['A not sufficient number of assigned initial conditions. ' ...
'Order ALPHA = %f requires %d initial conditions. See FDE12.'], ...
alpha,ic.m_alpha);
end
% Check compatibility size of the problem with size of the vector field
f_temp = f_vectorfield(t0,y0(:,1),Probl) ;
if Probl.problem_size ~= size(f_temp,1)
error('MATLAB:fde12:SizeNotCompatible', ...
['Size %d of the problem as obtained from initial conditions ' ...
'(i.e. the number of rows of Y0) not compatible with the ' ...
'size %d of the output of the vector field FDEFUN. ' ...
'See FDE12.'], Probl.problem_size,size(f_temp,1));
end
% Number of points in which to evaluate weights and solution
r = 16 ;
N = ceil((tfinal-t0)/h) ;
Nr = ceil((N+1)/r)*r ;
Q = ceil(log2(Nr/r)) - 1 ;
NNr = 2^(Q+1)*r ;
% Preallocation of some variables
y = zeros(Probl.problem_size,N+1) ;
fy = zeros(Probl.problem_size,N+1) ;
zn_pred = zeros(Probl.problem_size,NNr+1) ;
if mu > 0
zn_corr = zeros(Probl.problem_size,NNr+1) ;
else
zn_corr = 0 ;
end
% Evaluation of coefficients of the PECE method
nvett = 0 : NNr+1 ;
nalpha = nvett.^alpha ;
nalpha1 = nalpha.*nvett ;
PC.bn = nalpha(2:end) - nalpha(1:end-1) ;
PC.an = [ 1 , (nalpha1(1:end-2) - 2*nalpha1(2:end-1) + nalpha1(3:end)) ] ;
PC.a0 = [ 0 , nalpha1(1:end-2)-nalpha(2:end-1).*(nvett(2:end-1)-alpha-1)] ;
PC.halpha1 = h^alpha/gamma(alpha+1) ;
PC.halpha2 = h^alpha/gamma(alpha+2) ;
PC.mu = mu ; PC.mu_tol = mu_tol ;
% Initializing solution and proces of computation
t = t0 + (0 : N)*h ;
y(:,1) = y0(:,1) ;
fy(:,1) = f_temp ;
[y, fy] = Triangolo(1, r-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl ) ;
% Main process of computation by means of the FFT algorithm
ff = [0 2 ] ; nx0 = 0 ; ny0 = 0 ;
for q = 0 : Q
L = 2^q ;
[y, fy] = DisegnaBlocchi(L, ff, r, Nr, nx0+L*r, ny0, t, y, fy, ...
zn_pred, zn_corr, N, PC, Probl ) ;
ff = [ff ff] ; ff(end) = 4*L ;
end
% Evaluation solution in TFINAL when TFINAL is not in the mesh
if tfinal < t(N+1)
c = (tfinal - t(N))/h ;
t(N+1) = tfinal ;
y(:,N+1) = (1-c)*y(:,N) + c*y(:,N+1) ;
end
t = t(1:N+1) ; y = y(:,1:N+1) ;
end
% =========================================================================
% =========================================================================
% r : dimension of the basic square or triangle
% L : factor of resizing of the squares
function [y, fy] = DisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, ...
zn_pred, zn_corr, N , PC, Probl)
nxi = nx0 ; nxf = nx0 + L*r - 1 ;
nyi = ny0 ; nyf = ny0 + L*r - 1 ;
is = 1 ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;
i_triangolo = 0 ; stop = 0 ;
while ~stop
stop = nxi+r-1 == nx0+L*r-1 | (nxi+r-1>=Nr-1) ; % Ci si ferma quando il triangolo attuale finisce alla fine del quadrato
[zn_pred, zn_corr] = Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl) ;
[y, fy] = Triangolo(nxi, nxi+r-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl) ;
i_triangolo = i_triangolo + 1 ;
if ~stop
if nxi+r-1 == nxf % Il triangolo finisce dove finisce il quadrato -> si scende di livello
i_Delta = ff(i_triangolo) ;
Delta = i_Delta*r ;
nxi = s_nxf(is)+1 ; nxf = s_nxf(is) + Delta ;
nyi = s_nxf(is) - Delta +1; nyf = s_nxf(is) ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;
else % Il triangolo finisce prima del quadrato -> si fa un quadrato accanto
nxi = nxi + r ; nxf = nxi + r - 1 ; nyi = nyf + 1 ; nyf = nyf + r ;
is = is + 1 ;
s_nxi(is) = nxi ; s_nxf(is) = nxf ; s_nyi(is) = nyi ; s_nyf(is) = nyf ;
end
end
end
end
% =========================================================================
% =========================================================================
function [zn_pred, zn_corr] = Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl)
coef_beg = nxi-nyf ; coef_end = nxf-nyi+1 ;
funz_beg = nyi+1 ; funz_end = nyf+1 ;
% Evaluation convolution segment for the predictor
vett_coef = PC.bn(coef_beg:coef_end) ;
vett_funz = [fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
zzn_pred = real(FastConv(vett_coef,vett_funz)) ;
zn_pred(:,nxi+1:nxf+1) = zn_pred(:,nxi+1:nxf+1) + zzn_pred(:,nxf-nyf+1-1:end-1) ;
% Evaluation convolution segment for the corrector
if PC.mu > 0
vett_coef = PC.an(coef_beg:coef_end) ;
if nyi == 0 % Evaluation of the lowest square
vett_funz = [zeros(Probl.problem_size,1) , fy(:,funz_beg+1:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
else % Evaluation of any square but not the lowest
vett_funz = [ fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg+1) ] ;
end
zzn_corr = real(FastConv(vett_coef,vett_funz)) ;
zn_corr(:,nxi+1:nxf+1) = zn_corr(:,nxi+1:nxf+1) + zzn_corr(:,nxf-nyf+1:end) ;
else
zn_corr = 0 ;
end
end
% =========================================================================
% =========================================================================
function [y, fy] = Triangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, PC, Probl)
for n = nxi : min(N,nxf)
% Evaluation of the predictor
Phi = zeros(Probl.problem_size,1) ;
if nxi == 1 % Case of the first triangle
j_beg = 0 ;
else % Case of any triangle but not the first
j_beg = nxi ;
end
for j = j_beg : n-1
Phi = Phi + PC.bn(n-j)*fy(:,j+1) ;
end
St = StartingTerm(t(n+1),Probl.ic) ;
y_pred = St + PC.halpha1*(zn_pred(:,n+1) + Phi) ;
f_pred = f_vectorfield(t(n+1),y_pred,Probl) ;
if PC.mu == 0
y(:,n+1) = y_pred ;
fy(:,n+1) = f_pred ;
else
j_beg = nxi ;
Phi = zeros(Probl.problem_size,1) ;
for j = j_beg : n-1
Phi = Phi + PC.an(n-j+1)*fy(:,j+1) ;
end
Phi_n = St + PC.halpha2*(PC.a0(n+1)*fy(:,1) + zn_corr(:,n+1) + Phi) ;
yn0 = y_pred ; fn0 = f_pred ;
stop = 0 ; mu_it = 0 ;
while ~stop
yn1 = Phi_n + PC.halpha2*fn0 ;
mu_it = mu_it + 1 ;
if PC.mu == Inf
stop = norm(yn1-yn0,inf) < PC.mu_tol ;
if mu_it > 100 && ~stop
warning('MATLAB:fde12:NonConvegence',...
strcat('It has been requested to run corrector iterations until convergence but ', ...
'the process does not converge to the tolerance %e in 100 iteration'),PC.mu_tol) ;
stop = 1 ;
end
else
stop = mu_it == PC.mu ;
end
fn1 = f_vectorfield(t(n+1),yn1,Probl) ;
yn0 = yn1 ; fn0 = fn1 ;
end
y(:,n+1) = yn1 ;
fy(:,n+1) = fn1 ;
end
end
end
% =========================================================================
% =========================================================================
function z = FastConv(x, y)
r = length(x) ; problem_size = size(y,1) ;
z = zeros(problem_size,r) ;
X = fft(x,r) ;
for i = 1 : problem_size
Y = fft(y(i,:),r) ;
Z = X.*Y ;
z(i,:) = ifft(Z,r) ;
end
end
% =========================================================================
% =========================================================================
function f = f_vectorfield(t,y,Probl)
if isempty(Probl.param)
f = feval(Probl.fdefun,t,y) ;
else
f = feval(Probl.fdefun,t,y,Probl.param) ;
end
end
% =========================================================================
% =========================================================================
function ys = StartingTerm(t,ic)
ys = zeros(size(ic.y0,1),1) ;
for k = 1 : ic.m_alpha
ys = ys + (t-ic.t0)^(k-1)/ic.m_alpha_factorial(k)*ic.y0(:,k) ;
end
end
function f=LE_RF(~,x)
%Output data must be a column vector
f=zeros(size(x));
%variables allocated to the variational equations
X= [x(4), x(7), x(10);
x(5), x(8), x(11);
x(6), x(9), x(12)];
%RF equations
f(1)=x(2)*(x(3)-1+x(1)*x(1))+0.1*x(1);
f(2)=x(1)*(3*x(3)+1-x(1)*x(1))+0.1*x(2);
f(3)=-2*x(3)*(0.98+x(1)*x(2));
%Jacobian matrix
J=[2*x(1)*x(2)+0.1, x(1)*x(1)+x(3)-1, x(2);
-3*x(1)*x(1)+3*x(3)+1,0.1,3*x(1);
-2*x(2)*x(3),-2*x(1)*x(3),-2*(x(1)*x(2)+0.98)];
%Righthand side of variational equations
f(4:12)=J*X;
end
function [t,LE]=FO_Lyapunov(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,out)
ne=3;
ext_fcn=@LE_RF;
t_start=0;
h_norm=0.01;
t_end=10;
x_start=[0.1;0.1;0.1];
h=0.01;
q=0.75;
out=0;
% [t,LE]=FO_Lyapunov(3,@LE_RF,0,0.02,300,[0.1;0.1;0.1],0.02,0.98,1000);
%
% Program to compute LEs as function of time (t) of autonomous systems of commensurate Fractional Order (FO)
% defined via Caputo's derivative;
%
% ***********************************************************************************************
% August 2022: Improved to overcome the problem of plot function of the last
% variants of Matlab
% ***********************************************************************************************
%
% author: Marius-F. Danca, 2018
% http://danca.rist.ro/
% email: danca@rist.ro
% The program uses a fast variant of the predictor-corrector Adams-Bashforth-Moulton,
% FDE12.m for FDEs, by Roberto Garrappa: https://goo.gl/XScYmi
% m-files required: FDE12, FO_Lyapunov.m and the function containing the extended system
% (see e.g. LE_RF.m).
% FO_Lyapunov.m, FDE12.m and LE_FO_RF.m must be in the same folder.
% FO_Lyapunov.m was developed, by modifying the program Lyapunov.m, by V.N. Govorukhin:
% https://goo.gl/wZVCtg
%
% How to use it:
% [t,LE]=FO_Lyapunov(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,out);
%
% Input:
% ne - system dimension (eqs number);
% ext_fcn - function containing the extended system of FO;
% t_start, t_end - time span (FDE12);
% h - integration step;
% h_norm - step for Gram-Schmidt renormalization (Lyapunov algorithm), an important parameter:
% experimentaly determined as an optimal choice if it is multiple of h, or equal;
% x_start - initial condition;
% q - the fractional order;
% out - priniting step of LEs values;
% out==0 - no print.
%
% Output:
% t - time values;
% LE Lyapunov exponents to each time value printed every 'out' step.
%
% Example of use for the RF system:
% [t,LE]=FO_Lyapunov(3,@LE_RF,0,0.02,300,[0.1;0.1;0.1],0.02,0.98,1000);
%
% Cite as:
%
% Marius-F. Danca and N. Kuznetsov, Matlab code for Lyapunov exponents of
% fractional order systems, International Journal of Bifurcation and Chaos,
% 28(05)(2018), 1850067
%
% see also:
% Marius.-F. Danca, Matlab code for Lyapunov exponents of fractional-order systems,
% Part II: The non-commensurate case, International Journal of Bifurcation and Chaos,
% 31(12), 2150187, (2021)
%
hold on;
colors = 'grby';%can be expanded expand that to include any of the 8 color codes bcgkmrwy
% Memory allocation
x=zeros(ne*(ne+1),1);
x0=x;
c=zeros(ne,1);
gsc=c; zn=c;
n_it = round((t_end-t_start)/h_norm);
% Initial values
x(1:ne)=x_start;
i=1;
while i<=ne
x((ne+1)*i)=1.0;
i=i+1;
end
t=t_start;
% Main loop
it=1;
while it<=n_it
LE=zeros(ne,1);
% Solutuion of extended ODE system of FO using FDE12 routine
[~,Y] = fde12(q,ext_fcn,t,t+h_norm,x,h);
t=t+h_norm;
Y=transpose(Y);
x=Y(size(Y,1),:); %solution at t+h_norm
i=1;
while i<=ne
j=1;
while j<=ne
x0(ne*i+j)=x(ne*j+i);
j=j+1;
end
i=i+1;
end
% construct new orthonormal basis by gram-schmidt
zn(1)=0.0;
j=1;
while j<=ne
zn(1)=zn(1)+x0(ne*j+1)^2;
j=j+1;
end
zn(1)=sqrt(zn(1));
j=1;
while j<=ne
x0(ne*j+1)=x0(ne*j+1)/zn(1);
j=j+1;
end
j=2;
while j<=ne
k=1;
while k<=j-1
gsc(k)=0.0;
l=1;
while l<=ne
gsc(k)=gsc(k)+x0(ne*l+j)*x0(ne*l+k);
l=l+1;
end
k=k+1;
end
k=1;
while k<=ne
l=1;
while l<=j-1
x0(ne*k+j)=x0(ne*k+j)-gsc(l)*x0(ne*k+l);
l=l+1;
end
k=k+1;
end
zn(j)=0.0;
k=1;
while k<=ne
zn(j)=zn(j)+x0(ne*k+j)^2;
k=k+1;
end
zn(j)=sqrt(zn(j));
k=1;
while k<=ne
x0(ne*k+j)=x0(ne*k+j)/zn(j);
k=k+1;
end
j=j+1;
end
% update running vector magnitudes
k=1;
while k<=ne
c(k)=c(k)+log(zn(k));
k=k+1;
end
% normalize exponent
k=1;
while k<=ne
LE(k)=c(k)/(t-t_start);
k=k+1;
end
% print results
if (mod(it,out)==0)
fprintf('%10.2f %10.4f %10.4f %10.4f\n',[t,LE']);
end
i=1;
while i<=ne
j=1;
while j<=ne
x(ne*j+i)=x0(ne*i+j);
j=j+1;
end
i=i+1;
end
x=transpose(x);
it=it+1;
% display LEs
for i=1:ne
plot(t,LE(i),'.','Color',colors(i));
%plot(t,LE(i),'.','markersize',1,'Color',colors(i));%for thiny draw
end
end
xlabel('t','fontsize',16)
ylabel('LEs','fontsize',14)
set(gca,'fontsize',14)
box on;
line([0,t],[0,0],'color','k')
axis tight