DifferentialEquations.jl and systems with heterogenous units

diffeq

#1

Hi everybody,

I am trying to solve the following system of ODE with DifferentialEquations.jl

function newton!(t, y, dy, μ)
    r = norm(y[1:3])
    dy[1:3] = y[4:6]
    dy[4:6] = -μ * y[1:3] / r^3
end

Without units this works perfectly:

function test()
    r0 = [1131.340, -2282.343, 6672.423]
    v0 = [-5.64305, 4.30333, 2.42879]
    Δt = 86400.0*365
    mu = 398600.4418
    rv0 = [r0; v0]
    prob = ODEProblem((t, y, dy) -> newton!(t, y, dy, mu, rv0, (0.0, Δt))
    solve(prob, Vern8())[end]
end

The unitful version on the other hand does not:

function test2()
    r0 = [1131.340, -2282.343, 6672.423]u"km"
    v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
    Δt = 86400.0*365u"s"
    mu = 398600.4418u"km^3/s^2"
    rv0 = [r0; v0]

    prob = ODEProblem((t, y, dy) -> newton!(t, y, dy, mu), rv0, (0.0u"s", Δt))
    solve(prob, Vern8())[end]
end
TypeError: setindex!: in typeassert, expected Quantity{Float64, Dimensions:{}, Units:{}}, got Float64

The state vector rv0 with heterogenous units seems to be the problem.

I realise that I could avoid this problem by defining the system of ODE as separate functions in the non-mutating form but this is not an option because the above will be just a small part of the RHS in the complete system.

Is there a way I can have my cake and also eat it?


#2

Yes, but it’s not too simple. The ODE functions accept AbstractArrays, so we can build one for this purpose. As an example, DEDataArrays is one such supported array:

So let’s build one for your purpose. Let’s call it a HeterogenousArray. Since you have two different components, let’s make it have two parts:

type HeterogeneousArray{T,T2} <: AbstractVector
  r::T
  v::T2
end

Now we need to make it actually act like an array. To do this, we implement the array interface. See

http://docs.julialang.org/en/stable/manual/interfaces/#abstract-arrays

It should just need:

size(A::HeterogeneousArray)   = (length(A.x) + length(A.y),)
getindex( A::HeterogeniousArray,    i::Int) = (i <= length(A.x) ? A.x[i] : A.y[i-length(A.x)]))
setindex!(A::HeterogeniousArray, x, i::Int) = (i <= length(A.x) ? (A.x[i] = x) : A.y[i-length(A.x)] = x) )

Now you should be able to create this array:

    r0 = [1131.340, -2282.343, 6672.423]u"km"
    v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
    rv0 = HeterogeniousArray(r0,v0)

and then it should be that r0 == u0[1:3] and v0 == u0[4:6]. This means you can also simplify your ODE;

function newton!(t, y, dy, μ)
    r = norm(y[1:3])
    dy.r0 .= dy.v0
    dy.v0 .= -μ .* y[1:3] / r^3
end

Then you should be able to solve the ODE with this array:

prob = ODEProblem((t, y, dy) -> newton!(t, y, dy, mu), rv0, (0.0u"s", Δt))

That should work, but I didn’t debug this code so there may be a bug in this. But that should at least get you very close.

Caveats

I should explain the caveats here. Unitful.jl treats units as type information. This means that the eltype of your HeterogeneousArray is not concrete, and so accessing the values is not type stable. The way your newton! function is written here will be type stable, but using the A[i] iteration will not be. This will slow the solver down a bit.

Essentially, the solver does:

for i in eachindex(A)
  tmp[i] = tmp[i] + some_constant*dy[i]
end

Maybe there is a way to use broadcast:

tmp .= tmp .+ some_constant*dy

and have it iterate fast on such heterogenous arrays? Maybe @stevengj would know. That would be a good trick to find out!

If we can find out how to do some kind of “fast type-stable broadcast” here, and switch the internal for loops to broadcast (I’m waiting on stack-allocated views, as until then there is a performance issue with small broadcasts), then this approach should be fine performance-wise.


#3

Excellent! Thanks, Chris!

I was planning to do something along this line anyway since the r and v are already part of a State type. I will try this out and post here again if run into trouble.


#4

I you don’t care too much about performance, and easier way to do this would be to just not have a strictly typed array.

    rv0 = Vector{Number}()
    r0 = [1131.340, -2282.343, 6672.423]u"km"
    v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
    append!(rv0,r0); append!(rv0,v0)

But I noticed that

     Δt = 86400.0*365u"s"
    mu = 398600.4418u"km^3/s^2"
    prob = ODEProblem((t, y, dy) -> newton!(t, y, dy, mu), rv0, (0.0u"s", Δt))
    solve(prob, Vern8())[end]

this errors.

ERROR: DimensionError: 1//1000000 and 1.13134 km are not dimensionally compatible.
 in .+(::Rational{Int64}, ::Unitful.Quantity{Float64,Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)},Unitful.FreeUnits{(Unitful.Unit{:Meter,Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}}(3,1//1),),Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}}}) at .\operators.jl:152
 in collect(::Base.Generator{Array{Unitful.Quantity{Float64,D,U},1},Base.##143#145{Rational{Int64}}}) at .\array.jl:307
 in .+(::Rational{Int64}, ::Array{Unitful.Quantity{Float64,D,U},1}) at .\arraymath.jl:90
 in +(::Rational{Int64}, ::Array{Unitful.Quantity{Float64,D,U},1}) at .\arraymath.jl:130
 in ode_determine_initdt(::Array{Number,1}, ::Unitful.Quantity{Float64,Unitful.Dimensions{(Unitful.Dimension{:Time}(1//1),)},Unitful.FreeUnits{(Unitful.Unit{:Second,Unitful.Dimensions{(Unitful.Dimension{:Time}(1//1),)}}(0,1//1),),Unitful.Dimensions{(Unitful.Dimension{:Time}(1//1),)}}}, ::Float64, ::Unitful.Quantity{Float64,Unitful.Dimensions{(Unitful.Dimension{:Time}(1//1),)},Unitful.FreeUnits{(Unitful.Unit{:Second,Unitful.Dimensions{(Unitful.Dimension{:Time}(1//1),)}}(0,1//1),),Unitful.Dimensions{(Unitful.Dimension{:Time}(1//1),)}}}, ::Rational{Int64}, ::Float64, ::OrdinaryDiffEq.#ODE_DEFAULT_NORM, ::DiffEqBase.ODEProblem{Array{Number,1},Quantity{Float64, Dimensions:{�}, Units:{s}},true,##3#4}, ::Int64) at C:\Users\Chris\.julia\v0.5\OrdinaryDiffEq\src\initdt.jl:5

This makes me think that the previous solution won’t work either. This isn’t going to work because it doesn’t know what to type some of the internal constants, and it usually types them using eltype(u0). In this case, eltype(u0) is not concrete, so it dies. To handle this case, the solvers will need to generalize how it does the typing.

Would you mind opening an issue? It’ll take some work to make heterogenous arrays like this work. And I think that with a proper broadcast (making them an iterator in the right), they can even work well. But it’ll take some work.


#5

Done: https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/147