Usage and Applications
Example 1: Fractional nonlinear equation
\[D^\beta y(t) = \frac{40320}{\Gamma(9 - \beta)} t^{8 - \beta} - 3 \frac{\Gamma(5 + \beta / 2)}{\Gamma(5 - \beta / 2)} t^{4 - \beta / 2} + \frac{9}{4} \Gamma(\beta + 1) + \left( \frac{3}{2} t^{\beta / 2} - t^4 \right)^3 - \left[ y(t) \right]^{3 / 2}\]
For $0<\beta\leq1$ being subject to the initial condition $y(0)=0$, the exact solution is:
\[y(t)=t^8-3t^{4+\beta/2}+9/4t^\beta\]
# Inputs
+tSpan = [0, 1]; # [intial time, final time]
+y0 = 0; # initial value
+β = 0.9; # order of the derivative
+
+# ODE Model
+par = β;
+F(t, y, par) = (40320 ./ gamma(9 - par) .* t .^ (8 - par) .- 3 .* gamma(5 + par / 2)
+ ./ gamma(5 - par / 2) .* t .^ (4 - par / 2) .+ 9/4 * gamma(par + 1) .+
+ (3 / 2 .* t .^ (par / 2) .- t .^ 4) .^ 3 .- y .^ (3 / 2));
+
+## Numerical solution
+t, Yapp = FDEsolver(F, tSpan, y0, β, par);
+
+# Plot
+plot(t, Yapp, linewidth = 5, title = "Solution of a 1D fractional IVP",
+ xaxis = "Time (t)", yaxis = "y(t)", label = "Approximation");
+plot!(t, t -> (t.^8 - 3 * t .^ (4 + β / 2) + 9/4 * t.^β),
+ lw = 3, ls = :dash, label = "Exact solution");
Example 2: Lotka-volterra-predator-prey
\[\begin{align*} +D^{μ_1}u &= (\alpha - \beta v) u = \alpha u - \beta u v \\ +D^{μ_2}v &= (-\gamma + \delta u) v = -\gamma v + \delta u v +\end{align*}\]
# Inputs
+tSpan = [0, 25]; # [initial time, final time]
+y0 = [34, 6]; # initial values
+μ = [0.98, 0.99]; # order of derivatives
+par = [0.55, 0.028, 0.84, 0.026]; # model parameters
+
+# ODE Model
+function F(t, y, par)
+
+ α = par[1] # growth rate of the prey population
+ β = par[2] # rate of shrinkage relative to the product of the population sizes
+ γ = par[3] # shrinkage rate of the predator population
+ δ = par[4] # growth rate of the predator population as a factor of the product
+ # of the population sizes
+
+ u = y[1] # population size of the prey species at time t[n]
+ v = y[2] # population size of the predator species at time t[n]
+
+ F1 = α .* u .- β .* u .* v
+ F2 = - γ .* v .+ δ .* u .* v
+
+ [F1, F2]
+
+end
+
+## Solution
+t, Yapp = FDEsolver(F, tSpan, y0, μ, par);
+
+# Plot
+plot(t, Yapp, linewidth = 5, title = "Solution to LV model with 2 FDEs",
+ xaxis = "Time (t)", yaxis = "y(t)", label = ["Prey" "Predator"]);
+plot!(legendtitle = "Population of");
Example 3: SIR model
One application of using fractional calculus is taking into account effects of memory in modeling including epidemic evolution.
\[\begin{align*} +D^{α_1}S &= -\beta IS, \\ +D^{α_2}I &= \beta IS - \gamma I, \\ +D^{α_3}R &= \gamma I. +\end{align*}\]
By defining the Jacobian matrix, the user can achieve a faster convergence based on the modified Newton–Raphson method.
# Inputs
+I0 = 0.001; # intial value of infected
+tSpan = [0, 100]; # [intial time, final time]
+y0 = [1 - I0, I0, 0]; # initial values [S0,I0,R0]
+α = [1, 1, 1]; # order of derivatives
+h = 0.1; # step size of computation (default = 0.01)
+par = [0.4, 0.04]; # parameters [β, recovery rate]
+
+## ODE model
+function F(t, y, par)
+
+ # parameters
+ β = par[1] # infection rate
+ γ = par[2] # recovery rate
+
+ S = y[1] # Susceptible
+ I = y[2] # Infectious
+ R = y[3] # Recovered
+
+ # System equation
+ dSdt = - β .* S .* I
+ dIdt = β .* S .* I .- γ .* I
+ dRdt = γ .* I
+
+ return [dSdt, dIdt, dRdt]
+
+end
+
+## Jacobian of ODE system
+function JacobF(t, y, par)
+
+ # parameters
+ β = par[1] # infection rate
+ γ = par[2] # recovery rate
+
+ S = y[1] # Susceptible
+ I = y[2] # Infectious
+ R = y[3] # Recovered
+
+ # System equation
+ J11 = - β * I
+ J12 = - β * S
+ J13 = 0
+ J21 = β * I
+ J22 = β * S - γ
+ J23 = 0
+ J31 = 0
+ J32 = γ
+ J33 = 0
+
+ J = [J11 J12 J13
+ J21 J22 J23
+ J31 J32 J33]
+
+ return J
+
+end
+
+## Solution
+t, Yapp = FDEsolver(F, tSpan, y0, α, par, JF = JacobF, h = h);
+
+# Plot
+plot(t, Yapp, linewidth = 5, title = "Numerical solution of SIR model",
+ xaxis = "Time (t)", yaxis = "SIR populations", label = ["Susceptible" "Infectious" "Recovered"]);
Example 4: Dynamics of interaction of N species microbial communities
The impact of ecological memory on the dynamics of interacting communities can be quantified by solving fractional form ODE systems.
\[D^{β_i}X_i = X_i \left( b_i f_i(\{X_k\}) - k_i X_i \right), \quad +f_i(\{X_k\}) = \prod_{\substack{k=1 \\ k \neq i}}^N \frac{K_{ik}^n}{K_{ik}^n + X_k^n}.\]
tSpan = [0, 50]; # time span
+h = 0.1; # time step
+N = 20; # number of species
+β = ones(N); # order of derivatives
+X0 = 2 * rand(N); # initial abundances
+
+# parametrisation
+par = [5,
+ rand(N),
+ rand(N),
+ 2 * rand(N, N),
+ N];
+
+# ODE model
+function F(t, x, par)
+
+ l = par[1] # Hill coefficient
+ b = par[2] # growth rates
+ k = par[3] # death rates
+ K = par[4] # inhibition matrix
+ N = par[5] # number of species
+
+ Fun = zeros(N)
+
+ for i in 1:N
+
+ # inhibition functions
+ f = prod(K[i, 1:end .!= i] .^ l ./
+ (K[i, 1:end .!= i] .^ l .+ x[ 1:end .!= i] .^l))
+
+ # System of equations
+ Fun[i] = x[ i] .* (b[i] .* f .- k[i] .* x[ i])
+
+ end
+
+ return Fun
+
+end
+
+# Solution
+t, Xapp = FDEsolver(F, tSpan, X0, β, par, h = h, nc = 3, tol = 10e-9);
+
+# Plot
+plot(t, Xapp, linewidth = 5,
+ title = "Dynamics of microbial interaction model",
+ xaxis = "Time (t)");
+ yaxis!("Log abundance", :log10, minorgrid = true);
Example 5: Fitting orders of a COVID-19 model
Different methods are used to adjust the order of fractional differential equation models, which helps in analyzing systems across various fields.
We use Optim.jl to demonstrate how modifying system parameters and the order of derivatives in FdeSolver can enhance the fitting of COVID-19 data. The model is as follows:
\[\begin{align*} +{D}_t^{\alpha_S} S(t) =& -\beta \frac{I}{N} S - l\beta \frac{H}{N} S - \beta' \frac{P}{N} S, \\ +{D}_t^{\alpha_E} E(t) =& \beta \frac{I}{N} S + l\beta \frac{H}{N} S + \beta' \frac{P}{N} S - \kappa E, \\ +{D}_t^{\alpha_I} I(t) =& \kappa \rho_1 E - (\gamma_a + \gamma_i) I - \delta_i I, \\ +{D}_t^{\alpha_P} P(t) =& \kappa \rho_2 E - (\gamma_a + \gamma_i) P - \delta_p P, \\ +{D}_t^{\alpha_A} A(t) =& \kappa (1 - \rho_1 - \rho_2) E, \\ +{D}_t^{\alpha_H} H(t) =& \gamma_a (I + P) - \gamma_r H - \delta_h H, \\ +{D}_t^{\alpha_R} R(t) =& \gamma_i (I + P) + \gamma_r H, \\ +{D}_t^{\alpha_F} F(t) =& \delta_i I + \delta_p P + \delta_h H, +\end{align*}\]
- Model M1: fits one parameter and uses integer orders.
- Model Mf1: fits one parameter, but adjusts the derivative orders; however, all orders are equal, representing a commensurate fractional order.
- Model Mf8: fits one parameter and allows for eight distinct derivative orders, accommodating incommensurate orders for more flexibility in modeling.
# Dataset subset
+repo=HTTP.get("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"); # dataset of Covid from CSSE
+dataset_CC = CSV.read(repo.body, DataFrame); # all data of confirmed
+Confirmed=dataset_CC[dataset_CC[!,2].=="Portugal",45:121]; #comulative confirmed data of Portugal from 3/2/20 to 5/17/20
+C=diff(Float64.(Vector(Confirmed[1,:])));# Daily new confirmed cases
+
+#preprocessing (map negative values to zero and remove outliers)
+₋Ind=findall(C.<0);
+C[₋Ind].=0.0;
+outlier=findall(C.>1500);
+C[outlier]=(C[outlier.-1]+C[outlier.-1])/2;
+
+## System definition
+
+# parameters
+β=2.55; # Transmission coefficient from infected individuals
+l=1.56; # Relative transmissibility of hospitalized patients
+β′=7.65; # Transmission coefficient due to super-spreaders
+κ=0.25; # Rate at which exposed become infectious
+ρ₁=0.58; # Rate at which exposed people become infected I
+ρ₂=0.001; # Rate at which exposed people become super-spreaders
+γₐ=0.94; # Rate of being hospitalized
+γᵢ=0.27; # Recovery rate without being hospitalized
+γᵣ=0.5; # Recovery rate of hospitalized patients
+δᵢ=1/23; # Disease induced death rate due to infected class
+δₚ=1/23; # Disease induced death rate due to super-spreaders
+δₕ=1/23; # Disease induced death rate due to hospitalized class
+# Define SIR model
+function SIR(t, u, par)
+ # Model parameters.
+ N, β, l, β′, κ, ρ₁, ρ₂, γₐ, γᵢ, γᵣ, δᵢ, δₚ, δₕ=par
+
+ # Current state.
+ S, E, I, P, A, H, R, F = u
+
+# ODE
+ dS = - β * I * S/N - l * β * H * S/N - β′* P * S/N # susceptible individuals
+ dE = β * I * S/N + l * β * H * S/N + β′ *P* S/N - κ * E # exposed individuals
+ dI = κ * ρ₁ * E - (γₐ + γᵢ )*I - δᵢ * I #symptomatic and infectious individuals
+ dP = κ* ρ₂ * E - (γₐ + γᵢ)*P - δₚ * P # super-spreaders individuals
+ dA = κ *(1 - ρ₁ - ρ₂ )* E # infectious but asymptomatic individuals
+ dH = γₐ *(I + P ) - γᵣ *H - δₕ *H # hospitalized individuals
+ dR = γᵢ * (I + P ) + γᵣ* H # recovery individuals
+ dF = δᵢ * I + δₚ* P + δₕ *H # dead individuals
+ return [dS, dE, dI, dP, dA, dH, dR, dF]
+end;
+
+#initial conditions
+N=10280000/875; # Population Size
+S0=N-5; E0=0; I0=4; P0=1; A0=0; H0=0; R0=0; F0=0;
+X0=[S0, E0, I0, P0, A0, H0, R0, F0]; # initial values
+tspan=[1,length(C)]; # time span [initial time, final time]
+
+par=[N, β, l, β′, κ, ρ₁, ρ₂, γₐ, γᵢ, γᵣ, δᵢ, δₚ, δₕ]; # parameters
+
+## optimazation of β for integer order model
+
+function loss_1(b) # loss function
+ par[2]=b[1]
+ _, x = FDEsolver(SIR, tspan, X0, ones(8), par, h = .1)
+ appX=vec(sum(x[1:10:end,[3,4,6]], dims=2))
+ rmsd(C, appX; normalize=:true) # Normalized root-mean-square error
+end;
+
+p_lo_1=[1.4]; #lower bound for β
+p_up_1=[4.0]; # upper bound for β
+p_vec_1=[2.5]; # initial guess for β
+Res1=optimize(loss_1,p_lo_1,p_up_1,p_vec_1,Fminbox(BFGS()),# Broyden–Fletcher–Goldfarb–Shanno algorithm
+# Result=optimize(loss_1,p_lo_1,p_up_1,p_vec_1,SAMIN(rt=.99), # Simulated Annealing algorithm (sometimes it has better perfomance than (L-)BFGS)
+ Optim.Options(outer_iterations = 10,
+ iterations=1000,
+ show_trace=false, # turn it true to see the optimization
+ show_every=1));
+p1=vcat(Optim.minimizer(Res1));
+par1=copy(par); par1[2]=p1[1];
+
+## optimazation of β and order of commensurate fractional order model
+function loss_F_1(pμ)
+ par[2] = pμ[1] # infectivity rate
+ μ = pμ[2] # order of derivatives
+
+ _, x = FDEsolver(SIR, tspan, X0, μ*ones(8), par, h = .1)
+ appX=vec(sum(x[1:10:end,[3,4,6]], dims=2))
+ rmsd(C, appX; normalize=:true)
+end;
+
+p_lo_f_1=vcat(1.4,.5); # lower bound for β and orders
+p_up_f_1=vcat(4,1); # upper bound for β and orders
+p_vec_f_1=vcat(2.5,.9); # initial guess for β and orders
+ResF1=optimize(loss_F_1,p_lo_f_1,p_up_f_1,p_vec_f_1,Fminbox(LBFGS()), # LBFGS is suitable for large scale problems
+# Result=optimize(loss,p_lo,p_up,pvec,SAMIN(rt=.99),
+ Optim.Options(outer_iterations = 10,
+ iterations=1000,
+ show_trace=false, # turn it true to see the optimization
+ show_every=1));
+pμ=vcat(Optim.minimizer(ResF1));
+parf1=copy(par); parf1[2]=pμ[1]; μ1=pμ[2];
+
+## optimazation of β and order of incommensurate fractional order model
+function loss_F_8(pμ)
+ par[2] = pμ[1] # infectivity rate
+ μ = pμ[2:9] # order of derivatives
+ if size(X0,2) != Int64(ceil(maximum(μ))) # to prevent any errors regarding orders higher than 1
+ indx=findall(x-> x>1, μ)
+ μ[indx]=ones(length(indx))
+ end
+ _, x = FDEsolver(SIR, tspan, X0, μ, par, h = .1)
+ appX=vec(sum(x[1:10:end,[3,4,6]], dims=2))
+ rmsd(C, appX; normalize=:true)
+end;
+
+p_lo=vcat(1.4,.5*ones(8)); # lower bound for β and orders
+p_up=vcat(4,ones(8)); # upper bound for β and orders
+pvec=vcat(2.5,.9*ones(8)); # initial guess for β and orders
+ResF8=optimize(loss_F_8,p_lo,p_up,pvec,Fminbox(LBFGS()), # LBFGS is suitable for large scale problems
+# Result=optimize(loss,p_lo,p_up,pvec,SAMIN(rt=.99),
+ Optim.Options(outer_iterations = 10,
+ iterations=1000,
+ show_trace=false, # turn it true to see the optimization
+ show_every=1));
+pp=vcat(Optim.minimizer(ResF8));
+parf8=copy(par); parf8[2]=pp[1]; μ8=pp[2:9];
+
+## plotting
+DateTick=Date(2020,3,3):Day(1):Date(2020,5,17);
+DateTick2= Dates.format.(DateTick, "d u");
+
+t1, x1 = FDEsolver(SIR, tspan, X0, ones(8), par1, h = .1); # solve ode model
+_, xf1 = FDEsolver(SIR, tspan, X0, μ1*ones(8), parf1, h = .1); # solve commensurate fode model
+_, xf8 = FDEsolver(SIR, tspan, X0, μ8, parf8, h = .1); # solve incommensurate fode model
+
+X1=sum(x1[1:10:end,[3,4,6]], dims=2);
+Xf1=sum(xf1[1:10:end,[3,4,6]], dims=2);
+Xf8=sum(xf8[1:10:end,[3,4,6]], dims=2);
+
+Err1=rmsd(C, vec(X1)); # RMSD for ode model
+Errf1=rmsd(C, vec(Xf1)); # RMSD for commensurate fode model
+Errf8=rmsd(C, vec(Xf8)); # RMSD for incommensurate fode model
+
+plot(DateTick2,X1, ylabel="Daily new confirmed cases in Portugal",lw=5,
+ label="M1",xrotation=rad2deg(pi/3), linestyle=:dashdot)
+ plot!(Xf1, label="Mf1", lw=5)
+ plot!(Xf8, label="Mf8", linestyle=:dash, lw=5)
+ scatter!(C, label= "Real data",legendposition=(.85,1),legend=:false)
+ plPortugal=bar!(["M1" "Mf1" "Mf8"],[Err1 Errf1 Errf8], ylabel="Error (RMSD)",
+ legend=:false, bar_width=2,yguidefontsize=8,xtickfontsize=7,
+ inset = (bbox(0.04, 0.08, 70px, 60px, :right)),
+ subplot = 2,
+ bg_inside = nothing)
