CODEWITHSUNDEEP
juliascientific-computingdifferentialequations.jl
user12158459
user12158459

Reputation:

Terminal Velocity using Differential Equation

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.

Upvotes: 2

Views: 311

Answers (2)

Wolf
Wolf

Reputation: 505

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.

Upvotes: 1

Korsbo
Korsbo

Reputation: 724

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)

Upvotes: 2

Related Questions