Error using Differential Equations and Unitful

Specific Domains Numerics
, ,
1

I am using Differential Equations and Unitful together for the first time and I am running into the following error when trying to solve my OdeProblem: ArgumentError: zero(Unitful.Quantity{Float64}) not defined.

My ODE is pretty simple and everything has units:

function rocket!(du, u, p, t)
	ρ = ρ_altitude(u[1]) # kg/m^3
	drag = Fd(u[2], ρ) # N

	F = thrust - drag # N

	du[1] = u[2] # x' = ẋ
	du[2] = F / 300000 # ẍ = (thrust - drag) / mass	
end

The problem works great when I remove the units, I have confirmed that all of the units are valid, and all packages are fully up to date so I am not sure why the argument error is popping up.

If there is any more information needed let me know and I will add it to the post, but the full code is in a comment below: Error using Differential Equations and Unitful - #3 by MisterBiggs

2

Quantity{Float64} is not a completely specified type.

julia> typeof(Quantity{Float64})
UnionAll

help?> UnionAll

  UnionAll

 A union of types over all values of a type parameter. UnionAll is used to describe parametric types where the values of some parameters are not
  known.

Furthermore,

julia> isconcretetype(Quantity{Float64})
false

julia> isconcretetype(typeof(Quantity(5, u"m")))
true

So, it doesn’t make sense to define the zero (additive identity) for an incomplete type. Well, you could make some kind of incompletely specified zero. But, this is not done anywhere in Julia that I know of. And, at least in one case, an assumption is made and a type is chosen, Int.

julia> typeof(Rational)
UnionAll

julia> zero(Rational)
0//1

julia> typeof(zero(Rational))
Rational{Int64}

But, I’d argue this is poorly designed. We just had a long thread, where a few people, most notably @StefanKarpinski opined that a major strength of Julia (as opposed to R) is that it refuses to guess what the user might mean when the intent is not fully specified. It’s true in general, but this may have slipped through. It might have been removed if an issue had been raised prior to 1.0. Also, this is not nearly so surprising as the examples given for R.

Returning to the problem at hand.
This does make sense:

julia> typeof(Quantity(5, u"m"))
Quantity{Int64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}

julia> zero(typeof(Quantity(5, u"m")))
0 m

I guess that somewhere in your code the units have not been specified. I can’t tell from your supplied code snippet.

2 Likes
3

I really appreciate the detailed reply. I wasn’t aware of the concept of concrete types but that seems to be a very likely source of the issue. I don’t think a missing unit is the issue because I can strip the units from ever appearing in the ODEProblem and everything runs perfectly:

function rocket!(du, u, p, t)
	ρ = ρ_altitude(u[1]) # calculate air density at current altitude
	drag = Fd(u[2]*u"m/s", ρ) # calculate drag force at current air density

	F = thrust - drag

	du[1] = u[2] 
	du[2] = F / mass |> x -> ustrip(u"m/s^2", x) # ẍ = (thrust - drag) / mass	
end

Here is the full code with all the units in place that isn’t currently working:

using Unitful
using DifferentialEquations

function ρ_altitude(height)
	"""
    https://www.grc.nasa.gov/www/k-12/airplane/atmosmet.html
	Calculate ρ as a function of altitude
	"""

	h = ustrip(u"m", height) # The units on the following equations don't make sense
	
	if h < 11000
		T = 15.04 - .00649h
		P = 101.29 * ((T + 273.1) / 288.08)^5.256
	elseif h < 25000
		T = -56.46
		P = 22.65 * exp(1.73 - .000157h)
		
	else
		T = -131.21 + .00299 * h 
		P = 2.488 * ((T + 273.1) / 216.6)^-11.388
	end
	
	P = P
	T = (T + 273.15)
	
	return P / .2869T * u"kg/m^3"
end

d = 2.5u"m" # max diameter Soyuz 7k wikipedia

A = π * d^2 / 4

Cd = .82

Fd(v,ρ) = Cd * .5 * ρ * A * v^2  |> u"N"

u0 = [0.01u"m", 0.01u"m/s"]

tspan = (0.01u"s", 118.0u"s") # Time for first stage

thrust = 846000u"lbf" |> u"N"

function rocket!(du, u, p, t)
	ρ = ρ_altitude(u[1])
	D = Fd(u[2], ρ)

	F = thrust - D

	du[1] = u[2] # x' = ẋ
	du[2] = F / 300000u"kg" # ẍ = (thrust - drag) / mass	
end

prob = ODEProblem(rocket!, u0, tspan)

sol = solve(prob)

I’ve tried setting the initial conditions to something non-zero to see if that would help but it has not I am still getting the ArgumentError: zero(Unitful.Quantity{Float64}) not defined from my original post.

4

You may be making a mistake, but I don’t see it. It’s likely a bug. Since you are not getting an answer here quickly, I suggest opening an issue on github.

1 Like
5

Printing u0 and tspan like so

println("u0 ", u0)
println("tspan ", tspan)
prob = ODEProblem(rocket!, u0, tspan)

shows

u0 Quantity{Float64, D, U} where {D, U}[0.01 m, 0.01 m s^-1]
tspan (0.01 s, 118.0 s)

Are you sure u0 is correctly defined?

6

It seems that multiple Unitful units are a known issue in DifferentialEquations.jl https://github.com/SciML/DifferentialEquations.jl/issues/352

Glad to see it wasn’t something I was doing that was the issue but also disappointing that this use case is not yet possible. I’ll post my workaround code just in case anyone else stumbles on this thread in the future.

using Plots
using DifferentialEquations
using Unitful
using Plots


function ρ_altitude(h)
	"""https://www.grc.nasa.gov/www/k-12/airplane/atmosmet.html
	Calculate ρ as a function of altitude
	"""
	
	# Units can be ignored in this function since these equations are just a curve fit
	
	if h < 11000
		T = 15.04 - .00649h
		P = 101.29 * ((T + 273.1) / 288.08)^5.256
	elseif h < 25000
		T = -56.46
		P = 22.65 * exp(1.73 - .000157h)
		
	else
		T = -131.21 + .00299 * h 
		P = 2.488 * ((T + 273.1) / 216.6)^-11.388
	end
	
	P = P
	T = (T + 273.15)
	
	return P / .2869T * u"kg/m^3"
end

d = 2.62u"m"
A = π * d^2 / 4

Cd = .82

Fd(v,ρ) = Cd * .5 * ρ * A * v^2 |> u"N"

u0 = [0.0, 0.0]
tspan = (0.0, 118.0) # Time for first stage

thrust = 760600u"N"
mass  = 39372u"kg"

function rocket!(du, u, p, t)
	ρ = ρ_altitude(u[1])
	drag = Fd(u[2] * u"m/s", ρ)

	F = thrust - drag

	du[1] = u[2] 

	# ẍ = (thrust - drag) / mass	
	du[2] = F / mass |> x -> ustrip(u"m/s^2", x)
	#                        ^^^^^^
	# Strip units before returning to the ode function 
end


prob = ODEProblem(rocket!, u0, tspan, dtmax=.5)
sol = solve(prob)


plot(sol, layout=(2, 1), label=nothing,xlabel=["" "Time (s)"], ylabel=["Altitude (m)" "Velocity (m/s)"])