Table des matières

1 Etude du mouvement du pendule

In [1]:

%matplotlib inline
%autosave 300
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 14
from IPython.core.display import HTML,display
#from JSAnimation import IPython_display
css_file = 'style.css'
HTML(open(css_file, "r").read())

Autosaving every 300 seconds

Out[1]:

Etude du mouvement du pendule¶

Marc BUFFAT, dpt mécanique, Université Claude Bernard Lyon 1

Mise à disposition selon les termes de la Licence Creative Commons
Attribution - Pas d’Utilisation Commerciale - Partage dans les Mêmes Conditions 2.0 France.

Mouvement du pendule simple¶

In [2]:

%run ./draw_pendule.py

modèle mécanique¶

Le principe fondamental de la dynamique s’écrit:

$$ \frac{d}{dt}(m l \frac{d\theta}{dt}\overrightarrow{e_t})= -m g \overrightarrow{e_y} + T \overrightarrow{e_n} $$

dont la projection suivant $\overrightarrow{e_t}$ donne l’équation d’évolution de l’angle $\theta(t)$

$$ m l \ddot{\theta} = - m g \sin\theta $$

Les conditions initiales associées sont: $$ \theta(0)=\theta_0, \dot{\theta}(0) = 0 $$

propriété du système mécanique¶

La tension dans le fil étant perpendiculaire à la vitesse, son travail est nul. La gravité découle d’un potentiel, et donc le système est conservatif: $$ E_t = \frac{1}{2} m v^2 + m g y = cste$$ L’énergie totale se conserve au cours du temps

solution analytique¶

Dans le cas de petites oscillations: $\theta_0 \ll 1$: $$ \ddot{\theta} = -\omega^2 \theta \mbox{ avec } \omega=\sqrt{\frac{g}{l}}$$ dont la solution est périodique de la forme: $$\theta(t) = \theta_0\cos(\omega t)$$

Visualisation du mouvement et des forces¶

In [3]:

from anim_pendule import anim_pendule
anim=anim_pendule()

In [4]:

HTML(anim.to_html5_video())

Out[4]:

Objectif¶

Etude de la période d’oscillation en fonction de $\theta_0$

Analyse descendante¶

Pour chaque valeur de $\theta_0$ de $\epsilon=0.001$ radians à $\frac{pi}{2}$

calcul de la solution $\theta(t)$
détermination de la période $T$

calcul de la solution $\theta(t)$¶

transformation en un système d’EDO du 1er ordre $$ \dot{Y} = F(Y,t) $$
- définition de $Y(t)$:
$$\begin{eqnarray} Y(t) &=& [Y[0],Y[1]]\ &=& [\theta(t),\dot{\theta}(t)] \end{eqnarray}$$
- définition du second membre $F(Y,t)$:
$$ \begin{eqnarray} F(Y,t) &=& [\dot{Y}[0],\dot{Y}[1]]\ &=& [ Y[1], -\omega_0^2 \sin(Y[0] )) \end{eqnarray}$$
- définition de la condition initiale $Y(0)$
$$ Y(0) = [ \theta_0 , 0 ] $$
utilisation du solveur d’EDO de la bibliothèque numpy/scipy
- odeint https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

calcul de la période¶

a partir de la courbe $\theta(t)$ donnée par les coupes de points $\{t_i, \theta_i\}_{i=0,n}$, on détermine les 2 premiers passages par 0 ce qui donne la demi période.
détermination du passage par zero: on détermine le premier point k tel que $\theta_k<0$ et donc le passage par zero se fait entre les points k-1 et k.
interpolation linéaire avec les points k-1 et k pour déterminer le passage par zero.

In [5]:

# parametres
g = 10.
L = 1.
omega0 = np.sqrt(g/L)
# fonctions de base
def F(Y,t):
    """second membre de l'EDO"""
    global omega0
    dY = np.array([Y[1],-omega0**2*np.sin(Y[0])])
    return dY

def Energie(Y):
    """Calcul l'energie totale du pendule en fct de la solution Y"""
    E = 0.5*L**2*Y[:,1]**2 + g*L*(1-np.cos(Y[:,0]))
    return E

In [6]:

# integration avec une CI theta0
from scipy.integrate import odeint
def solution(theta0):
    T = np.linspace(0.,2*(2*np.pi)/omega0,101)
    # condition initiale
    Y0 = np.array([theta0,0.])
    # integration
    Y = odeint(F, Y0, T)
    return T,Y

In [7]:

# tracer de la solution  et vérification conservation energie
# sauveargarde de l'image si image_png est définie)
def trace(T,Y,image_png=None):
    Theta=Y[:,0]
    Thetalin=Theta[0]*np.cos(omega0*T)
    
    plt.figure(figsize=(14,7))
    # trace theta(t)
    plt.subplot(1,2,1)
    plt.plot(T,Theta,'x',label="$\\theta$",lw=2)
    plt.plot(T,Thetalin,label="$\\theta_{lin}$",lw=2)
    plt.title('Angle $\\theta(t)$')
    plt.xlabel('t')
    plt.legend()
    # et l'energie Et
    plt.subplot(1,2,2)
    Et=Energie(Y)
    plt.plot(T,(Et-Et[0])/Et[0],lw=2)
    plt.ylim(-1.e-5,1.e-5)
    plt.xlabel('t')
    plt.title("Erreur relative sur Et")
    if image_png!=None: plt.savefig(image_png)
    return

In [8]:

# calcul de la période de la solution numérique 
# determination des 2 intervalles k1-1,k1  et k2-1,k2
# correspondant au passage par zero de theta
def periode(T,Theta):
    k1=0
    while Theta[k1] > 0 : k1=k1+1
    k2=k1
    while Theta[k2]<0: k2=k2+1
    # interpolation lineaire
    t1 = T[k1-1] - (T[k1]-T[k1-1])*Theta[k1-1]/(Theta[k1]-Theta[k1-1])
    t2 = T[k2-1] - (T[k2]-T[k2-1])*Theta[k2-1]/(Theta[k2]-Theta[k2-1])
    # d'ou la periode
    return 2*(t2-t1)

validation cas petites oscillations¶

In [9]:

# verification dans le cas lineaire
T,Y=solution(0.01)
per =periode(T,Y[:,0])
print("periode lineaire ",per," err:",per-2*np.pi/omega0)
trace(T,Y,"caslin.png")

periode lineaire  1.9869313449  err: 1.36917407347e-05

calcul en grandes oscillations¶

In [10]:

T,Y = solution(np.pi/2)
per =periode(T,Y[:,0])
print("periode lineaire ",per," err:",per-2*np.pi/omega0)
trace(T,Y,"casnonlin.png")

periode lineaire  2.34527849973  err: 0.358360846573

Evolution de la période en fonction de $\theta_0$¶

comparaison avec le développement asymptotique de Borda

$$T = 2\pi \sqrt{\frac{l}{g}} \left(1+(\theta_0)^{\frac{2}{16}}\right)$$

In [11]:

# analyse parametrique de 0.01 a pi/2
Per0 = 2*np.pi/omega0
N  = 21
Theta0 = np.linspace(0.01,np.pi/2,N)
Per    = np.zeros(N)

for k,theta0 in enumerate(Theta0):
    # integration
    T,Y = solution(theta0)
    # periode
    Per[k] = periode(T,Y[:,0])

In [12]:

# interpolation lagrange avec 1 pt sur 5
from scipy.interpolate import lagrange

print("periode T0:",Per0)
Pi = lagrange(Theta0[0::5],Per[0::5])
print("polynome interpolation:\n", Pi)

periode T0: 1.98691765316
polynome interpolation:
          4          3          2
0.01366 x - 0.0126 x + 0.1324 x - 0.00176 x + 1.987

In [13]:

# tracer de la periode
plt.figure(figsize=(10,8))
plt.plot(Theta0,Per,'x',label="T non lin.",ms=18)
plt.plot(Theta0,Pi(Theta0),label="interp.",lw=2)
plt.plot(Theta0,Per0*(1+Theta0**2/16.),label="Borda",lw=2)
plt.plot(Theta0,Per0*np.ones(N),'--',label="T0")
plt.legend()
plt.title("Periode fct de $\\theta_0$")
plt.xlabel("$\\theta_0$")
plt.savefig("periode.png")

Fin¶

In [ ]:

Cours Python scientifique