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:

import RecursiveArrayTools
A = RecursiveArrayTools.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

import RecursiveArrayTools
A = RecursiveArrayTools.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:

import Unitful, RecursiveArrayTools, DifferentialEquations as DE
import 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)

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

retcode: Success
Interpolation: specialized 8th order lazy interpolation
t: 15105-element Vector{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}:
    0.0 s
    0.5010872865302903 s
    5.511960151833193 s
   49.65749501743387 s
  331.9280101507043 s
  900.7742864610404 s
 1599.5825166217742 s
 2385.82928614952 s
 3217.9099994243315 s
 4168.9716554101815 s
    ⋮
    3.1517664278677177e7 s
    3.1519529327302657e7 s
    3.1521625344985675e7 s
    3.1524260873042308e7 s
    3.1527931889868982e7 s
    3.1530880809101906e7 s
    3.1532902114655253e7 s
    3.153469351671248e7 s
    3.1536e7 s
u: 15105-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.7205005381535 km, -2065.668459043897 km, 6783.977807956757 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.696876655500711 km s^-1, 4.421468136877816 km s^-1, 2.0631714570331234 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-771.8465767290205 km, -746.3087626831779 km, 7064.478273234513 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.709042215173237 km s^-1, 4.858657387750038 km s^-1, -0.08965408508251303 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-3699.911845487239 km, 1985.8920466193579 km, 5806.693713054252 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-4.283854081052035 km s^-1, 4.464462596763785 km s^-1, -4.196201407573931 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5506.853642558792 km, 4339.86955569903 km, 1662.0926494816074 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-0.6695275734710813 km s^-1, 1.9823495072687172 km s^-1, -7.134343750504316 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-4271.59569088355 km, 4389.731023817644 km, -3871.733146093691 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[3.6268393274665396 km s^-1, -1.8517676697886847 km s^-1, -6.17063313436589 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-150.45517312706122 km, 1533.8743931296285 km, -7092.298525715554 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[5.673020060289419 km s^-1, -4.591757419674905 km s^-1, -1.1023435702637705 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[4497.667994006312 km, -2845.1669071977008 km, -4884.629475714578 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[3.3106507670897765 km s^-1, -3.875386833834127 km s^-1, 5.387084334972934 km s^-1])
 ⋮
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5633.980947079544 km, 5474.600474171505 km, -3516.987878625412 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.545753809591989 km s^-1, -0.6787668394028044 km s^-1, -7.463284053377219 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[2199.506355656813 km, 490.3629859399043 km, -11878.817819557236 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[4.571463421817369 km s^-1, -3.5786464964964475 km s^-1, -1.5010905127742338 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[10046.809887721993 km, -6418.869850453851 km, -10591.520881812608 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.7660685862917713 km s^-1, -2.7738114246206176 km s^-1, 2.160387005627397 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[13953.783720758354 km, -11338.914450954462 km, -2486.049192445895 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[0.2269663193308073 km s^-1, -0.9091779482385521 km s^-1, 3.614619314358947 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[8674.556874555921 km, -9298.945588703531 km, 9801.584049171322 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-3.055129195966569 km s^-1, 2.1092213694362814 km s^-1, 2.42741542914114 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-2967.8517306703266 km, 540.2526077481538 km, 9966.871518602411 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-4.039473773447032 km s^-1, 4.142800770518627 km s^-1, -3.6190335268945595 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5687.55193011198 km, 5489.001604799489 km, -3360.529355664488 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.4696281144160914 km s^-1, -0.6052752674182666 km s^-1, -7.508381440733178 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[1760.37975716595 km, 832.7540673850431 km, -11727.730134829206 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[4.613576553464361 km s^-1, -3.562484895615996 km s^-1, -1.762686683765369 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[7213.338612210131 km, -3728.992720697248 km, -12040.371567766482 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[3.6258714711068567 km s^-1, -3.276450803636394 km s^-1, 1.018521743512115 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.