FdeSolver: A Julia Package for Solving Fractional Differential Equations

TL;DR

FdeSolver addresses solving Caputo fractional differential equations with memory by providing a Julia-based, open-source solver that combines product-integration discretization with predictor-corrector and Newton–Raphson schemes, accelerated by FFT. It presents the Volterra integral form $\mathbf{X}(t)=T_{m-1}[\mathbf{X};t_0]+\frac{1}{\Gamma(\boldsymbol{\beta})}\int_{t_0}^t (t-\tau)^{\boldsymbol{\beta}-1}\mathbf{F}(\tau,\mathbf{X}(\tau))d\tau$ and describes the software architecture, reproducible development practices, and benchmarking against Matlab and FractionalDiffEq across one- and multi-dimensional problems. It demonstrates applications to microbial community dynamics and Covid-19 transmission modeling with both commensurate and incommensurate orders, showing improved speed and accuracy and illustrating memory effects as a modeling feature. The code and data are openly available on GitHub and Zenodo, enabling reuse, verification, and extension within the Julia ecosystem.

Abstract

Implementing and executing numerical algorithms to solve fractional differential equations has been less straightforward than using their integer-order counterparts, posing challenges for practitioners who wish to incorporate fractional calculus in applied case studies. Hence, we created an open-source Julia package, FdeSolver, that provides numerical solutions for fractional-order differential equations based on product-integration rules, predictor-corrector algorithms, and the Newton-Raphson method. The package covers solutions for one-dimensional equations with orders of positive real numbers. For high-dimensional systems, the orders of positive real numbers are limited to less than (and equal to) one. Incommensurate derivatives are allowed and defined in the Caputo sense. Here, we summarize the implementation for a representative class of problems, provide comparisons with available alternatives in Julia and Matlab, describe our adherence to good practices in open research software development, and demonstrate the practical performance of the methods in two applications; we show how to simulate microbial community dynamics and model the spread of Covid-19 by fitting the order of derivatives based on epidemiological observations. Overall, these results highlight the efficiency, reliability, and practicality of the FdeSolver Julia package.

Figures (7)

• Figure 1: A schematic of FdeSolver package v1.0.7. This package can solve one-dimensional fractional models for the order derivative $\beta>0$ and multi-dimensional ones for the order derivatives $0<\beta \leq 1$ of a class of problem \ref{['eq:scalar']}. The name of the solver function is FDEsolver. There are ten arguments where the four positional arguments are F: the equation function, tSpan: the domain of the problem, $\beta$: the orders of derivatives, and X0: the initial conditions and the four optional arguments are h: step-size of computation, nc: the number of correction for the predictor-corrector method, itmax: the maximum number of iterations for the Newton-Raphson method, and tol: the tolerance of each iteration or correction. The input par is any additional parameters related to the problem and JF is a Jacobian function needed for the Netwon-Raphson method. The predictor-corrector method \ref{['Sec: pc']} is encoded in main.jl and SubFuns.jl, and the Newton-Raphson method \ref{['Sec: NR']} in main_Jacob.jl and SubFuns_Jacob.jl. Structs.jl leverages multiple dispatch to choose a proper numerical algorithm based on the presence or the lack of a Jacobian function and tune its functionality based on arguments' types. The FFTW, LinearAlgebra, and SpecialFunctions Julia packages support the FdeSolver package.
• Figure 2: Error versus execution time of FdeSolver and Matlab codes for solving the one-dimensional differential equation of (a) Example \ref{['sec:nonstiff']}, (b) Example \ref{['sec: stiff']}, and (c) Example \ref{['sec: high-order']}. The axes are visualized on a logarithmic scale. The unit of the X-axes is second (Sec) and the Y-axes show the Euclidean norm of the difference between the approximations and exact values. The notations used for the methods of FdeSolver are; J1: the PC method, and J2: the NR method. The notations used for the methods of the Matlab routines are; M1: the PC method, M2: the NR method, M3: the NR method but with the PI rectangular rule, and M4: the explicit PI rectangular rule without PC.
• Figure 3: Error versus execution time of Julia and Matlab codes for solving for multi-dimensional systems of (a) Example \ref{['Sec: SIR']} and (b) Example \ref{['Sec:LVmulti']}. The axes are visualized on a logarithmic scale. The unit of the X-axes is second (Sec) and the Y-axes show the accuracy by measuring the Euclidean norm of the difference between the approximations and the results with fine step size $2^{-10}$. The notations used for the methods of FdeSolver are; J1: the PC method, and J2: the NR method. The notations used for the methods of FractionalDiffEq are; J3: the PC method, J4: the explicit PI rectangular rule without PC, J5: according to the FOTF Toolbox, J6, J7, and J8: the fractional linear multistep methods. The notations used for the methods of the Matlab routines are; M1: the PC method, M2: the NR method, M3: the NR method but with the PI rectangular rule, and M4: the explicit PI rectangular rule without PC.
• Figure 4: This is the replication of Fig. 2 of Ref. KhalighiPlosCB using the mature FdeSolver Julia package to illustrate the influence of commensurate fractional derivatives on resistance and resilience of the 3-species community model \ref{['Eq: microb']}. Let us suppose $X_B$, $X_R$, and $X_G$ as the abundance of blue, red, and gray species, where indexes $B, R$, and $G$ indicate the colors and present through all parameters. We solve the model with initial conditions $X_B(0)=0.99,$$X_R(0)=0.01,$ and $X_B(0)=0.01$, and parameters $K_{i,j}=0.1, \forall{i \neq j}$, $n=2$, $k_i=1$. The order of the derivatives is set to integer one and fractional values $\beta_B=\beta_R=\beta_G=0.9$. (a) The growth rates of the system during relaxation are $b_B=1$, $b_R=0.95$, and $b_G=1.05$. A pulse perturbation is applied to the system, during a time window indicated via the grey background, by lowering the growth rate of the blue species ($b_B=0.5$) and raising that of the gray species ($b_G=2$). (b) The perturbation temporarily moves away the system from the original stable state. However, fractional orders increase resistance to perturbation since the community (solid lines) is not displaced as far from its initial state compared to the dynamics of the system with integer order (dashed lines). (c) A slightly stronger pulse perturbation is applied to the growth rate of the gray species ($b_G=2.2$). (d) It triggers a shift toward an alternative stable state dominated by the gray species in the system with integer order (dashed lines). However, fractional orders can also entirely prevent a state shift (solid lines) by increasing both resistance and resilience to perturbation.
• Figure 5: The comparison between real data on the daily new confirmed cases as retrieved from CSSE Data-Covid and the estimation of $I+P+H$ from the system equation \ref{['eq: Covid']} with integer order derivatives (model M1 shown by the dash-dotted blue line), commensurate fractional derivatives (model Mf1 shown by the solid orange line), and incommensurate fractional derivatives (model Mf8 shown by the green dashed line). The circles indicate the real data for (a) Spain from 25 February to 7 June and (b) Portugal from 3 March to 17 May. The infection rate $\beta$ is fitted for the M1 model, the commensurate order for the Mf1 model together with $\beta$ are optimized, and finally, the eight incommensurate orders and $\beta$ for model Mf8 are optimized to fit real data. The specified error bars are based on root mean square deviation (RMSD) such that the fewer values indicate the less residual variance. The fitted values are listed in Table \ref{['tab:fitValues']}.