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: 14946-element Vector{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}:
    0.0 s
    0.5010872865302903 s
    5.805979586113739 s
   49.989903364138584 s
  342.1028180992968 s
  915.869325716335 s
 1617.2063447343583 s
 2406.7339075637497 s
 3240.2704452544294 s
 4194.3222362074885 s
    ⋮
    3.151729124094261e7 s
    3.1519365480169553e7 s
    3.152191390951259e7 s
    3.1525503041856922e7 s
    3.1528881409946427e7 s
    3.153095060645922e7 s
    3.1532841062493805e7 s
    3.1534964837375272e7 s
    3.1536e7 s
u: 14946-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}}[1098.555905880957 km, -2257.316020038803 km, 6686.401376068454 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.650132441313096 km s^-1, 4.317748551324979 km s^-1, 2.3863606250269203 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[847.826759865609 km, -2064.1986013034525 km, 6784.663213307144 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.697185330213727 km s^-1, 4.422219094190918 km s^-1, 2.0607041935462735 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-829.8902424751218 km, -696.831585467906 km, 7063.166585793368 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-5.700140085490734 km s^-1, 4.866678078701444 km s^-1, -0.16817272380932363 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-3764.1175728249627 km, 2053.0354582011378 km, 5742.637970968562 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-4.222839175522549 km s^-1, 4.431446291305348 km s^-1, -4.290611795217702 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5517.741654767647 km, 4374.086285347134 km, 1536.0902028261482 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-0.5660481066594456 km s^-1, 1.9005582696498062 km s^-1, -7.164362749812324 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-4194.80611096071 km, 4350.01875437514 km, -3999.8318042908422 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[3.7195384180363926 km s^-1, -1.9474596293494506 km s^-1, -6.084447109851014 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-23.575817425478952 km, 1430.8098498782806 km, -7115.096676930759 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[5.675048855398431 km s^-1, -4.6263185023195605 km s^-1, -0.9367195916641606 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[4580.056621961463 km, -2942.431853521859 km, -4746.419631625624 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[3.188942380161787 km s^-1, -3.7977904132393756 km s^-1, 5.5162115531788505 km s^-1])
 ⋮
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[2055.3576079642003 km, 592.1726743202313 km, -11775.841629115048 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[4.62440870536565 km s^-1, -3.6107222856571077 km s^-1, -1.565735832974777 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[9950.496538592994 km, -6339.042031282196 km, -10582.244530015407 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.8311564631764563 km s^-1, -2.82464129980446 km s^-1, 2.138397345894639 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[13985.8198619486 km, -11299.50213521092 km, -2821.8117764094077 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[0.36918337373658 km s^-1, -1.0247058263268203 km s^-1, 3.5890877651307402 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[9556.738197133409 km, -9961.964021263384 km, 9372.851075156326 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-2.781133062032331 km s^-1, 1.8359596348346023 km s^-1, 2.633843278131105 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-3280.301826398145 km, 808.4024713814631 km, 9950.665627251161 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[-3.893430082528038 km s^-1, 4.056686823667068 km s^-1, -3.8092704653613163 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[-5565.736613827789 km, 5477.222166951774 km, -3822.0445427547893 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.7244604019558754 km s^-1, -0.8556535174058069 km s^-1, -7.335409641259811 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[2528.592725381743 km, 220.85860168396263 km, -11926.922847104552 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[4.569815913989794 km s^-1, -3.618988889062963 km s^-1, -1.2905889511083415 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[10369.637143990009 km, -6759.894583807411 km, -10252.16577713465 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[2.6797213754706637 km s^-1, -2.7266967415709824 km s^-1, 2.2920250523489005 km s^-1])
 (Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}[12613.32042371958 km, -9221.432050040974 km, -7432.7621034178355 km], Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(km, s^-1), 𝐋 𝐓^-1, nothing}}[1.6596628362611132 km s^-1, -2.0188728276648806 km s^-1, 3.0844008995511976 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.