# DifferentialEquations.jl and systems with heterogenous units

**URL:** <https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952>\
**Category:** Numerics\
**Tags:** diffeq, unitful\
**Created:** [March 29, 2017, 5:28pm UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952 "2017-03-29T17:28:02Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![helgee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgee/32/2022_2.png) [@helgee](https://discourse.julialang.org/u/helgee)\
**Post date:** [March 29, 2017, 5:28pm UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/1 "2017-03-29T17:28:02Z")

</div>

Hi everybody,

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

```julia
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:

```julia
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:

```julia
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

```

```julia
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 29, 2017, 8:30pm UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/2 "2017-03-29T20:30:17Z")

</div>

> [@helgee](#):
>
> Is there a way I can have my cake and also eat it?

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

[https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/data\_array.jl](https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/data_array.jl)

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:

```julia
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](http://docs.julialang.org/en/stable/manual/interfaces/#abstract-arrays)

It should just need:

```julia
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:

```julia
    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;

```julia
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:

```julia
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:

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

```

Maybe there is a way to use broadcast:

```julia
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.

---

<div class="post-metadata">

**Author:** ![helgee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgee/32/2022_2.png) [@helgee](https://discourse.julialang.org/u/helgee)\
**Post date:** [March 30, 2017, 4:09am UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/3 "2017-03-30T04:09:00Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 30, 2017, 4:26am UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/4 "2017-03-30T04:26:45Z")

</div>

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.

```julia
    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

```nohighlight
     Δ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.

```julia
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.

---

<div class="post-metadata">

**Author:** ![helgee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/helgee/32/2022_2.png) [@helgee](https://discourse.julialang.org/u/helgee)\
**Post date:** [March 30, 2017, 4:37am UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/5 "2017-03-30T04:37:34Z")

</div>

Done: [https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/147](https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/147)

---

<div class="post-metadata">

**Author:** ![ggggggggg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ggggggggg/32/265_2.png) [@ggggggggg](https://discourse.julialang.org/u/ggggggggg)\
**Post date:** [July 13, 2020, 6:16pm UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/6 "2020-07-13T18:16:55Z")

</div>

The example in the issue no longer works, I was able to adapt it to work as follows

```julia
using Unitful, RecursiveArrayTools, DiffEqBase, OrdinaryDiffEq, LinearAlgebra

function newton(u, p, t)
    mu = 398600.4418u"km^3/s^2"
    r = norm(u.x[1])
    ArrayPartition(u.x[2], -mu .* u.x[1] / r^3)
end

r0 = [1131.340, -2282.343, 6672.423]u"km"
v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
Δt = 86400.0*365u"s"
rv0 = ArrayPartition(r0,v0)

prob = ODEProblem(newton, rv0, (0.0u"s", Δt))
sol = solve(prob, Vern8(),dt=100u"s",adaptive=false)

```

It’s probably inefficient, but I wasn’t able to get it to use the function signature that modifies `du`.

---

<div class="post-metadata">

**Author:** ![jonniedie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonniedie/32/12842_2.png) [@jonniedie](https://discourse.julialang.org/u/jonniedie)\
**Post date:** [July 13, 2020, 8:34pm UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/7 "2020-07-13T20:34:38Z")

</div>

This also works now with ComponentArrays:

```julia
using ComponentArrays
using OrdinaryDiffEq
using LinearAlgebra
using Unitful

function newton(du, u, p, t)
    mu = 398600.4418u"km^3/s^2"
    r = norm(u.r)
    du.r = u.v
    du.v = -mu .* u.r / r^3
end

r0 = [1131.340, -2282.343, 6672.423]u"km"
v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
Δt = 86400.0*365u"s"
rv0 = ComponentArray(r=r0, v=v0)

prob = ODEProblem(newton, rv0, (0.0u"s", Δt))
sol = solve(prob, Vern8(), dt=100u"s", adaptive=false)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 14, 2020, 3:16am UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/8 "2020-07-14T03:16:18Z")

</div>

This is the recommended solution.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 14, 2020, 3:16am UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/9 "2020-07-14T03:16:38Z")

</div>

I think we have an open issue on this. Can you post in there and we’ll close it?

---

<div class="post-metadata">

**Author:** ![jonniedie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonniedie/32/12842_2.png) [@jonniedie](https://discourse.julialang.org/u/jonniedie)\
**Post date:** [July 14, 2020, 1:22pm UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/10 "2020-07-14T13:22:26Z")

</div>

[Is this the one?](https://github.com/SciML/DifferentialEquations.jl/issues/393)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 14, 2020, 2:05pm UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/11 "2020-07-14T14:05:14Z")

</div>

yes

---

<div class="post-metadata">

**Author:** ![pengwyn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pengwyn/32/6961_2.png) [@pengwyn](https://discourse.julialang.org/u/pengwyn)\
**Post date:** [July 15, 2020, 7:54am UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/12 "2020-07-15T07:54:39Z")

</div>

It should be pointed out that this relies on a bit of type twiddling underneath the hood. The `ComponentArray` type in this case will not have a concrete type for any of the elements.

The issue that I came across that led me here was that the algorithms in `OrdinaryDiffEq` rely on `zero(rate_prototype)`. Surprisingly, `zero(rv0)` returns a ComponentVector with no units, because it resolves to `similar.(rv0) .= 0`. (To be clear, it doesn’t return elements of type `Float64`, but rather elements with a `Unitful.Quantity` type with `NoUnits`)

So this should work, but will likely suffer from a speed penalty. Just pointing this out in case others get stuck in the loop I was in.

---

<div class="post-metadata">

**Author:** ![jonniedie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonniedie/32/12842_2.png) [@jonniedie](https://discourse.julialang.org/u/jonniedie)\
**Post date:** [July 15, 2020, 11:59am UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/13 "2020-07-15T11:59:13Z")

</div>

> [@pengwyn](#):
>
> Surprisingly, `zero(rv0)` returns a ComponentVector with no units.

Huh, that shouldn’t be the case. I’ll fix that.

---

<div class="post-metadata">

**Author:** ![jonniedie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonniedie/32/12842_2.png) [@jonniedie](https://discourse.julialang.org/u/jonniedie)\
**Post date:** [July 15, 2020, 12:16pm UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/14 "2020-07-15T12:16:41Z")

</div>

Alright, it’s fixed now.

```julia
julia> rv0.r
3-element view(::Array{Quantity{Float64,D,U} where U where D,1}, 1:3) with eltype Quantity{Float64,D,U} where U where D:
   1131.34 km
 -2282.343 km
  6672.423 km

julia> rv0.v
3-element view(::Array{Quantity{Float64,D,U} where U where D,1}, 4:6) with eltype Quantity{Float64,D,U} where U where D:
 -5.64305 km s^-1
  4.30333 km s^-1
  2.42879 km s^-1

```

Now, there is still going to be some speed penalty due to the fact that the `ComponentArray` is wrapping a heterogeneous array. @ChrisRackauckas, this might be a use case where `ArrayPartition`s are the better choice.

```julia
ap_rv0 = ArrayPartition(r0, v0)

@btime $rv0 + $rv0 # 6.360 μs (46 allocations: 1.78 KiB)
@btime $ap_rv0 + $ap_rv0 # 67.111 ns (4 allocations: 272 bytes)

```

But weirdly enough, ComponentArrays outperform their underlying array…

```julia
julia> a_rv0 = getdata(rv0)
6-element Array{Quantity{Float64,D,U} where U where D,1}:
       1131.34 km
     -2282.343 km
      6672.423 km
 -5.64305 km s^-1
  4.30333 km s^-1
  2.42879 km s^-1

julia> @btime $a_rv0 + $a_rv0;
  14.499 μs (71 allocations: 3.81 KiB)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 15, 2020, 1:16pm UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/15 "2020-07-15T13:16:39Z")

</div>

> [@jonniedie](#):
>
> But weirdly enough, ComponentArrays outperform their underlying array…

Probably the broadcast overload + constant prop is helping

---

<div class="post-metadata">

**Author:** ![pengwyn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pengwyn/32/6961_2.png) [@pengwyn](https://discourse.julialang.org/u/pengwyn)\
**Post date:** [July 15, 2020, 11:12pm UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/16 "2020-07-15T23:12:10Z")

</div>

> [@jonniedie](#):
>
> Alright, it’s fixed now.

Nice!

---

<div class="post-metadata">

**Author:** ![pengwyn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pengwyn/32/6961_2.png) [@pengwyn](https://discourse.julialang.org/u/pengwyn)\
**Post date:** [July 15, 2020, 11:30pm UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/17 "2020-07-15T23:30:53Z")

</div>

I think the underlying issue here is how to play nice with DifferentialEquations.jl underlying optimisations in preparing cache arrays, etc…

I think this is a discussion for a different thread, but I still want to mention what I’ve tried for this specific situation. The possibilities I’ve fiound are: (I’ll use `Q{D}` as a shorthand for `Unitful.Quantity{Float64,D,U}` which is not concrete)

1. Plain Vector with non-concrete type `Q{D}` to represent quantities with different units. This breaks because `zero(Q{D})` is (rightly) not defined.

2. ComponentArray with non-concrete parameter `Q{D}`, but each component itself is of a definite type. Works fine and produces correct results, but is slow because the derivative function receives an object with non-concrete type.

3. Hand-rolling my own struct:

```julia
struct Object{Tr,Tv}
r::Tr 
v::Tv 
end

```

I don’t think this works, because it is really difficult to make the internal `DifferentialEquations.jl` variables the correct type for units (e.g. handling the types returned by `recursive_unitless_bottom_eltype`, etc) in general. I think this approach would be fine if units were not an issue. (But then I would recommend `ComponentArrays` instead)

1. Mimicking a vector with:

```julia
struct Object{Tr,Tv} <: AbstractVector
r::Tr 
v::Tv 
end

```

This doesn’t work as it is invalid inheritance. The AbstractVector needs an element type and it is impossible to put a concrete type there.

My conclusion is that currently `DifferentialEquations.jl` requires a **single** element type, but it can happily work with a non-concrete type, it will just be slow.

Perhaps the frustrating bit here is that it isn’t possible to override the internal `DifferentialEquations.jl` choices for types of derivatives. Would there be anything wrong if I were to attempt a pull request that extracts out those lines of code to functions like `rate_type(u,t)`? This might also reduce the size of functions like `__init`.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 16, 2020, 1:31am UTC](https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/18 "2020-07-16T01:31:08Z")

</div>

> [@pengwyn](#):
>
> Would there be anything wrong if I were to attempt a pull request that extracts out those lines of code to functions like `rate_type(u,t)` ?

that would be helpful.

> [@pengwyn](#):
>
> My conclusion is that currently `DifferentialEquations.jl` requires a **single** element type, but it can happily work with a non-concrete type, it will just be slow.

More correctly, Julia in general requires strict typing on containers to be fast. I don’t think that holding units in the type domain should be the final answer, and instead it should be `isbits` but (value,symbol) pairs IMO. But that’s a bigger story.