In [1]:
import numpy as np
import matplotlib.pylab as plt

%matplotlib inline

Calculons le lieu d'impact du boulet¶

Soit un bon vieux canon qui tire avec un angle $\alpha$ un boulet de poids $m$ avec une force boum pendant un temps $t_b$.

  • Dessiner la trajectoire du boulet avec une approximation par la méthode des différences finies à l'ordre 1 et à l'ordre 2.
  • Calculer le point d'impact exact.
  • Ajouter le frottement de l'air en supposant que le frottement est une force égale à -k fois la vitesse du boulet (k=0.05 est bien).
In [2]:
# votre long programme ici :

Solution¶

A tout moment $t$, la position du boulet est $(x,y)(t)$, sa vitesse est $\vec{v}(t) = (v_x,v_y)(t)$ et son accélération est $\vec{a}(t) = (a_x,a_y)(t)$.

L'équation fondamentale ici est $\vec{f} = m \; \vec{a} $.

Il y a deux forces, celle qui propulse le boulet hors du canon et la gravité qui fait retomber le boulet sur terre.

  • la force d'explosion boum dure le temps que le boulet est dans le canon. Supposons que cela soit 0,1 seconde alors la vitesse de sortie du boulet est $\vec{v_0}$ qui vient de l'intégration de notre équation fondamentale en 0 et 0,1 seconde : $$ \int_0^{t_b} \mathtt{boum} \; dt = m \; (\vec{v_0} - 0) = m\; \vec{v_0}$$ soit $$\vec{v_0} = \vec{v}(0) = \frac{t_b}{m} \mathtt{boum} $$ avec la direction du vecteur étant celle du canon (et celle de boum).

  • la gravité s'applique tout le temps où le boulet est en l'air sur la composante verticale de la vitesse : $$a_y = g = - 9.81 \quad \forall \; t \; / \; y > 0$$

Pour calculer la trajectoire avec la méthode des différences finies, on utilise le fait que la position est la dérivée de la vitesse, ainsi : $$ \vec{v}(t) = \frac{\partial{(x,y)}}{\partial{t}}(t)$$ On approche la dérivée à l'aide des développements limité :

$$ x(t+\epsilon) = x(t) + \epsilon \; x'(t) + \frac{\epsilon^2}{2} x''(t) + ...$$

donc en introduisant la vitesse et l'accélération :

$$x(t+\epsilon) = x(t) + \epsilon \; v_x(t) + \frac{\epsilon^2}{2} a_x(t) + ...$$

$$y(t+\epsilon) = y(t) + \epsilon \; v_y(t) + \frac{\epsilon^2}{2} a_y(t) + ...$$

On peut ainsi calculer la position en $t+\epsilon$ connaissant la position en $t$. Bien sûr si on "oublie" des termes, comme l'accélération on risque d'avoir une erreur. C'est ce qui se passe lorsque qu'on ne fait le développement limité qu'à l'ordre 1.

Notons aussi que l'on aurait pu faire le développement limité en $x(t-\epsilon)$.

Enfin on a vu que la seule accélération qui s'applique une fois le boulet parti est la gravité, donc $(a_x,a_y) = (0,-9,81)\;$. Donc en prenant en compte cette accélération le résultat sera meilleur (approximation à l'ordre 2 ci-dessous).

Finalement, comme il n'y a pas de terme de niveau 3 ou plus qui s'applique ici, le jerk est nul donc les $...$ valent 0 et donc l'approximation d'ordre 2 donne le résultat exact. On note en effet que la courbe d'ordre 2 touche $y=0$ au point d'impact calculé analytiquement.

Calcul du point d'impact

$$ \exists \; T > 0 \; / \; y(T) = 0$$ $$ \int_0^T v_y(t)\; dt = 0 - 0 = 0$$

$ \displaystyle \int_0^T v_y(0) - 9,81 \; t \quad dt = 0 \quad$ avec le calcul précédent de $v_y(t)$

$$ v_y(0) \; T - 9,81 \; \frac{T^2}{2} = 0 $$ donc

$\displaystyle T = \frac{2\; v_y(0)}{9,81} = \frac{2 \; v_0 \sin(\alpha)}{9.81}\quad $ avec $\alpha$ l'angle de tir.

On en conclut la position de chute du boulet : $x(T) = 0 + T \; v_x = T \; v_0 \; cos(\alpha)$.

Frottements

L'ajout des frottements d'air, tel que suggéré, revient à appliquer l'accélération $\displaystyle \vec{a} = \frac{-k\; \vec{v}}{m}$

In [3]:
boum = 200
m = 2
t_boum = 0.1
a0 = boum/m
v0 = a0 * t_boum # v0 est a vitesse à la fin de t_boum
alpha = np.pi/4     # angle de tir
print("Vitesse initiale : %4.2f m/s\nAngle de tir     : %4.2f°" % (v0,alpha*180/np.pi))
delta = 0.1      # pas de temps de la discrétisation
                 #   passer à 0.01 pour voir la différence
plt.figure(figsize=(13,3), dpi=100)
plt.axis([0,13,0,3])

# approximations au degre 1
# (on utilise la vitesse du pas de temps suivant)
vx=[v0*np.cos(alpha)]
vy=[v0*np.sin(alpha)]
x=[0]
y=[0]
for _ in np.arange(delta,2,delta):
    vx.append(vx[-1])                    # [-1] est le dernier élément de la liste qui grandit
    vy.append(vy[-1]-delta*9.81)
    x.append(x[-1]+delta*vx[-1])
    y.append(y[-1]+delta*vy[-1])
plt.plot(x,y, label = "1 suiv")

# approximations au degre 1
# (on utilise la vitesse du pas de temps précédant
vx=[v0*np.cos(alpha)]
vy=[v0*np.sin(alpha)]
x=[0]
y=[0]
for _ in np.arange(delta,2,delta):
    x.append(x[-1]+delta*vx[-1])
    y.append(y[-1]+delta*vy[-1])
    vx.append(vx[-1])
    vy.append(vy[-1]-delta*9.81)
plt.plot(x,y, label="1 prec")

# approximations au degre 2
vx=[v0*np.cos(alpha)]
vy=[v0*np.sin(alpha)]
x=[0]
y=[0]
for _ in np.arange(delta,2,delta):
    x.append(x[-1] + delta*vx[-1])
    y.append(y[-1] + delta*vy[-1] - 0.5 * delta**2 * 9.81)
    vx.append(vx[-1])
    vy.append(vy[-1] - delta*9.81)
plt.plot(x,y, label="ordre 2")

# approximations au degre 2 avec frottements
k=0.05 # coef de frottement
vx=[v0*np.cos(alpha)]
vy=[v0*np.sin(alpha)]
x=[0]
y=[0]
for _ in np.arange(delta,2,delta):
    x.append(x[-1] + delta*vx[-1] - 0.5 * delta**2 * k * vx[-1] / m)
    y.append(y[-1] + delta*vy[-1] - 0.5 * delta**2 * (9.81 + k * vy[-1]/m))
    vx.append(vx[-1] - delta * k * vx[-1] / m)
    vy.append(vy[-1] - delta*(9.81 + k*vy[-1]/m))
plt.plot(x,y, label="2 + air")

# calcul exact
# y = 0 au temps t 
t = 2 * v0 * np.sin(alpha) / 9.81
xx = t * v0 * np.cos(alpha)
print("Calcul exact sans frottement : le projectile retombe en x = %5.2f" % xx)
plt.plot([xx],[0],'ko')
plt.legend()
    

Vitesse initiale : 10.00 m/s
Angle de tir     : 45.00°
Calcul exact sans frottement : le projectile retombe en x = 10.19
Out[3]:
<matplotlib.legend.Legend at 0x7fe3b81b82e0>
No description has been provided for this image

On note qu'en baissant la masse les frottements d'air sont plus importants.

Optimisons¶

Nous avons travaillé avec des listes ce qui n'est pas très performant. Il est préférable d'utiliser les tableaux de Numpy mais cela doit être fait correctement...

Voici 3 fonctions qui calculent la même chose,

  • la première avec des listes
  • la seconde avec les tableaux mais naïvement
  • la troisième avec les tableaux mieux utilisés

On va comparer leur vitesse d'exécution avec la commande magique %timeit de iPython (commande qui appelle simplement timeit de Python).

In [4]:
def calcul_liste(delta):
    vx=[v0*np.cos(alpha)]
    vy=[v0*np.sin(alpha)]
    x=[0]
    y=[0]
    for _ in np.arange(delta,2,delta):
        vx.append(vx[-1])
        vy.append(vy[-1]-delta*9.81)
        x.append(x[-1]+delta*vx[-1])
        y.append(y[-1]+delta*vy[-1])
    return (x,y)

%timeit calcul_liste(0.01)
(x1,y1) = calcul_liste(0.1)

130 µs ± 3.68 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [5]:
def calcul_array_naive(delta):
    nb_iteration = round(2/delta)
    x = np.zeros(nb_iteration)
    y = np.zeros(nb_iteration)
    vx=np.empty(nb_iteration)
    vy=np.empty(nb_iteration)
    vx[0]=v0*np.cos(alpha)
    vy[0]=v0*np.sin(alpha)
    for i in np.arange(1,nb_iteration):
        vx[i] = vx[i-1]
        vy[i] = vy[i-1]-delta*9.81
        x[i] = x[i-1]+delta*vx[i-1]
        y[i] = y[i-1]+delta*vy[i-1]
    return(x,y)

%timeit calcul_array_naive(0.01)

415 µs ± 22.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [6]:
def calcul_array(delta):
    nb_iteration = round(2/delta)
    vx=np.empty(nb_iteration)
    vy=np.empty(nb_iteration)
    vx.fill(v0*np.cos(alpha))
    vy.fill(-delta*9.81)
    vy[0] = v0*np.sin(alpha)
    vy = np.cumsum(vy)   # somme cumulative des éléments de vy 
    x = delta*vx
    x[0] = 0
    x = np.cumsum(x)
    y = delta*vy
    y[0] = 0
    y = np.cumsum(y)
    return (x,y)
    
%timeit calcul_array(0.01)
(x2,y2) = calcul_array(0.1)
print(x1-x2)  # pour vérifier qu'on a bien le même résultat
print(y1-y2)

14.7 µs ± 121 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

On constate que

  • la première méthode est 3 fois plus rapide que la deuxième qui utilise pourtant les tableaux, oups !
  • la méthode naïve est près de 28 fois plus lente que la méthode avec les cumsum
  • la meilleure méthode est 9 fois plus rapide que celle avec les listes

Moralité : l'optimisation a du bon

Notons enfin, que la dernière méthode est parallélisable (la parallélisation de cumsum est celle du problème dit du préfixe).

In [ ]: