DifferentialEquations.jl is a package for solving differential equations in Julia. It covers discrete equations (function maps, discrete stochastic (Gillespie/Markov) simulations), ordinary differential equations, stochastic differential equations, algebraic differential equations, delay differential equations, hybrid differential equations, jump diffusions, and (stochastic) partial differential equations. Through extensive use of multiple dispatch, metaprogramming, plot recipes, foreign function interfaces (FFI), and call-overloading, DifferentialEquations.jl offers a unified user interface to solve and analyze various forms of differential equations while not sacrificing features or performance. Many modern features are integrated into the solvers, such as allowing arbitrary user-defined number systems for high-precision and arithmetic with physical units, built-in multithreading and parallelism, and symbolic calculation of Jacobians. Integrated into the package is an algorithm testing and benchmarking suite to both ensure accuracy and serve as an easy way for researchers to develop and distribute their own methods. Together, these features build a highly extendable suite which is feature-rich and highly performant.

Differential equations are fundamental components of many scientific models; they are used to describe large-scale physical phenomena like planetary systems [

However, these software packages contain many limitations which stem from their implementation and the time when they were developed. Since the time of their inception, many other forms of differential equations have become commonplace tools not only for mathematicians, but throughout the sciences. Stochastic differential equations (SDEs), have become more commonplace not only in mathematical finance [

Other design limitations stem from the programming languages used in the implementation. Many of these algorithms, being developed in early C/Fortran, do not have abstractions for generalized array formats. In order to use these algorithms, one must provide the solver with a vector. In cases where a matrix or a higher dimensional tensor are the natural representation of the differential equation, the user is required to transform their equation into a vector equation for use in these solvers. Also, these solvers are limited to using 64-bit floating point calculations. The numerical precision limits their use in high-precision applications, requiring specialized codes when precision lower than 10^{–16} is required. Lastly, many times these programs are interfaced via a scripting language where looping is not optimized and where “vectorized” codes provide the most efficient solution. However, vectorized coding in the style of MATLAB or NumPy results in temporary allocations and can lack compiler optimizations which require type inference. This increases the computational burden of the user-defined functions which degrades the efficiency of the solver.

The goal of DifferentialEquations.jl is build off of the foundation created by these previous differential equation libraries and modernize them using Julia. Julia is a scripting language, used in-place of languages like R, Python, MATLAB, but offers the performance one would associate with low-level compiled languages. This allows users to start prototypes in Julia, but also solve their large-scale models within the same language, instead of resorting to two language solutions when performance is needed. The language achieves this goal by extensive utilization of multiple dispatch and metaprogramming to design a language that is both easy for a compiler to understand and easy for a programmer to use [

We start by describing the innovations in usability. In Section 1.1 we show how multiple dispatch is used to consolidate the functions the user needs to into simple descriptive commands like

DifferentialEquations.jl uses multiple dispatch on specialized types to arrive at a unified user-API for the different types of equations. To use the package, one follows the steps:

Define a problem.

Solve the problem.

Plot the solution.

This standardization of the API makes complicated solvers accessible to less programming-inclined individuals, giving a good framework for future development and allows for the latest research in numerical differential equations to be utilized without complications.

To define a problem, a user must call the constructor for the appropriate problem object. Since ordinary differential equations (ODEs) are represented in the general form as

the _{0}

Many other examples are provided in the documentation

By using a dispatch architecture on

For most other packages, one would normally have to define

The solver returns a solution object which holds all of the information about the solution. Dispatches to array functions are provided on the

The solution can be plotted using the provided plot recipes through Plots.jl

These defaults are deliberately made so that a standard user does not need to dig further into the manual and understand the differences between all of the algorithms. However, an extensive set of functionality is available if the user wishes. All of these functions can be modified via additional arguments. For example, to change the solver algorithm to a highly efficient Order 7 method due to Verner [

The output of this command is shown in Figure

Example of the ODE plot recipe. This plot was created using the PyPlot backend through Plots.jl. Shown is the solution to the 4 × 2 ODE with

Lastly, these solvers tie into Julia integrated development environments (IDEs) to further enhance the ease of use. Users of the Juno IDE [

By using multiple-dispatch, the same user API is offered for other types of equations. For example, if one wishes to solve a stochastic differential equation (SDE):

then one builds an SDEProblem object by specifying the initial condition and now the two functions,

While this user interface is simple, the default methods these algorithms can call are efficient high-order solvers with adaptive timestepping [

Again, the same user API is offered for the available stochastic PDE solvers. Instead, one builds a

Additional keyword arguments can be supplied to

Most differential equations packages require that the user understands some details about the implementation of the library. However, the DifferentialEquations.jl ecosystem implements various Domain-Specific Languages (DSLs) via macros in order to give more natural options for defining mathematical constructs. In this section we will demonstrate the DSL for defining ODEs. For demonstrations related to other types of equations, please see the documentation.

The famous Lorenz system is mathematically defined as

A user must re-write this function in a “computer friendly format”, defining

While this format is accepted by DifferentialEquations.jl, additional usability macros are provided which will automatically translate user input from a more mathematical format. For ODEs,

Since Julia allows for the use of Unicode within code, this format matches the style one would expect to see in a TeX’d publication. The macro takes in this definition, finds the values for the left-hand side of the form “d__”, and uses a dictionary in order to find/replace these values to write a function which is in the format of the other scripting language libraries. Thus the translation to a vector system can be done by DifferentialEquations.jl, allowing the users to have more readable scripts while not sacrificing performance. In addition, the macro produces a function which updates an input

A unique feature from this form of function definition is that the parameters are built into the function type itself. The actual implementation involves creating a type

Also, since the code is analyzed by the program at the expression level, silent optimizations are able to be performed. For example, during the construction of the function, the code is transformed into a symbolic form for use in the high-performance CAS SymEngine.jl

The full functionality of DifferentialEquations.jl is defined in the more than 40 packages in the JuliaDiffEq Github organization

DifferentialEquations.jl is designed to use multiple-dispatch in order to allow for the different solver methods to be defined in separate packages. In Julia, the command:

calls a different “method” depending on the type of

OrdinaryDiffEq.jl

Sundials.jl

ODE.jl

ODEInterface.jl

LSODA.jl

Some of the packages, like OrdinaryDiffEq.jl and ODE.jl, are native Julia libraries, whereas Sundials.jl, ODEInterface.jl, and LSODA.jl are all interfaces to popular C and Fortran codes. Through this interface, users can switch between libraries by switching only the algorithm choice, and new packages can extend what’s available as needed. Note that the information in this structure is unidirectional: anyone can add a new package of solvers to the

The add-on packages are a set of functionality which use the common solve interface. For example, parameter estimation functionality is provided by defining algorithms which only use the abstraction of the

Julia’s base library defines its standard numeric types, _{0} (the initial condition) and timespan (i.e. separate types are allowed for the dependent and independent variables). It then calls a type-dependent integration function which is optimized via JIT compilation for the numeric types given to the function. Different dispatches are given for subtypes of Number and AbstractArray since arrays are mutable and heap allocated, meaning that when numbers are treated directly instead of as arrays of size 1 a large speedup can occur. This allows for the internal integration algorithms to achieve C/Fortran speeds, while allowing for the generic numerical types and the readability of being in Julia itself. The following subsections highlight two important examples.

One advantage of this design beyond speed is that it allows the user of DifferentialEquations.jl to use any type which is a subclass of Number as the number type for the equations. This includes not only the basic types like

The combination of high-performance number systems with high order Runge-Kutta methods such as the Order 14 methods due to Feagin allows for fast solving with high accuracy. For an example showing this combination, see the “Feagin’s Order 10, 12, and 14 methods” notebook in the examples folder

This design also allows DifferentialEquations.jl to be compatible with number systems which have physical units. SIUnits.jl

The attentive reader should realize that this will correctly throw an error: the output of the function in an ODE must be a rate, and therefore must have units of

This will produce an output whose units are in terms of Newtons, and with time in terms of seconds.

DifferentialEquations.jl also includes parallelism whenever possible. One area where parallelism is currently being employed is via “within-method” parallelism for Runge-Kutta methods. Using Julia’s experimental multithreading, DifferentialEquations.jl provides a multithreaded version of the

Also, DifferentialEquations.jl provides methods for performing parallel Monte Carlo simulations. Using Julia’s pmap construct, one is able to specify for a problem to be solved N times, and DifferentialEquations.jl will distribute this automatically across multiple nodes of a cluster. A vector of results along with summary statistics is returned for the solution. This functionality has been tested on the local UC Irvine cluster (using SGE) and the XSEDE Comet cluster (using Slurm).

DifferentialEquations.jl includes a suite specifically designed for researchers interested in developing new methods for differential equations (like the authors themselves). This includes functionality for easy integration of new methods, extensive testing, and a benchmarking suite.

The design of DifferentialEquations.jl allows for users to add new integration methods by adding new dispatches. One way to add new methods is to simply create a new package which extends the _{n}_{n+1}

Additionally a new

The DifferentialEquations.jl suite includes a large number of testing functions to ensure correctness of all of the algorithms. Many premade problems with analytical solutions are provided and convergence testing functionality is included to be able to test the order of accuracy and plot the results. All of the DifferentialEquations.jl algorithms are tested using the Travis and AppVoyer Continuous Integration (CI) testing services to ensure correctness.

Lastly, a benchmarking suite is included to test the efficiency of different algorithms. Two forms of benchmarking are included: the

As of this publication, the benchmarks show that native methods from OrdinaryDiffEq.jl (which were developed in tandem with DifferentialEquations.jl) achieve an order of magnitude speedup on nonstiff problems when achieving the same error over the classic Hairer Runge-Kutta implementations and the ODE.jl implementations.

While DifferentialEquations.jl already offers many new features and high performance, the package is still under heavy development and will be for the foreseeable future. Currently, most of the methods for stiff equations are wrapped methods (only the Rosenbrock with/without local extrapolation, and the order 1/2 BDF methods exist in native forms). While these methods, such as CVODE (provided by Sundials.jl

The native Julia methods have far less limitations because they work on the general abstract types

However, while the genericness of the implementation makes it very flexible, one limitation of the design is that the full extent of the compatibility is not able to be easily documented or known. The practice of “duck typing” means that the generic functions are left open ended, and functionality will work on available types depending on whether certain traits or operations are defined. For example, Julia-defined numbers systems are compatible with the if certain operations (+,–,*,/) are defined. However, different solver algorithms can have slightly different compatibility requirements. Adaptive timestepping also requires that the number system has a well-defined

Also, DifferentialEquations.jl is currently limited on the types of PDEs it natively supports, and the mesh generation tools are still in their infancy. To address these issues and more, planned functionality includes (but is not limited to):

Finite Difference Methods for common elliptic, parabolic, and hyperbolic PDEs, including high order methods for SPDEs

Highly parallel accelerated solvers using GPGPUs and Xeon Phi cards (prototypes have already been developed [

High order methods for stiff SDEs

Check the repository issues for the most up to date roadmap.

Continuous Integration testing with the latest versions of Julia on Mac, Linux, and Windows are provided via Travis and AppVoyer. These tests check most of the features of DifferentialEquations.jl, including the convergence of each algorithm, the ability to plot, the number types used in the computations, and more. Coveralls and Coverage badges are provided on the repository for test coverage analysis. As with other Julia packages, a user can check to see if these functionalities are working on their local machine via the command Pkg.test (“DifferentialEquations”). Benchmarks in Jupyter notebooks are provided to test the differences between the integrator implementations. Additionally, each Base library, component solver, and add-on package contains its own set of tests for the functionality that it implements. All have the same continuous integration setup.

Another method for quality control is user feedback. DifferentialEquations.jl receives bug reports and feature requests through the Julialang Discourse

DifferentialEquations.jl is CI tested on MacOSX and Linux via Travis CI, and Windows via AppVoyer.

Julia v0.5+

Dependencies are split into two groups. The direct dependencies of DifferentialEquations.jl are the packages of JuliaDiffEq which are built around the common interface and developed in tandem with DifferentialEquations.jl. The indirect dependencies are the dependencies of the direct dependencies, which are packages which are not actively developed as part of JuliaDiffEq activity.

The direct dependencies of DifferentialEquations.jl are the packages of JuliaDiffEq. These are:

DiffEqBase.jl

StochasticDiffEq.jl

FiniteElementDiffEq.jl

DiffEqDevTools.jl

OrdinaryDiffEq.jl

AlgebraicDiffEq.jl

StokesDiffEq.jl

DiffEqParamEstim.jl

DiffEqSensitivity.jl

Sundials.jl

ODEInterfaceDiffEq.jl

ParameterizedFunctions.jl

DiffEqPDEBase.jl

DelayDiffEq.jl

DiffEqCallbacks.jl

DiffEqMonteCarlo.jl

DiffEqJump.jl

DiffEqFinancial.jl

DiffEqBiological.jl

MultiScaleArrays.jl

Indirect dependencies include:

RecipesBase.jl

Optim.jl

Parameters.jl

ForwardDiff.jl

IterativeSolvers.jl

GenericSVD.jl

Compat.jl

InplaceOps.jl

SymEngine.jl

All of these dependencies will automatically install upon

Optional dependencies of DifferentialEquations.jl include the additional solver packages:

ODEInterface.jl

ODE.jl

LSODA.jl

For information on how to install these libraries, see their respective repositories.

Christopher Rackauckas, Lead Developer of JuliaDiffEq

Mauro Werder, contributions to DiffEqBase.jl

Scott P. Jones, contributions to DiffEqBase.jl

Virgile Andreani, contributed the enhanced plotting functionality

Ethan Levien, contributions to DiffEqMonteCarlo.jl and DiffEqJump.jl

Michael Fiano, contributions to the documentation

David Barton, contributions to the dense output in OrdinaryDiffEq.jl

English

Differential equations form the bedrock of many scientific fields. Therefore, there is no question as to whether numerical differential equation solvers will be used, rather the question is which ones will be used. Julia is a relatively young language which is seeing rapid adoption in the fields of data science and scientific computing due to the performance and productivity that it offers. Because of this, many scientists using Julia will need these tools either as a means to analyze models themselves, or as intermediate tools in more complex methods. With its applicability to many classes of differential equations, its included analysis tools for performing parameter estimation and sensitivity analysis, and the rapid pace at which this software is being developed, DifferentialEquations.jl looks to be a viable choice for many Julians looking for a differential equations library. Lastly, as an open-source software with a modular structure, it is easily extendable. For information on how to extend the functionality of DifferentialEquations.jl please see the Contributor’s Guide at DiffEqDevDocs.jl.

We would like to thank the Julia community for the support they have given me throughout this project. Special thanks goes out to Ismael Venegas Castello (@Ismael-VC), Lyndon White (@oxinabox), Fengyang Wang (@TotalVerb), and Scott P. Jones (@ScottPJones) from whom CR learned many languages tricks in our long Gitter discussions. CR would also like to acknowledge Tom Breloff (@tbreloff) for his work on Plots.jl and the development of the recipe concept which allowed DifferentialEquations.jl to so easily meld with the plot functionalities. Lastly, CR would like to thank @finmod for repeatedly helping to identify documentation which needs updates.

The authors have no competing interests to declare.