Unit Checked Arithmetic via SIUnits

Units and dimensional analysis are standard tools across the sciences for checking the correctness of your equation. However, most ODE solvers only allow for the equation to be in dimensionless form, leaving it up to the user to both convert the equation to a dimensionless form, punch in the equations, and hopefully not make an error along the way.

DifferentialEquations.jl allows for one to use Unitful.jl to have unit-checked arithmetic natively in the solvers. Given the dispatch implementation of the Unitful, this has little overhead.

Using Unitful

To use SIUnits, you need to have the package installed. Then you can add units to your variables. For example:

In [1]:
using Unitful
t = 1.0u"s"
INFO: Precompiling module Unitful.
Out[1]:
1.0 s

Notice that t is a variable with units in seconds. If we make another value with seconds, they can add

In [2]:
t2 = 1.02u"s"
t+t2
Out[2]:
2.02 s

and they can multiply:

In [3]:
t*t2
Out[3]:
1.02 s^2

You can even do rational roots:

In [4]:
sqrt(t)
Out[4]:
1.0 s^1/2

Many operations work.

These operations will check to make sure units are correct, and will throw an error for incorrect operations:

In [5]:
t + sqrt(t)
DimensionError: 1.0 s and 1.0 s^1/2 are not dimensionally compatible.

Stacktrace:
 [1] +(::Quantity{Float64, Dimensions:{𝐓}, Units:{s}}, ::Quantity{Float64, Dimensions:{𝐓^1/2}, Units:{s^1/2}}) at /home/crackauc/.julia/v0.6/Unitful/src/Unitful.jl:474
 [2] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/crackauc/.julia/v0.6/IJulia/src/execute_request.jl:160
 [3] eventloop(::ZMQ.Socket) at /home/crackauc/.julia/v0.6/IJulia/src/eventloop.jl:8
 [4] (::IJulia.##11#14)() at ./task.jl:335

Using Unitful with DifferentialEquations.jl

Just like with other number systems, you can choose the units for your numbers by simply specifying the units of the initial condition and the timestep. For example, to solve the linear ODE where the variable has units of Newton's and t is in Seconds, we would use:

In [10]:
using DifferentialEquations
f = (t,y) -> 0.5*y
u0 = 1.5u"N"
prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))
sol = solve(prob,Tsit5())
DimensionError: 1.5 N and 0.015 N s are not dimensionally compatible.

Stacktrace:
 [1] ode_determine_initdt(::Quantity{Float64, Dimensions:{𝐋 𝐌 𝐓^-2}, Units:{N}}, ::Quantity{Float64, Dimensions:{𝐓}, Units:{s}}, ::Float64, ::Quantity{Float64, Dimensions:{𝐓}, Units:{s}}, ::Quantity{Float64, Dimensions:{𝐋 𝐌 𝐓^-2}, Units:{N}}, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.ODEProblem{Quantity{Float64, Dimensions:{𝐋 𝐌 𝐓^-2}, Units:{N}},Quantity{Float64, Dimensions:{𝐓}, Units:{s}},false,##5#6,Void,UniformScaling{Int64}}, ::Int64) at /home/crackauc/.julia/v0.6/OrdinaryDiffEq/src/initdt.jl:68
 [2] #init#134(::Int64, ::Array{Quantity{Float64, Dimensions:{𝐓}, Units:{s}},1}, ::Array{Quantity{Float64, Dimensions:{𝐓}, Units:{s}},1}, ::Array{Quantity{Float64, Dimensions:{𝐓}, Units:{s}},1}, ::Void, ::Bool, ::Void, ::Bool, ::Bool, ::Bool, ::Quantity{Float64, Dimensions:{𝐓}, Units:{s}}, ::Bool, ::Rational{Int64}, ::Void, ::Void, ::Int64, ::Rational{Int64}, ::Rational{Int64}, ::Bool, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Quantity{Float64, Dimensions:{𝐓}, Units:{s}}, ::Quantity{Float64, Dimensions:{𝐓}, Units:{s}}, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN, ::DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE, ::Void, ::Void, ::Bool, ::Bool, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Quantity{Float64, Dimensions:{𝐋 𝐌 𝐓^-2}, Units:{N}},Quantity{Float64, Dimensions:{𝐓}, Units:{s}},false,##5#6,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/crackauc/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:137
 [3] init(::DiffEqBase.ODEProblem{Quantity{Float64, Dimensions:{𝐋 𝐌 𝐓^-2}, Units:{N}},Quantity{Float64, Dimensions:{𝐓}, Units:{s}},false,##5#6,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/crackauc/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:46
 [4] #solve#133(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Quantity{Float64, Dimensions:{𝐋 𝐌 𝐓^-2}, Units:{N}},Quantity{Float64, Dimensions:{𝐓}, Units:{s}},false,##5#6,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/crackauc/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [5] solve(::DiffEqBase.ODEProblem{Quantity{Float64, Dimensions:{𝐋 𝐌 𝐓^-2}, Units:{N}},Quantity{Float64, Dimensions:{𝐓}, Units:{s}},false,##5#6,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.Tsit5) at /home/crackauc/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [6] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/crackauc/.julia/v0.6/IJulia/src/execute_request.jl:160
 [7] eventloop(::ZMQ.Socket) at /home/crackauc/.julia/v0.6/IJulia/src/eventloop.jl:8
 [8] (::IJulia.##11#14)() at ./task.jl:335

Notice that we recieved a unit mismatch error. This is correctly so! Remember that for an ODE:

$$ \frac{dy}{dt} = f(t,y) $$

we must have that f is a rate, i.e. f is a change in y per unit time. So we need to fix the units of f in our example to be N/s. Notice that we then do not receive an error if we do the following:

In [11]:
f = (t,y) -> 0.5*y/3.0u"s"
prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))
sol = solve(prob,Tsit5())
Out[11]:
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 3-element Array{Quantity{Float64, Dimensions:{𝐓}, Units:{s}},1}:
      0.0 s
 0.143116 s
      1.0 s
u: 3-element Array{Quantity{Float64, Dimensions:{𝐋 𝐌 𝐓^-2}, Units:{N}},1}:
     1.5 N
 1.53621 N
 1.77204 N

This gives a a normal solution object. Notice that the values are all with the correct units:

In [12]:
print(sol[:])
Quantity{Float64, Dimensions:{𝐋 𝐌 𝐓^-2}, Units:{N}}[1.5 N, 1.53621 N, 1.77204 N]

We can plot the solution using the plot recipe.

In [20]:
#Pkg.clone("https://github.com/ajkeller34/UnitfulPlots.jl")
using UnitfulPlots, Plots
gr()
plot(sol.t,sol[:],lw=3)
Out[20]:
0.0s 0.2s 0.4s 0.6000000000000001s 0.8s 1.0s 1.5N 1.6N 1.7000000000000002N y1

Notice that here we pulled the units for the label directly from the solution. Thus if the units change, the labels will change automatically.

Other Solvers and Future Developments

As of right now, the ODE solvers are the only solvers which are fully compatible with units. The SDE solvers will come shortly. For the FEM PDE solvers, the tooling is all compatible with units (i.e. you can make meshes with units, and most of the functions will work). However, the solver step is what's not compatible with units. The reason is because \ uses CHOLMOD which does not work with units, and sparse multiplication also is undefined for units. These facts are major stop gaps in development here. I believe it's not worth it to simply "turn off" units at that point because that is there area where one would wish to have units checked.