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: 14925-element Vector{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}:
                  0.0 s
   0.5010872865302903 s
    5.511960151833193 s
    49.65749620961972 s
   344.09099029529807 s
    919.0533764696446 s
   1620.9797864665975 s
    2411.297174606073 s
   3245.7215806731824 s
    4200.321308253043 s
                      ⋮
   3.15169528456449e7 s
 3.1519082125625703e7 s
  3.152177180598552e7 s
 3.1525492491336733e7 s
 3.1528468809713744e7 s
  3.153052301989648e7 s
 3.1532337301919535e7 s
 3.1534417755378913e7 s
             3.1536e7 s
u: 14925-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.7204937464173 km, -2065.668453772683 km, 6783.977810416442 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.696876656609012 km s^-1, 4.421468139572094 km s^-1, 2.0631714481847045 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-841.2213025895506 km, -687.1542939115907 km, 7062.8169801172 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.698325368328203 km s^-1, 4.868180996349745 km s^-1, -0.18351261265070093 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-3777.542619913166 km, 2067.134112110032 km, 5728.944948604166 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-4.209837923943154 km s^-1, 4.424343430911492 km s^-1, -4.3103879373428 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5519.835765354315 km, 4381.22474932109 km, 1509.0443248053725 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-0.5438725039197423 km s^-1, 1.8829679700495678 km s^-1, -7.170480733822731 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-4177.78719105982 km, 4341.084576300797 km, -4027.5530248936416 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[3.739544894847063 km s^-1, -1.9682271895725276 km s^-1, -6.065265528516156 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[7.35981319465464 km, 1405.5691509791045 km, -7120.092589718263 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[5.675094956563403 km s^-1, -4.634379554130762 km s^-1, -0.8962629329643912 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[4599.099903766753 km, -2965.1588129394418 km, -4713.237400575327 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[3.15981000502549 km s^-1, -3.779041198598221 km s^-1, 5.5462342618575216 km s^-1])
 ⋮
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[2528.887705991007 km, 224.8647124452155 km, -11948.388051998214 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[4.565660151224991 km s^-1, -3.6151467050732093 km s^-1, -1.2921946144027088 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[10384.533896362598 km, -6767.8485899879415 km, -10275.75516681439 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.678418597403076 km s^-1, -2.7248288063545396 km s^-1, 2.288175710976672 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[14080.937816554622 km, -11580.7537979118 km, -1810.1478676632807 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[0.10840454520311822 km s^-1, -0.8113345142380913 km s^-1, 3.6281800929569052 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[8381.542244077935 km, -9165.23132596768 km, 10380.256545910703 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-3.105339044707374 km s^-1, 2.182439338506913 km s^-1, 2.2728733445532305 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-3267.518318071178 km, 791.6638201660726 km, 9980.416213912104 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-3.8971982694778275 km s^-1, 4.0563437739246835 km s^-1, -3.791426397696102 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5641.111005389567 km, 5507.423526867451 km, -3652.031907848906 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.6378241293739917 km s^-1, -0.7725964758708339 km s^-1, -7.383801340663307 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[2043.766577548585 km, 607.1767706461003 km, -11801.94386418947 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[4.6183267730878175 km s^-1, -3.603740687636387 km s^-1, -1.5749373843593577 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[9957.692625433827 km, -6336.485222314852 km, -10625.911924688176 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.831943936623866 km s^-1, -2.8233353973874222 km s^-1, 2.12844387805993 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[13207.331596820251 km, -9964.830981267367 km, -6223.808334445178 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[1.28639473314743 km s^-1, -1.740551531917572 km s^-1, 3.2769755293028027 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.