Getting Started with Differential Equations in Julia

This tutorial will introduce you to the functionality for solving ODEs. Additionally, a video tutorial walks through this material.

Example 1 : Solving Scalar Equations

In this example, we will solve the equation

\[\frac{du}{dt} = f(u,p,t)\]

on the time interval $t\in[0,1]$ where $f(u,p,t)=αu$. Here, $u$ is the current state variable, $p$ is our parameter variable (containing things like a reaction rate or the constant of gravity), and $t$ is the current time.

(In our example, we know by calculus that the solution to this equation is $u(t)=u₀\exp(αt)$, but we will use DifferentialEquations.jl to solve this problem numerically, which is essential for problems where a symbolic solution is not known.)

The general workflow is to define a problem, solve the problem, and then analyze the solution. The full code for solving this problem is:

using DifferentialEquations
f(u, p, t) = 1.01 * u
u0 = 1 / 2
tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8)

using Plots
plot(sol, linewidth = 5, title = "Solution to the linear ODE with a thick line",
    xaxis = "Time (t)", yaxis = "u(t) (in μm)", label = "My Thick Line!") # legend=false
plot!(sol.t, t -> 0.5 * exp(1.01t), lw = 3, ls = :dash, label = "True Solution!")
Example block output

where the pieces are described below.

Step 1: Defining a Problem

To solve this numerically, we define a problem type by giving it the equation, the initial condition, and the timespan to solve over:

using DifferentialEquations
f(u, p, t) = 1.01 * u
u0 = 1 / 2
tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan)
ODEProblem with uType Float64 and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 0.5

Note that DifferentialEquations.jl will choose the types for the problem based on the types used to define the problem type. For our example, notice that u0 is a Float64, and therefore this will solve with the dependent variables being Float64. Since tspan = (0.0,1.0) is a tuple of Float64's, the independent variables will be solved using Float64's (note that the start time and end time must match types). You can use this to choose to solve with arbitrary precision numbers, unitful numbers, etc. Please see the notebook tutorials for more examples.

The problem types include many other features, including the ability to define mass matrices and hold callbacks for events. Each problem type has a page which details its constructor and the available fields. For ODEs, the appropriate page is here. In addition, a user can specify additional functions to be associated with the function in order to speed up the solvers. These are detailed at the performance overloads page.

Step 2: Solving a Problem

Controlling the Solvers

After defining a problem, you solve it using solve.

sol = solve(prob)
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 5-element Vector{Float64}:
 0.0
 0.09964258706516003
 0.3457024247583422
 0.6776921908052249
 1.0
u: 5-element Vector{Float64}:
 0.5
 0.552938681151017
 0.7089376245893467
 0.9913594502399238
 1.3728004409033037

The solvers can be controlled using the available options are described on the Common Solver Options manual page. For example, we can lower the relative tolerance (in order to get a more correct result, at the cost of more timesteps) by using the command reltol:

sol = solve(prob, reltol = 1e-6)
retcode: Success
Interpolation: specialized 7th order lazy interpolation, specialized 4rd order "free" stiffness-aware interpolation
t: 6-element Vector{Float64}:
 0.0
 0.08395921283331977
 0.29459485215346937
 0.5701065300378511
 0.8686054589334371
 1.0
u: 6-element Vector{Float64}:
 0.5
 0.5442490221301345
 0.6732716570549966
 0.889283154476355
 1.2021893261622514
 1.3728005075542231

There are many controls for handling outputs. For example, we can choose to have the solver save every 0.1 time points by setting saveat=0.1. Chaining this with the tolerance choice looks like:

sol = solve(prob, reltol = 1e-6, saveat = 0.1)
retcode: Success
Interpolation: 1st order linear
t: 11-element Vector{Float64}:
 0.0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1.0
u: 11-element Vector{Float64}:
 0.5
 0.553138320881742
 0.6119240040558909
 0.6769572322220972
 0.7489019734808645
 0.8284927602387266
 0.9165421885254866
 1.0139492143572488
 1.1217083317951
 1.240919726340709
 1.3728005075542231

More generally, saveat can be any collection of time points to save at. Note that this uses interpolations to keep the timestep unconstrained to speed up the solution. In addition, if we only care about the endpoint, we can turn off intermediate saving in general:

sol = solve(prob, reltol = 1e-6, save_everystep = false)
retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
 0.0
 1.0
u: 2-element Vector{Float64}:
 0.5
 1.3728005075542231

which will only save the final time point.

Choosing a Solver Algorithm

DifferentialEquations.jl has a method for choosing the default solver algorithm, which will find an efficient method to solve your problem. To help users receive the right algorithm, DifferentialEquations.jl offers a method for choosing algorithms through hints. This default chooser utilizes the precision of the number types and the keyword arguments (such as the tolerances) to select an algorithm. Additionally, one can provide alg_hints to help choose good defaults using properties of the problem and necessary features for the solution. For example, if we have a stiff problem where we need high accuracy, but don't know the best stiff algorithm for this problem, we can use:

sol = solve(prob, alg_hints = [:stiff], reltol = 1e-8, abstol = 1e-8)
retcode: Success
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 14-element Vector{Float64}:
 0.0
 0.012407826196308189
 0.04943554396660186
 0.09925889132467855
 0.15847071468568555
 0.2285503638356638
 0.30734723422488125
 0.3945147615835462
 0.48838261679111594
 0.5880567114707754
 0.6923232042782961
 0.800231609063282
 0.9109623745010172
 1.0
u: 14-element Vector{Float64}:
 0.5
 0.5063053789114713
 0.5255987021154772
 0.552724440835384
 0.5867879554585146
 0.62982623626939
 0.6819994183534772
 0.7447644602267912
 0.8188284026934183
 0.9055526107898983
 1.006117894091279
 1.1219707582414566
 1.25473553046225
 1.372800507495401

You can also explicitly choose the algorithm to use. DifferentialEquations.jl offers a much wider variety of solver algorithms than traditional differential equations libraries. Many of these algorithms are from recent research and have been shown to be more efficient than the “standard” algorithms. For example, we can choose a 5th order Tsitouras method:

sol = solve(prob, Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 5-element Vector{Float64}:
 0.0
 0.09964258706516003
 0.3457024247583422
 0.6776921908052249
 1.0
u: 5-element Vector{Float64}:
 0.5
 0.552938681151017
 0.7089376245893467
 0.9913594502399238
 1.3728004409033037

Note that the solver controls can be combined with the algorithm choice. Thus we can for example solve the problem using Tsit5() with a lower tolerance via:

sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8)
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 17-element Vector{Float64}:
 0.0
 0.012407826196308189
 0.042501278333560696
 0.0817804940926822
 0.12887385246570865
 0.18409797482126694
 0.24627458684331965
 0.31479299098506436
 0.38859636482100657
 0.46686178626184555
 0.5487161506239446
 0.633434723668386
 0.7203630271543068
 0.8089580249850934
 0.8987655385980757
 0.9894162242718207
 1.0
u: 17-element Vector{Float64}:
 0.5
 0.5063053789114713
 0.5219304750950854
 0.5430527156531716
 0.5695067808068426
 0.6021743690740426
 0.6412025714747711
 0.6871475333575261
 0.7403258498418478
 0.8012223528078949
 0.8702768771165198
 0.9480214886604926
 1.0350186821189897
 1.1319031558272872
 1.239373504325514
 1.3582039555461738
 1.3728005076225747

In DifferentialEquations.jl, some good “go-to” choices for ODEs are:

  • AutoTsit5(Rosenbrock23()) handles both stiff and non-stiff equations. This is a good algorithm to use if you know nothing about the equation.
  • AutoVern7(Rodas5()) handles both stiff and non-stiff equations in a way that's efficient for high accuracy.
  • Tsit5() for standard non-stiff. This is the first algorithm to try in most cases.
  • BS3() for fast low accuracy non-stiff.
  • Vern7() for high accuracy non-stiff.
  • Rodas4() or Rodas5() for small stiff equations with Julia-defined types, events, etc.
  • KenCarp4() or TRBDF2() for medium-sized (100-2000 ODEs) stiff equations
  • RadauIIA5() for really high accuracy stiff equations
  • QNDF() for large stiff equations

For a comprehensive list of the available algorithms and detailed recommendations, please see the solver documentation. Every problem type has an associated page detailing all the solvers associated with the problem.

Step 3: Analyzing the Solution

Handling the Solution Type

The result of solve is a solution object. We can access the 5th value of the solution with:

sol[5]
0.5695067808068426

or get the time of the 8th timestep by:

sol.t[8]
0.31479299098506436

Convenience features are also included. We can build an array using a comprehension over the solution tuples via:

[t + u for (u, t) in tuples(sol)]
17-element Vector{Float64}:
 0.5
 0.5187132051077795
 0.564431753428646
 0.6248332097458538
 0.6983806332725513
 0.7862723438953096
 0.8874771583180907
 1.0019405243425905
 1.1289222146628544
 1.2680841390697404
 1.4189930277404645
 1.5814562123288787
 1.7553817092732964
 1.9408611808123806
 2.13813904292359
 2.3476201798179943
 2.372800507622575

or more generally

[t + 2u for (u, t) in zip(sol.u, sol.t)]
17-element Vector{Float64}:
 1.0
 1.0250185840192507
 1.0863622285237315
 1.1678859253990255
 1.2678874140793939
 1.3884467129693523
 1.5286797297928618
 1.6890880577001166
 1.8692480645047023
 2.069306491877635
 2.289269904856984
 2.5294777009893714
 2.790400391392286
 3.0727643366396675
 3.377512547249104
 3.7058241353641685
 3.7456010152451493

allows one to use more parts of the solution type. The object that is returned by default acts as a continuous solution via an interpolation. We can access the interpolated values by treating sol as a function, for example:

sol(0.45) # The value of the solution at t=0.45
0.7876927465687831

Note the difference between these: indexing with [i] is the value at the ith step, while (t) is an interpolation at time t!

If in the solver dense=true (this is the default unless saveat is used), then this interpolation is a high order interpolation and thus usually matches the error of the solution time points. The interpolations associated with each solver is detailed at the solver algorithm page. If dense=false (unless specifically set, this only occurs when save_everystep=false or saveat is used) then this defaults to giving a linear interpolation.

For more details on handling the output, see the solution handling page.

Plotting Solutions

While one can directly plot solution time points using the tools given above, convenience commands are defined by recipes for Plots.jl. To plot the solution object, simply call plot:

#]add Plots # You need to install Plots.jl before your first time using it!
using Plots
#plotly() # You can optionally choose a plotting backend
plot(sol)
Example block output

The plot function can be formatted using the attributes available in Plots.jl. Additional DiffEq-specific controls are documented at the plotting page.

For example, from the Plots.jl attribute page, we see that the line width can be set via the argument linewidth. Additionally, a title can be set with title. Thus we add these to our plot command to get the correct output, fix up some axis labels, and change the legend (note we can disable the legend with legend=false) to get a nice-looking plot:

plot(sol, linewidth = 5, title = "Solution to the linear ODE with a thick line",
    xaxis = "Time (t)", yaxis = "u(t) (in μm)", label = "My Thick Line!") # legend=false
Example block output

We can then add to the plot using the plot! command:

plot!(sol.t, t -> 0.5 * exp(1.01t), lw = 3, ls = :dash, label = "True Solution!")
Example block output

Example 2: Solving Systems of Equations

In this example, we will solve the Lorenz equations:

\[\begin{aligned} \frac{dx}{dt} &= σ(y-x) \\ \frac{dy}{dt} &= x(ρ-z) - y \\ \frac{dz}{dt} &= xy - βz \\ \end{aligned}\]

Defining your ODE function to be in-place updating can have performance benefits. What this means is that, instead of writing a function which outputs its solution, you write a function which updates a vector that is designated to hold the solution. By doing this, DifferentialEquations.jl's solver packages are able to reduce the amount of array allocations and achieve better performance.

The way we do this is we simply write the output to the 1st input of the function. For example, our Lorenz equation problem would be defined by the function:

function lorenz!(du, u, p, t)
    du[1] = 10.0 * (u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end
lorenz! (generic function with 1 method)

and then we can use this function in a problem:

using DifferentialEquations
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz!, u0, tspan)
sol = solve(prob)
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 1263-element Vector{Float64}:
   0.0
   3.5678604836301404e-5
   0.0003924646531993154
   0.0032624077544510573
   0.009058075635317072
   0.01695646895607931
   0.02768995855685593
   0.04185635042021763
   0.06024041165841079
   0.08368541255159562
   ⋮
  99.30760258626904
  99.39665422328268
  99.49536147459878
  99.58822928767293
  99.68983993598462
  99.77864535713971
  99.85744078539504
  99.93773320913628
 100.0
u: 1263-element Vector{Vector{Float64}}:
 [1.0, 0.0, 0.0]
 [0.9996434557625105, 0.0009988049817849058, 1.781434788799208e-8]
 [0.9961045497425811, 0.010965399721242457, 2.146955365838907e-6]
 [0.9693591634199452, 0.08977060667778931, 0.0001438018342266937]
 [0.9242043615038835, 0.24228912482984957, 0.0010461623302512404]
 [0.8800455868998046, 0.43873645009348244, 0.0034242593451028745]
 [0.8483309847495312, 0.6915629321083602, 0.008487624590227805]
 [0.8495036669651213, 1.0145426355349096, 0.01821208962127994]
 [0.9139069574560097, 1.4425599806525806, 0.03669382197085303]
 [1.088863826836895, 2.052326595543049, 0.0740257368585531]
 ⋮
 [4.669609096878053, 3.061564434452441, 25.1424735017959]
 [4.188801916573263, 4.617474401440693, 21.09864175382292]
 [5.559603854699961, 7.905631612648314, 18.79323210016923]
 [8.556629716266505, 12.533041060088328, 20.6623639692711]
 [12.280585075547771, 14.505154761545633, 29.332088452699942]
 [11.736883151600804, 8.279294641640229, 34.68007510231878]
 [8.10973327066804, 3.2495066495235854, 31.97052076740117]
 [4.958629886040755, 2.194919965065022, 26.948439650907677]
 [3.8020065515435855, 2.787021797920187, 23.420567509786622]

Using the plot recipe tools defined on the plotting page, we can choose to do a 3D phase space plot between the different variables:

using Plots
plot(sol, idxs = (1, 2, 3))
Example block output

Note that the default plot for multidimensional systems is an overlay of each timeseries. We can plot the timeseries of just the second component using the variable choices interface once more:

plot(sol, idxs = (0, 2))
Example block output

Note that here “variable 0” corresponds to the independent variable (“time”).

Defining Parameterized Functions

Often, you want to explicitly have parameters associated with your differential equations. This can be used by things like parameter estimation routines. In this case, you use the p values via the syntax:

function parameterized_lorenz!(du, u, p, t)
    du[1] = p[1] * (u[2] - u[1])
    du[2] = u[1] * (p[2] - u[3]) - u[2]
    du[3] = u[1] * u[2] - p[3] * u[3]
end
parameterized_lorenz! (generic function with 1 method)

and then we add the parameters to the ODEProblem:

u0 = [1.0, 0.0, 0.0]
tspan = (0.0, 1.0)
p = [10.0, 28.0, 8 / 3]
prob = ODEProblem(parameterized_lorenz!, u0, tspan, p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 3-element Vector{Float64}:
 1.0
 0.0
 0.0

We can make our functions look nicer by doing a few tricks. For example:

function parameterized_lorenz!(du, u, p, t)
    x, y, z = u
    σ, ρ, β = p
    du[1] = dx = σ * (y - x)
    du[2] = dy = x * (ρ - z) - y
    du[3] = dz = x * y - β * z
end
parameterized_lorenz! (generic function with 1 method)

Note that the type for the parameters p can be anything: you can use arrays, static arrays, named tuples, etc. to enclose your parameters in a way that is sensible for your problem.

Since the parameters exist within the function, functions defined in this manner can also be used for sensitivity analysis, parameter estimation routines, and bifurcation plotting. This makes DifferentialEquations.jl a full-stop solution for differential equation analysis which also achieves high performance.

Example 3: Solving Nonhomogeneous Equations using Parameterized Functions

Parameterized functions can also be used for building nonhomogeneous ordinary differential equations (these are also referred to as ODEs with nonzero right-hand sides). They are frequently used as models for dynamical systems with external (in general time-varying) inputs. As an example, consider a model of a pendulum consisting of a slender rod of length l and mass m:

\[\begin{aligned} \frac{\mathrm{d}\theta(t)}{\mathrm{d}t} &= \omega(t)\\ \frac{\mathrm{d}\omega(t)}{\mathrm{d}t} &= - \frac{3}{2}\frac{g}{l}\sin\theta(t) + \frac{3}{ml^2}M(t) \end{aligned},\]

where θ and ω are the angular deviation of the pendulum from the vertical (hanging) orientation and the angular rate, respectively, M is an external torque (developed, say, by a wind or a motor), and finally, g stands for gravitational acceleration.

using DifferentialEquations
using Plots

l = 1.0                             # length [m]
m = 1.0                             # mass [kg]
g = 9.81                            # gravitational acceleration [m/s²]

function pendulum!(du, u, p, t)
    du[1] = u[2]                    # θ'(t) = ω(t)
    du[2] = -3g / (2l) * sin(u[1]) + 3 / (m * l^2) * p(t) # ω'(t) = -3g/(2l) sin θ(t) + 3/(ml^2)M(t)
end

θ₀ = 0.01                           # initial angular deflection [rad]
ω₀ = 0.0                            # initial angular velocity [rad/s]
u₀ = [θ₀, ω₀]                       # initial state vector
tspan = (0.0, 10.0)                  # time interval

M = t -> 0.1sin(t)                    # external torque [Nm]

prob = ODEProblem(pendulum!, u₀, tspan, M)
sol = solve(prob)

plot(sol, linewidth = 2, xaxis = "t", label = ["θ [rad]" "ω [rad/s]"], layout = (2, 1))
Example block output

Note how the external time-varying torque M is introduced as a parameter in the pendulum! function. Indeed, as a general principle the parameters can be any type; here we specify M as time-varying by representing it by a function, which is expressed by appending the dependence on time (t) to the name of the parameter.

Note also that, in contrast with the time-varying parameter, the (vector of) state variables u, which is generally also time-varying, is always used without the explicit dependence on time (t).

Example 4: Using Other Types for Systems of Equations

DifferentialEquations.jl can handle many different dependent variable types (generally, anything with a linear index should work!). So instead of solving a vector equation, let's let u be a matrix! To do this, we simply need to have u0 be a matrix, and define f such that it takes in a matrix and outputs a matrix. We can define a matrix of linear ODEs as follows:

using DifferentialEquations
using Plots
A = [1.0 0 0 -5
    4 -2 4 -3
    -4 0 0 1
    5 -2 2 3]
u0 = rand(4, 2)
tspan = (0.0, 1.0)
f(u, p, t) = A * u
prob = ODEProblem(f, u0, tspan)
ODEProblem with uType Matrix{Float64} and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 4×2 Matrix{Float64}:
 0.4854    0.983144
 0.481319  0.362922
 0.279091  0.0245504
 0.334704  0.829293

Here our ODE is on a 4x2 matrix, and the ODE is the linear system defined by multiplication by A. To solve the ODE, we do the same steps as before.

sol = solve(prob)
plot(sol)
Example block output

We can instead use the in-place form by using Julia's in-place matrix multiplication function mul!:

using LinearAlgebra
f(du, u, p, t) = mul!(du, A, u)
f (generic function with 2 methods)

Additionally, we can use non-traditional array types as well. For example, StaticArrays.jl offers immutable arrays which are stack-allocated, meaning that their usage does not require any (slow) heap-allocations that arrays normally have. This means that they can be used to solve the same problem as above, with the only change being the type for the initial condition and constants:

using StaticArrays
A = @SMatrix [1.0 0.0 0.0 -5.0
    4.0 -2.0 4.0 -3.0
    -4.0 0.0 0.0 1.0
    5.0 -2.0 2.0 3.0]
u0 = @SMatrix rand(4, 2)
tspan = (0.0, 1.0)
f2(u, p, t) = A * u
prob = ODEProblem(f2, u0, tspan)
sol = solve(prob)
plot(sol)
Example block output

Note that the analysis tools generalize over to systems of equations as well.

sol[4]
4×2 StaticArraysCore.SMatrix{4, 2, Float64, 8} with indices SOneTo(4)×SOneTo(2):
 0.24258   -0.0591023
 0.646122   0.594281
 0.361672   0.342085
 0.933108   1.56089

still returns the solution at the fourth timestep. It also indexes into the array as well. The last value is the timestep, and the beginning values are for the component. This means

sol[5, 3]
0.5133782877492716

is the value of the 5th component (by linear indexing) at the 3rd timepoint, or

sol[2, 1, :]
10-element Vector{Float64}:
  0.3611711044854693
  0.5454537040676509
  0.6912996674706529
  0.6461221743527011
  0.35879138615607237
 -0.17814461785658348
 -0.7722102486231299
 -0.8741980923154496
  0.3637358592202831
  1.983571313217542

is the timeseries for the component, which is the 2nd row and 1 column.

Going Beyond ODEs: How to Use the Documentation

Not everything can be covered in the tutorials. Instead, this tutorial will end by pointing you in the directions for the next steps.

Common API for Defining, Solving, and Plotting

One feature of DifferentialEquations.jl is that this pattern for solving equations is conserved across the different types of differential equations. Every equation has a problem type, a solution type, and the same solution handling (+ plotting) setup. Thus the solver and plotting commands in the Basics section applies to all sorts of equations, like stochastic differential equations and delay differential equations. Each of these different problem types are defined in the Problem Types section of the docs. Every associated solver algorithm is detailed in the Solver Algorithms section, sorted by problem type. The same steps for ODEs can then be used for the analysis of the solution.

Additional Features and Analysis Tools

In many cases, the common workflow only starts with solving the differential equation. Many common setups have built-in solutions in DifferentialEquations.jl. For example, check out the features for:

Many more are defined in the relevant sections of the docs. Please explore the rest of the documentation, including tutorials for getting started with other types of equations. In addition, to get help, please either file an issue at the main repository or come have an informal discussion at the Julia Zulip chatroom.