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: 15057-element Vector{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}:
                  0.0 s
   0.5010872865302903 s
    5.511960151833193 s
    49.65740927941105 s
    338.4010306376774 s
    910.5445167643159 s
   1611.0318505574191 s
   2399.4869536144943 s
   3231.9860061580493 s
     4185.20951819496 s
                      ⋮
 3.1514491609223112e7 s
 3.1518170449560165e7 s
 3.1521121758305743e7 s
 3.1523149203188542e7 s
 3.1524943862757143e7 s
  3.152700008881182e7 s
 3.1529521278852012e7 s
 3.1533086565582093e7 s
             3.1536e7 s
u: 15057-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.7209889770912 km, -2065.6688381318254 km, 6783.977631064484 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.696876575795444 km s^-1, 4.421467943114541 km s^-1, 2.0631720933819193 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-808.7833764623475 km, -714.8417327609528 km, 7063.736257511328 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.703453416231231 km s^-1, 4.863823732841863 km s^-1, -0.1396087045695675 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-3741.5740892102035 km, 2029.407456276776 km, 5765.3963786746335 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-4.244479901848722 km s^-1, 4.443216961937508 km s^-1, -4.257431342225192 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5514.134607253488 km, 4362.262586088522 km, 1580.295029455917 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-0.6023199440638227 km s^-1, 1.9292826139225157 km s^-1, -7.15411634324551 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-4221.645811548964 km, 4364.012137590654 km, -3955.6285241256796 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[3.68759664275205 km s^-1, -1.9143884988177136 km s^-1, -6.114639326794478 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-70.5889661186245 km, 1469.0846191427606 km, -7107.082014499345 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[5.674642132596068 km s^-1, -4.613793657965086 km s^-1, -0.9981460834626812 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[4550.796069534674 km, -2907.6950545652585 km, -4796.477986929869 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[3.232954718797914 km s^-1, -3.8259890404465 km s^-1, 5.470212943412284 km s^-1])
 ⋮
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[13977.179033859908 km, -11405.904718614733 km, -2248.2472630222637 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[0.18566281718228972 km s^-1, -0.8758069983751592 km s^-1, 3.6229488103529475 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[8560.908195826289 km, -9239.29530769904 km, 9986.752001492026 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-3.075677615924994 km s^-1, 2.136856962945943 km s^-1, 2.3759142472541592 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-3067.4010941794672 km, 626.990453235967 km, 9955.137207985596 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-3.9934857956920657 km s^-1, 4.116676489607163 km s^-1, -3.683942030851453 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5661.563522361361 km, 5486.260996661042 km, -3457.8426543937308 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.53201821102877 km s^-1, -0.6654178339960191 km s^-1, -7.47186840404052 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[1857.0248540023833 km, 752.2450345556621 km, -11734.989625134802 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[4.619643422045183 km s^-1, -3.580427850141074 km s^-1, -1.6981399886837407 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[9740.795297541099 km, -6121.505554490888 km, -10782.582412634165 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.881818238002766 km s^-1, -2.8503938819986856 km s^-1, 2.0516276089566516 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[13894.034051997038 km, -11139.281104631034 km, -3237.3387490486666 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[0.43969233802947877 km s^-1, -1.0811757489675375 km s^-1, 3.572359884225057 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[9690.366780913819 km, -10000.415356786063 km, 8995.331261836007 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-2.7528667314283193 km s^-1, 1.796193672492318 km s^-1, 2.7135165171673807 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-1298.6744450709762 km, -1122.1420509898878 km, 11212.781337995435 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-4.288454578298883 km s^-1, 4.090147230433325 km s^-1, -2.2887667416896518 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.