DiffEq-Specific Array Types

In many cases, a standard array may not be enough to fully hold the data for a model. Many of the solvers in DifferentialEquations.jl (only the native Julia methods) allow you to solve problems on AbstractArray types, which allow you to extend the meaning of an array. This page describes some AbstractArray types which can be helpful for modeling differential equations problems.

ArrayPartitions

ArrayPartitions in DiffEq are used for heterogeneous arrays. For example, DynamicalODEProblem 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 x is an array of arrays. Then, A will act like a single array, and its broadcast will be type stable, allowing for it to be efficiently used inside the native Julia DiffEq solvers. This is a good way to generate an array which has different units for different parts, or different amounts of precision.

Usage

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.

Example: Dynamics Equations

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

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)

(Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[1131.34 km, -2282.343 km, 6672.423 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.64305 km s^-1, 4.30333 km s^-1, 2.42879 km s^-1])

Here, r0 is the initial positions, and v0 are the initial velocities. rv0 is the ArrayPartition initial condition. We now write our update function in terms of the ArrayPartition:

function f(dy, y, μ, t)
    r = norm(y.x[1])
    dy.x[1] .= y.x[2]
    dy.x[2] .= -μ .* y.x[1] / r^3
end

f (generic function with 1 method)

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(f, rv0, (0.0u"s", Δt), μ)
sol = solve(prob, Vern8())

retcode: Success
Interpolation: specialized 8th order lazy interpolation
t: 14944-element Vector{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}:
    0.0 s
    0.5010872865302903 s
    5.511960151833193 s
   49.657408270019076 s
  342.0794963810424 s
  916.0512373944414 s
 1617.4717343252378 s
 2407.137427116769 s
 3240.787644052227 s
 4194.93521388548 s
    ⋮
    3.1516661045151204e7 s
    3.1518709694033787e7 s
    3.1520518728919588e7 s
    3.1522593272173803e7 s
    3.1525142048289914e7 s
    3.1528731743771587e7 s
    3.1532110540768426e7 s
    3.1534180023415577e7 s
    3.1536e7 s
u: 14944-element Vector{RecursiveArrayTools.ArrayPartition{Unitful.Quantity{Float64}, Tuple{Vector{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}}, Vector{Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}}}}}:
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[1131.34 km, -2282.343 km, 6672.423 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.64305 km s^-1, 4.30333 km s^-1, 2.42879 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[1128.5121841206073 km, -2280.186342651942 km, 6673.639119236024 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.6436694625389 km s^-1, 4.304580664105494 km s^-1, 2.4251316255351654 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[1100.217102647565 km, -2258.5854152557135 km, 6685.699423475585 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.649778787886409 km s^-1, 4.317022208036712 km s^-1, 2.3885114106020557 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[849.7209947274717 km, -2065.668842594817 km, 6783.97762898194 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.696876574857083 km s^-1, 4.42146794083337 km s^-1, 2.0631721008736457 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-829.7573051670589 km, -696.9450845563044 km, 7063.170505771802 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.700161226692357 km s^-1, 4.866660324261165 km s^-1, -0.1679927773238569 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-3764.885688994574 km, 2053.84155310996 km, 5741.8573557461605 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-4.222097610369185 km s^-1, 4.431041786662769 km s^-1, -4.291742954513994 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5517.891670905848 km, 4374.590509382225 km, 1534.1887985802746 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-0.5644886852632467 km s^-1, 1.89932201315544 km s^-1, -7.164796604459822 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-4193.304847938871 km, 4349.2325457435845 km, -4002.286653361141 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[3.721310875919539 km s^-1, -1.9492978327292696 km s^-1, -6.0827562160652775 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-20.640690800630747 km, 1428.4169279464832 km, -7115.580153101964 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[5.675060779556315 km s^-1, -4.6270894761379 km s^-1, -0.9328823788287788 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[4582.010449638985 km, -2944.7592161833245 km, -4743.037387101988 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[3.185971415515786 km s^-1, -3.7958813883094744 km s^-1, 5.519288692025722 km s^-1])
 ⋮
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-3276.404855752255 km, 803.8463110851006 km, 9956.978726518842 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-3.893804431320847 km s^-1, 4.056219407059085 km s^-1, -3.8053123375426847 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5624.904844611834 km, 5495.083271163315 km, -3659.09984963823 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.650104565630201 km s^-1, -0.7826519280315944 km s^-1, -7.385604422957212 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[2054.6051154998827 km, 593.070575673022 km, -11777.152064833706 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[4.6240928041891936 km s^-1, -3.610339419697076 km s^-1, -1.5663158226753748 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[9950.699336925276 km, -6338.736618068762 km, -10584.652039362652 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.831238262794342 km s^-1, -2.824595296177967 km s^-1, 2.13781553805108 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[13987.155449447837 km, -11300.133086903978 km, -2824.3411435318444 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[0.36952324508114337 km s^-1, -1.0249321893471133 km s^-1, 3.5887759698903303 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[9558.672876843153 km, -9963.453058542756 km, 9372.087332843888 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-2.7808135460980687 km s^-1, 1.835639977949965 km s^-1, 2.6340890320278527 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-3279.594695985543 km, 807.5427383880175 km, 9951.97755268997 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-3.8937162912882055 km s^-1, 4.056751970751892 km s^-1, -3.8083750982447 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5566.682192214814 km, 5477.903571131761 km, -3821.4374288888907 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.7236621461768773 km s^-1, -0.8549733168846072 km s^-1, -7.335426458686554 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[2203.0435753224933 km, 477.7053688764 km, -11830.108745694055 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[4.607048025473749 km s^-1, -3.613115820725251 km s^-1, -1.4794247847055464 km s^-1])

MultiScaleArrays

The multi-scale modeling functionality is provided by MultiScaleArrays.jl. It allows for designing a multi-scale model as an extension of an array, which in turn can be directly used in the native Julia solvers of DifferentialEquations.jl.

For more information, please see the MultiScaleArrays.jl README.