Solving differential equations with automated unit checking

Jakob Peder Pettersen

June 10, 2026

Why do we need quantity calculus?

  • Computers work with numbers
  • Scientists work with measurements
  • With pen and paper: Show your calculations!
  • With computers: Comments in code, no enforcement

A horror example

  • Mars Climate Orbiter 1998
  • Thrust force calculations
  • SI (metric) units mixed up with US customary units
  • Orbiter crashed

NASA/JPL-Caltech

Unitful.jl

  • A Julia package for handling measurements

  • Templated types

  • Unit checking reduced to type checking

Examples

using Unitful

Compatiable, but different units

length_sum = 3u"m" + 6u"ft"
length_sum |> float |> println
4.8288 m

Unit composition and conversion

speed = 15u"m" / 3u"s"
println(speed)
speed |> u"km/hr" |> println
5.0 m s^-1
18.0 km hr^-1

Mixing up multiplication and division

try
  speed = 15u"m" * 3u"s"
  speed |> u"km/hr" |> println
catch e
  e
end
Unitful.DimensionError(km hr^-1, m s)

Adding incompatible quantities

try
  length_area_sum = 2u"m" + 8u"m"^2
catch e
  e
end
Unitful.DimensionError(2 m, 8 m^2)

How does this work for differential equations?

Orbital mechanics ODE:

\(\dfrac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} = \mathbf{v}\)

\(\dfrac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} = - \frac{\mu\cdot \mathbf{r}}{\left\|\mathbf{v}\right\|_2^3}\)

  • \(\left[\mathbf{r}\right] = \mathrm{km}\)
  • \(\left[\mathbf{v}\right] = \mathrm{km/s}\)
  • \(\left[\mu\right] = \mathrm{km^3/s^2}\)
  • \(\left[t\right] = \mathrm{s}\)

Which problems must be addressed?

Performance

  • The process of checking the units should only happen once
  • Independent on the size of \(\mathbf{r}\) and \(\mathbf{v}\)
  • Must be statically inferable from the types of the provided values

Handling heterogenity

  • The system consists of two vectors with different units of measurement
  • Still the ODE solver expects a single array with the system variables

\(\left[\mathbf{r}\right] = \mathrm{km}\) and \(\left[\mathbf{v}\right] = \mathrm{km/s}\), but what is \(\left[\begin{pmatrix}\mathbf{r} \\ \mathbf{v}\end{pmatrix}\right]\)?

Propagation of units through the solver

  • The solver must know that \(\left[\dfrac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}\right] = \mathrm{km/s^2}\)
  • Operations must either be dimensionaly consistent such as \(\mathbf{u}_{\mathrm{res}}= k_1\mathbf{u}_1 + k_2\mathbf{u}_2\)…
  • …or units must be stripped such as when computing \(\left\|\mathbf{u}\right\|\)
  • Everything which touches the units must be able to support it

How is it solved?

  • The package DifferentialEquations supports Unitful for some ODE solvers

Tricks of the trade

  • Element-wise operations are not performed by explicit loops (u_res[i] = k_1 * u_1[i] + k_2 * u_2[i]), but instead by broadcasting u .= k_1 .* u_1 .+ k_2 .* u_2
  • No explicit constructurs, but general overloadable functions such as zero()
  • Some functions are explicitly overloaded to behave differently when encountering Unitful quantities. Examples:
    • Determining the element type of the derivative
    • Computing vector norm

RecursiveArrayTools.jl

Supports both heterogeneous element types and type-stable broadcasting

using DifferentialEquations, Unitful
using RecursiveArrayTools
using LinearAlgebra
r0 = [1131.340, -2282.343, 6672.423]Unitful.u"km"
v0 = [-5.64305, 4.30333, 2.42879]Unitful.u"km/s"
Δt = 86400.0 * 365Unitful.u"s"
μ = 398600.4418Unitful.u"km^3/s^2"
rv0 = RecursiveArrayTools.ArrayPartition(r0, v0)
function f_orbital!(dy, y, μ, t)
    r = norm(y.x[1])
    dy.x[1] .= y.x[2]
    dy.x[2] .= -μ .* y.x[1] / r^3
end

prob = ODEProblem(f_orbital!, rv0, (0.0Unitful.u"s", Δt), μ)
sol = solve(prob, Vern8())

But we can do better

  • ArrayPartition is supposed to provide type-stable broadcasting, but sometimes the compiler gives up on inference, causing considerable slowdowns
  • ArrayPartition does not allow the components to have names, only a numeric index

HeterogeneousArrays.jl

Broadcast unpacking

Consider arrays y, x_1, x_2 with components u and v. We decompose the the heterogeneous broadcast y .= a .* x_1 .+ b .* x_2 into the homogeneous broadcasts y.u .= a .* x_1.u .+ b .* x_2.u and y.v .= a .* x_1.v .+ b .* x_2.v

The example revisited

using HeterogeneousArrays
function f_orbital_het!(dy, y, μ, t)
    r_mag = norm(y.r)
    dy.r .= y.v
    dy.v .= -μ .* y.r ./ r_mag^3
    dy
    return dy
end
u0 = HeterogeneousVector(r=r0, v=v0)
prob = ODEProblem(f_orbital_het!, u0, (0.0Unitful.u"s", Δt), μ)
sol = solve(prob, Vern8())

Remaining challenges and future directions

  • Can we develop an array type which reinterprets a plain numeric Array as having Unitful components?
  • Implicit solvers do still not support Unitful: Can we construct a wrapper to make any solver work?