I am new to Juia lang and trying to solve the following differential equations to find the terminal velocity of a ball using Julia.
F = - m * g - 1/2 rho * v² Cd * A
This is the code that I wrote:
# Termal velocity of a falling ball
using DifferentialEquations
using Plots
g = 9.8 # Accelaration of gravity
p = 1.2 # Density of air
m = 0.100 # A 100 g ball
r = 0.10 # 10 cm radius
Cd = 0.5 # Drag coeficient for a small spherical object
y0 = 1000.0 # Initial height of the body (1000 m)
v0 = 10.0 # Initial velocity of the body (10 m/s^2, going up)
A = pi*r^2; # Cross-section area of the body;
u0 = [v0;y0] # Initial Conditions
tspan = (0.0,5.0) # Time span to solve for
p = [g;p;m;Cd;A]
function Terminal_Velocity(du,u,p,t)
du[1] = u[1] # velocity
du[2] = -1.0 * p[1] - 0.5 * (p[2]/p[3]) * (u[1]^2) * p[4] * p[5] # acceleration
end
prob = ODEProblem(Terminal_Velocity,u0,tspan,p)
sol = solve(prob)
plot(sol,vars=(0,1))
I think the Problem is that I am giving y0 as the initial condition for the acceleration and not for the height. But I can do not understand the Syntax well enough yet.
My starting point was this article : https://nbviewer.jupyter.org/github/JuliaLang/ODE.jl/blob/master/examples/Terminal_Velocity.ipynb
Thanks for your help in advance.
There are several errors in your example. Most of them are not very related to programming but to physics and mathematics.
You are disregarding a sign-change in the drag term. Also, the drag term you specified in your F equation has an additional error (an extra 1/m).
You seem to be mixing up what velocity and acceleration is. du[2] is acceleration as it is the derivative of the velocity (u[2]). You are using u[1] as velocity.
du[1] = u[1] gives an exponential increase of u[1], what you want is du[1] = u[2] which is saying that the position affected by the velocity.
The order of u0 = [v0;y0] is flipped, u[1] is the y coordinate while u[2] is the velocity.
The only programming error that I can see is in using 0-based indexing when choosing which variables to plot.
After fixing the above points, you get:
using DifferentialEquations
using Plots
g = 9.8 # Accelaration of gravity
p = 1.2 # Density of air
m = 0.100 # A 100 g ball
r = 0.10 # 10 cm radius
Cd = 0.5 # Drag coeficient for a small spherical object
y0 = 1000.0 # Initial height of the body (1000 m)
v0 = 10.0 # Initial velocity of the body (10 m/s^2, going up)
A = pi*r^2; # Cross-section area of the body;
u0 = [y0;v0] # Initial Conditions
tspan = (0.0,5.0) # Time span to solve for
p = [g;p;m;Cd;A]
function Terminal_Velocity(du,u,p,t)
du[1] = u[2] # velocity
du[2] = - p[1] - sign(u[2]) * 0.5 * (p[2]/p[3]) * (u[2]^2) * p[4] * p[5] # acceleration
end
prob = ODEProblem(Terminal_Velocity,u0,tspan,p)
sol = solve(prob)
plt1 = plot(sol; vars=1)
plt2 = plot(sol; vars=2)
plot(plt1, plt2)
One could go even further and use callbacks to ensure that the sign-change does not cause numerical errors.
To do so, replace the solve line with
cond(u, t, i) = u[2]
callback = ContinuousCallback(cond, nothing)
sol = solve(prob; callback=callback)