2

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.

2 Answers2

2

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)
Korsbo
  • 724
  • 5
  • 5
  • I copied you solution and it gave back the following Error: UndefVarError: callback not defined Stacktrace: [1] top-level scope at In[62]:23 –  Oct 03 '19 at 11:05
  • Yeah, there was an error but I have already edited it out, try again. – Korsbo Oct 03 '19 at 11:07
  • It is working now. But regarding the physics the 1/m term is correct since F=m* d²x/dt² . I think you should fix that in you answer! –  Oct 03 '19 at 11:47
  • The `1/m` factor should indeed exist in your `du[2]` equation since this is an acceleration. `F = - m * g - 1/2 rho/m * v² Cd * A` is, however, wrong. Have a look at the notebook you referred to. – Korsbo Oct 03 '19 at 11:59
  • That one does not have a problem with the masses, but the equation you type out before the code does. – Korsbo Oct 03 '19 at 12:06
  • I understood that now and fixed it in my initial question! –  Oct 03 '19 at 12:14
1

I think the main error came from the flipped sign in:

du[2] = -1.0 * p[1] - 0.5 * (p[2]/p[3]) * (u[1]^2) * p[4] * p[5]

Which should be:

du[2] = +1.0 * p[1] - 0.5 - sign(u[2]) * (p[2]/p[3]) * (u[1]^2) * p[4] * p[5]

However, p and rho are also easy to confuse because you reassign it when setting up the ODE parameters.

I changed the setup of the ODE slightly (i.e. u[1] is the displacement now). This should work:

# Termal velocity of a falling ball
using DifferentialEquations
using Plots

g  = 9.8                # Accelaration of gravity
rho  = 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 
u0 = [y0, v0]            # Initial Conditions
tspan = (0.0,5.0)       # Time span to solve for
p = [g rho m Cd A]

function Terminal_Velocity(du,u,p,t)
 (g, rho, m, Cd, A) = p
 du[1] = u[2]                                                      # velocity 
 du[2] = -g - 0.5 * sign(u[2]) * (rho/m) * (u[2]^2) * Cd * A  # acceleration
end



prob = ODEProblem(Terminal_Velocity,u0,tspan,p)
sol = solve(prob)

p1 = plot(sol, vars=(1), label="Displacement")
p2 = plot(sol, vars=(2), label="Velocity")

plot(p1, p2)

Edit: Fixed sign error.

Wolf
  • 485
  • 4
  • 8