Type handling of DiffEq.jl

Chris Rackauckas, Mosè Giordano

Solving Equations in With Julia-Defined Types

One of the nice things about DifferentialEquations.jl is that it is designed with Julia's type system in mind. What this means is, if you have properly defined a Number type, you can use this number type in DifferentialEquations.jl's algorithms!

[Note that this is restricted to the native algorithms of OrdinaryDiffEq.jl. The other solvers such as ODE.jl, Sundials.jl, and ODEInterface.jl are not compatible with some number systems.]

DifferentialEquations.jl determines the numbers to use in its solvers via the types that are designated by tspan and the initial condition of the problem. It will keep the time values in the same type as tspan, and the solution values in the same type as the initial condition. [Note that adaptive timestepping requires that the time type is compaible with sqrt and ^ functions. Thus dt cannot be Integer or numbers like that if adaptive timestepping is chosen].

Let's solve the linear ODE first define an easy way to get ODEProblems for the linear ODE:

using DifferentialEquations

f = (u,p,t) -> (p*u)

prob_ode_linear = ODEProblem(f,1/2,(0.0,1.0),1.01);

First let's solve it using Float64s. To do so, we just need to set u0 to a Float64 (which is done by the default) and dt should be a float as well.

prob = prob_ode_linear

sol = solve(prob,Tsit5())

Notice that both the times and the solutions were saved as Float64. Let's change the time to use rational values. Rationals are not compatible with adaptive time stepping since they do not have an L2 norm (this can be worked around by defining internalnorm, but rationals already explode in size!). To account for this, let's turn off adaptivity as well:

prob = ODEProblem(f, 1/2, (0//1,1//1), 101//100);

sol = solve(prob, RK4(), dt=1//2^(6), adaptive=false)

println(sol)

Now let's do something fun. Let's change the solution to use Rational{BigInt} and print out the value at the end of the simulation. To do so, simply change the definition of the initial condition.

prob = ODEProblem(f,BigInt(1)//BigInt(2),(0//1,1//1),101//100);

sol =solve(prob,RK4(),dt=1//2^(6),adaptive=false)

println(sol[end])

That's one huge fraction!

Other Compatible Number Types

BigFloats

prob_ode_biglinear = ODEProblem(f,big(1.0)/big(2.0),(big(0.0),big(1.0)),big(1.01))

sol =solve(prob_ode_biglinear,Tsit5())

println(sol[end])

DoubleFloats.jl

There's are Float128-like types. Higher precision, but fixed and faster than arbitrary precision.

using Pkg

pkg"add DoubleFloats"

using DoubleFloats

prob_ode_doublelinear = ODEProblem(f,Double64(1)/Double64(2),(Double64(0),Double64(1)),Double64(1.01))

sol =solve(prob_ode_doublelinear,Tsit5())

println(sol[end])

ArbFloats

These high precision numbers which are much faster than Bigs for less than 500-800 bits of accuracy.

using Pkg

pkg"add ArbNumerics"

using ArbNumerics

prob_ode_arbfloatlinear = ODEProblem(f,ArbFloat(1)/ArbFloat(2),(ArbFloat(0.0),ArbFloat(1.0)),ArbFloat(1.01))

sol =solve(prob_ode_arbfloatlinear,Tsit5())

println(sol[end])

Incompatible Number Systems

DecFP.jl

Next let's try DecFP. DecFP is Intel Decimal Floating-Point Math Library, a fixed-precision decimals library which is made to give both performance but known decimals of accuracy.

using Pkg

pkg"add DecFP"

using DecFP

prob_ode_decfplinear = ODEProblem(f,Dec128(1)/Dec128(2),(Dec128(0.0),Dec128(1.0)),Dec128(1.01))

# StackOverflowError

sol =solve(prob_ode_decfplinear,Tsit5())

println(sol[end])

println(typeof(sol[end]))

Decimals.jl

using Pkg

pkg"add Decimals"

using Decimals

prob_ode_decimallinear = ODEProblem(f,[decimal("1.0")]./[decimal("2.0")],(0//1,1//1),decimal(1.01))

# MethodError: Decimal(::Rational{Int64}) is ambiguous

sol =solve(prob_ode_decimallinear,RK4(),dt=1/2^(6)) #Fails

println(sol[end])

println(typeof(sol[end]))

At the time of writing this, Decimals are not compatible. This is not on DifferentialEquations.jl's end, it's on partly on Decimal's end since it is not a subtype of Number. Thus it's not recommended you use Decimals with DifferentialEquations.jl

Part 1 Conclusion

As you can see, DifferentialEquations.jl can use arbitrary Julia-defined number systems in its arithmetic. If you need 128-bit floats, i.e. a bit more precision but not arbitrary, DoubleFloats.jl is a very good choice! For arbitrary precision, ArbNumerics are the most feature-complete and give great performance compared to BigFloats, and thus I recommend their use when high-precision (less than 512-800 bits) is required. DecFP is a great library for high-performance decimal numbers and works well as well. Other number systems could use some modernization.

Numbers with Uncertainties

By Mosè Giordano, Chris Rackauckas

The result of a measurement should be given as a number with an attached uncertainties, besides the physical unit, and all operations performed involving the result of the measurement should propagate the uncertainty, taking care of correlation between quantities.

There is a Julia package for dealing with numbers with uncertainties: Measurements.jl. Thanks to Julia's features, DifferentialEquations.jl easily works together with Measurements.jl out-of-the-box.

This notebook will cover some of the examples from the tutorial about classical Physics.

Caveat about Measurement type

Before going on with the tutorial, we must point up a subtlety of Measurements.jl that you should be aware of:

using Pkg

pkg"add Measurements"

using Measurements

5.23 ± 0.14 === 5.23 ± 0.14

(5.23± 0.14) - (5.23 ± 0.14)

(5.23 ± 0.14) / (5.23 ± 0.14)

x = 5.23 ± 0.14

x === x

x - x

x / x

Radioactive Decay of Carbon-14

The rate of decay of carbon-14 is governed by a first order linear ordinary differential equation inline_formula not implementedwhere inline_formula not implemented is the mean lifetime of carbon-14, which is related to the half-life inline_formula not implementedyears by the relation inline_formula not implemented.

using DifferentialEquations, Measurements, Plots

# Half-life and mean lifetime of radiocarbon, in years

t_12 = 5730 ± 40

τ = t_12 / log(2)

#Setup

u₀ = 1 ± 0

tspan = (0.0, 10000.0)

#Define the problem

radioactivedecay(u,p,t) = - u / τ

#Pass to solver

prob = ODEProblem(radioactivedecay, u₀, tspan)

sol = solve(prob, Tsit5(), reltol = 1e-8)

# Analytic solution

u = exp.(- sol.t / τ)

plot(sol.t, sol.u, label = "Numerical", xlabel = "Years", ylabel = "Fraction of Carbon-14")

plot!(sol.t, u, label = "Analytic")

The two curves are perfectly superimposed, indicating that the numerical solution matches the analytic one. We can check that also the uncertainties are correctly propagated in the numerical solution:

println("Quantity of carbon-14 after ",  sol.t[11], " years:")

println("Numerical: ", sol[11])

println("Analytic:  ", u[11])

Both the value of the numerical solution and its uncertainty match the analytic solution within the requested tolerance. We can also note that close to 5730 years after the beginning of the decay (half-life of the radioisotope), the fraction of carbon-14 that survived is about 0.5.

Simple pendulum

Small angles approximation

The next problem we are going to study is the simple pendulum in the approximation of small angles. We address this simplified case because there exists an easy analytic solution to compare.

The differential equation we want to solve is

inline_formula not implemented

where g = 9.97 ± 0.02 is the gravitational acceleration measured where the experiment is carried out, and L = 1.00 ± 0.01 is the length of the pendulum.

When you set up the problem for DifferentialEquations.jl remember to define the measurements as variables, as seen above.

using DifferentialEquations, Measurements, Plots

g = 9.79 ± 0.02; # Gravitational constants

L = 1.00 ± 0.01; # Length of the pendulum

#Initial Conditions

u₀ = [0 ± 0, π / 60 ± 0.01] # Initial speed and initial angle

tspan = (0.0, 6.3)

#Define the problem

function simplependulum(du,u,p,t)

    θ  = u[1]

    dθ = u[2]

    du[1] = dθ

    du[2] = -(g/L)*θ

end

#Pass to solvers

prob = ODEProblem(simplependulum, u₀, tspan)

sol = solve(prob, Tsit5(), reltol = 1e-6)

# Analytic solution

u = u₀[2] .* cos.(sqrt(g / L) .* sol.t)

plot(sol.t, getindex.(sol.u, 2), label = "Numerical")

plot!(sol.t, u, label = "Analytic")

Also in this case there is a perfect superimposition between the two curves, including their uncertainties.

We can also have a look at the difference between the two solutions:

plot(sol.t, getindex.(sol.u, 2) .- u, label = "")

Arbitrary amplitude

Now that we know how to solve differential equations involving numbers with uncertainties we can solve the simple pendulum problem without any approximation. This time the differential equation to solve is the following:

inline_formula not implemented

g = 9.79 ± 0.02; # Gravitational constants

L = 1.00 ± 0.01; # Length of the pendulum

#Initial Conditions

u₀ = [0 ± 0, π / 3 ± 0.02] # Initial speed and initial angle

tspan = (0.0, 6.3)

#Define the problem

function simplependulum(du,u,p,t)

    θ  = u[1]

    dθ = u[2]

    du[1] = dθ

    du[2] = -(g/L) * sin(θ)

end

#Pass to solvers

prob = ODEProblem(simplependulum, u₀, tspan)

sol = solve(prob, Tsit5(), reltol = 1e-6)

plot(sol.t, getindex.(sol.u, 2), label = "Numerical")

We note that in this case the period of the oscillations is not constant.

ArrayPartitions and Unitful

Heterogenous array and number with units.

ArrayPartitions in DiffEq are used for heterogeneous arrays. For example, PartitionedODEProblem solvers use them internally to turn the separate parts into a single array. You can construct an ArrayPartition using RecursiveArrayTools.jl

using RecursiveArrayTools

A = ArrayPartition(x::AbstractArray...)

where is x a list of arrays. The resulting A will act like a single array, and its broadcast will be type stable, allowing for it to be used inside of the native Julia DiffEq solvers in an efficient way. This is a good way to generate an array which has different units for different parts, or different amounts of precision.

An ArrayPartition acts like a single array. A[i] indexes through the first array, then the second, etc. all linearly. But A.x is where the arrays are stored. Thus for

using RecursiveArrayTools

A = ArrayPartition(y,z)

We would have A.x[1]==y and A.x[2]==z. Broadcasting like f.(A) is efficient.

In this example we will show using heterogeneous units in dynamics equations. Our arrays will be:

using Pkg

pkg"add RecursiveArrayTools Unitful"

using Unitful, RecursiveArrayTools, DifferentialEquations

using LinearAlgebra

r0 = [1131.340, -2282.343, 6672.423]u"km"

v0 = [-5.64305, 4.30333, 2.42879]u"km/s"

Δt = 86400.0*365u"s"

μ = 398600.4418u"km^3/s^2"

rv0 = ArrayPartition(r0,v0)

a

function f2(dy, y, μ, t)

    r = norm(y.x[1])

    dy.x[1] .= y.x[2]

    dy.x[2] .= -μ .* y.x[1] / r^3

end

Notice that y.x[1] is the r part of y, and y.x[2] is the v part of y. Using this kind of indexing is type stable, even though the array itself is heterogeneous. Note that one can also use things like 2y or y.+x and the broadcasting will be efficient.

Now to solve our equations, we do the same thing as always in DiffEq

prob = ODEProblem(f2, rv0, (0.0u"s", Δt), μ)

sol = solve(prob, Vern8())