Here is the way I found to do it but that's quiet heavy.
Plots it's so easy, but I don't convert ODE properly in julia.
const g = 9.8 # Accelaration of gravity
const p = 1.2 # Density of air
# Caracteristics of the problem
const m = 0.100 # A 100 g ball
const r = 0.10 # 10 cm radius
const Cd = 0.5 # Drag coeficient for a small spherical object
const y0 = 1000.0 # Initial height of the body (1000 m)
const v0 = 10.0 # Initial velocity of the body (10 m/s^2, going up)
const A = pi*r^2; # Cross-section area of the body;
function gm(t, f)
(y, v) = f # Extract y and v (i.e., dy/dt) from the f mapping
dy_dt = v # The differential equations
dv_dt = -1.0*g - sign(v)*(1./2.)*(p/m)*Cd*A*v^2.0
[dy_dt; dv_dt] # Return the derivatives
end;
# Initial conditions (position and velocity)
const start = [y0; v0]
# Time span (from 0 to 5 secs)
ts = [0.0; 5.0];
t, res = ode45(gm, start, ts)
y = map(x -> x[1], res)
v = map(x -> x[2], res);
using PyPlot
fig, ax = subplots(1, 2, sharex=true,
figsize=(10,4), dpi=80)
ax[1][:plot](t, v)
ax[1][:set_axis_bgcolor]("red")
ax[1][:set_axis_bgcolor]((0.5, 0.5, 0.5))
ax[1][:set_title]("Velocity over time");
ax[1][:set_xlabel]("Time (sec)")
ax[1][:set_ylabel]("Velocity (m/sec)")
ax[2][:plot](t, y)
ax[2][:set_title]("Height over time");
ax[2][:set_xlabel]("Time (sec)")
ax[2][:set_ylabel]("Height (m)");
Le vendredi 24 juin 2016 12:56:42 UTC+2, Henri Girard a écrit :