Résolution numérique d'une équation différentielle du second ordre

Le script aborde la résolution numérique d'une équation différentielle du second ordre du type :

\begin{equation} \mathsf{ \frac{d^2u}{dt^2} = g(\frac{du}{dt}) + h(u) } \end{equation}

Bibliothèques utilisées

La bibliothèque numpy est utilisée pour sa collection de fonctions (sin par exemple) et pou générer une liste d'abscisses avec la fonction linspace.

La bibliothèque matplolib est utilisée pour tracer des graphiques.

Le module integrate de la bibliothèque scipypermet d'intégrer numériquement une équation différentielle.

Principe de résolution numérique

L'intégration de l'équation différentielle est ici réalisée grâce à la fonction odeint qui se trouve dans le module integratede la bibliothèque scipy. Lorsque l'équation est du second degré en une variable $\mathsf{u}$, la résolution nécessite de créer un vecteur contenant la grandeur $\mathsf{u}$ et sa dérivée première $\mathsf{\frac{du}{dt}}$ :

\begin{equation} \mathsf{ X = \begin{pmatrix} \u \\ \frac{du}{dt} \end{pmatrix} } \end{equation}

Il faut ensuite définir une fonction $\mathsf{f}$ qui contient la dérivée première du vecteur $\mathsf{X}$ :

\begin{equation} \mathsf{ f : \begin{pmatrix} \u \\ \frac{du}{dt} \end{pmatrix} \mapsto \begin{pmatrix} \frac{du}{dt} \\ \frac{d^2u}{{dt}^2} \end{pmatrix} } \end{equation}

Lors de la définition de la fonction $\mathsf{f}$, le terme de dérivée seconde est remplacé par l'expression issue de l'équation différentielle.

Exemple avec l'équation différentielle précédemment décrite : \begin{equation} \mathsf{ f : \begin{pmatrix} u \\ \frac{du}{dt} \end{pmatrix} \mapsto \begin{pmatrix} \frac{du}{dt} \\ g(\frac{du}{dt}) + h(u) \end{pmatrix} } \end{equation}

Trois arguments doivent être saisis dans la commande odeint pour que l'intégration soit réalisée :

odeint(f, X0 , t)

où :

  • X0désigne la liste des conditions initiales $\mathsf{{u}_0}$ et $\mathsf{({\frac{du}{dt}})_0}$
  • tdésigne la liste des instants t. Celle-ci peut créée grâce à la fonction linspace de la bibliothèque numpy sous la forme np.linspace(borne_inf,borne_sup,nombre_valeurs).

odeint renvoie un tableau dont la première colonne est la liste des valeurs successives de $\mathsf{u}$ et la seconde, celles des valeurs succesives de la dérivée première $\mathsf{\frac{du}{dt}}$.

Pour extraire une partie des données :

  • liste des valeurs de $\mathsf{u}$ : ajouter [:,0] (1ère colonne)
  • liste des valeurs de $\mathsf{\frac{du}{dt}}$ : ajouter [:,1] (2ème colonne)

Premier exemple : chute libre sans frottement

Soit l'équation différentielle :

\begin{equation} \mathsf{ \frac{d^2z}{dt^2} = -9.8 } \end{equation}

Dans cette situation, la fonction $\mathsf{f}$ doit être définie comme suit : \begin{equation} \mathsf{ f : \begin{pmatrix} z \\ zpoint \end{pmatrix} \mapsto \begin{pmatrix} zpoint \\ -9.8 \end{pmatrix} } \end{equation}

In [22]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Durée d'intégration
t=np.linspace(0,4,30)

# Fonctions utiles à la résolution
def f(X,t):
    z,zpoint = X
    return zpoint,-9.8

def alt(z0,zpoint0):
    return odeint(f,[z0,zpoint0],t)[:,0]
    # Retourne la liste des angles (1ère colonne)

def vit(z0,zpoint0):
    return odeint(f,[z0,zpoint0],t)[:,1]
    # Retourne la liste des vitesses (2ème colonne)
   
# Tracé
z0=50
zpoint0=0


plt.subplot(211)
plt.title(f'Chute libre \n Vitesse initiale = {zpoint0} m/s \n Hauteur initiale = {z0} m')
plt.plot(t,alt(z0,zpoint0),'.')
plt.xlabel('temps en s')
plt.ylabel('altitude en m')
plt.grid()
         
plt.subplot(212)
plt.plot(t,vit(z0,zpoint0),'.')
plt.xlabel('temps en s')
plt.ylabel('vitesse en m')
plt.grid()
         
plt.show()

Second exemple : chute libre avec frottements

Soit l'équation différentielle :

\begin{equation} \mathsf{ \frac{d^2z}{dt^2} = -9.8 - 0.9 \cdot \frac{dz}{dt} } \end{equation}

Dans cette situation, la fonction $\mathsf{f}$ doit être définie comme suit : \begin{equation} \mathsf{ f : \begin{pmatrix} z \\ zpoint \end{pmatrix} \mapsto \begin{pmatrix} zpoint \\ -9.8 - 0.9 \cdot zpoint \end{pmatrix} } \end{equation}

In [27]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Durée d'intégration
t=np.linspace(0,6,30)

# Fonctions utiles à la résolution
def f(X,t):
    z,zpoint = X
    return zpoint, -9.8-.9*zpoint

def alt(z0,zpoint0):
    return odeint(f,[z0,zpoint0],t)[:,0]
    # Retourne la liste des angles (1ère colonne)

def vit(z0,zpoint0):
    return odeint(f,[z0,zpoint0],t)[:,1]
    # Retourne la liste des vitesses (2ème colonne)
   
# Tracé
z0=50
zpoint0=0

plt.subplot(211)
plt.title(f'Chute libre avec frottements \n Vitesse initiale = {zpoint0} m/s \n Hauteur initiale = {z0} m')
plt.plot(t,alt(z0,zpoint0),'.')
plt.xlabel('temps en s')
plt.ylabel('altitude en m')
plt.grid()
         
plt.subplot(212)
plt.plot(t,vit(z0,zpoint0),'.')
plt.xlabel('temps en s')
plt.ylabel('vitesse en m')
plt.grid()
         
plt.show()

Troisième exemple : Régime libre d'un oscillateur harmonique amorti

Soit l'équation différentielle :

\begin{equation} \mathsf{ \frac{d^2z}{dt^2} = - 9 \cdot z - \frac{3}{Q} \cdot \frac{dz}{dt} } \end{equation}

Plusieurs valeurs de Q vont être testées.

La fonction $\mathsf{f}$ doit être définie comme suit : \begin{equation} \mathsf{ f : \begin{pmatrix} z \\ zpoint \end{pmatrix} \mapsto \begin{pmatrix} zpoint \\ -9 \cdot z - \frac{3}{Q} \cdot zpoint \end{pmatrix} } \end{equation}

In [31]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint


# Durée d'intégration
t=np.linspace(0,6,100)


# Fonctions utiles à la résolution
def f(X,t):
    z,zpoint = X
    return zpoint, -9*z - 3/Q*zpoint

def alt(z0,zpoint0):
    return odeint(f,[z0,zpoint0],t)[:,0]
    # Retourne la liste des angles (1ère colonne)

    
# Tracé
z0=5
zpoint0=0

plt.figure()
plt.title('Oscillateur amorti')
for Q in [.5 , 1, 5] :
    plt.plot(t,alt(z0,zpoint0),label=f'Q = {Q}')
plt.xlabel('temps en s')
plt.ylabel('altitude en m')
plt.legend()
plt.grid()         
plt.show()

Quatrième exemple : Oscillateur harmonique amorti en régime sinusoïdal

Soit l'équation différentielle :

\begin{equation} \mathsf{ \frac{d^2z}{dt^2} + 10 \cdot \frac{dz}{dt} + 9 \cdot z = 8 \cdot cos(5 \cdot t) } \end{equation}

Ce tracé permet d'illustrer qu'au-delà d'un régime transitoire, le régime devient sinusoïdal.

La fonction $\mathsf{f}$ doit être définie comme suit : \begin{equation} \mathsf{ f : \begin{pmatrix} z \\ zpoint \end{pmatrix} \mapsto \begin{pmatrix} zpoint \

        - 9 \cdot z - 10 \cdot zpoint + 8 \cdot cos(5 \cdot t)
\end{pmatrix}

} \end{equation}

In [75]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint


# Durée d'intégration
t=np.linspace(0,12,5000)


# Fonctions utiles à la résolution
def f(X,t):
    z,zpoint = X
    return zpoint, - 9*z - 10*zpoint + 8*np.cos(5*t)

def sol(z0,zpoint0):
    return odeint(f,[z0,zpoint0],t)[:,0]

    
# Tracé
z0=-0
zpoint0=8

plt.figure()
plt.plot(t,sol(z0,zpoint0))
plt.xlabel('temps en s')
plt.ylabel('z')
plt.grid()         
plt.show()

Cinquième exemple : Equation différentielle non linéaire

Soit l'équation différentielle :

\begin{equation} \mathsf{ \frac{d^2z}{dt^2} = 9 \cdot sin(z) - 3 \cdot \left( \frac{dz}{dt} \right)^2 } \end{equation}

La fonction $\mathsf{f}$ doit être définie comme suit : \begin{equation} \mathsf{ f : \begin{pmatrix} z \\ zpoint \end{pmatrix} \mapsto \begin{pmatrix} zpoint \\ 9 \cdot sin(z) - 3 \cdot \left( \frac{dz}{dt} \right)^2 \end{pmatrix} } \end{equation}

In [86]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint


# Durée d'intégration
t=np.linspace(0,.5,1000)


# Fonctions utiles à la résolution
def f(X,t):
    z,zpoint = X
    return zpoint, 9*np.sin(z) - 3*zpoint**2

def sol(z0,zpoint0):
    return odeint(f,[z0,zpoint0],t)[:,0]

    
# Tracé
z0=10
zpoint0=20

plt.figure()
plt.plot(t,sol(z0,zpoint0))
plt.xlabel('temps en s')
plt.ylabel('z')
plt.grid()         
plt.show()