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.