In [4]:
# Pensez à évaluer cette cellule ! (Shift-Entrée)

%xmode Plain

import numpy as np
import math

import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.style.use('fivethirtyeight')  # https://matplotlib.org/3.1.1/gallery/style_sheets/style_sheets_reference.html
%matplotlib inline

from scipy.integrate import odeint
Exception reporting mode: Plain

Conventions

Attention, dans l'énoncé original, certaines conventions sont différentes...

Pour simplifier les choses et coller aux conventions de la fonction odeint :

  • On considère des équations différentielles de la forme $y'(t) = F(y(t), t)$. Arguments : y d'abord, puis t
  • On définira euler(F, y0, temps) (et idem pour les suivantes) où les temps seront donnés dans une liste (ou array) créé à l'avance.

Méthode d'Euler

Schéma d'Euler explicite : $$y_{k+1} = y_k + (t_{k+1}-t_k) \times F(y_k, t_k)$$

In [15]:
def euler(F, y0, temps):
    y = np.array(y0)              # On travaille avec des np.array : Vectorialisation
    res = [y]
    for k in range(len(temps)-1):
        t = temps[k]
        h = temps[k+1]-temps[k]
        pente = np.array(F(y, t)) # Approximation du schéma d'Euler
        y = y + h * pente
        res.append(y)
    return np.array(res)          # Pour avoir un tableau à la fin
In [19]:
def F(y, t):
    """ EquaDiff : y' = y """
    return y

x = np.linspace(0, 1, 10)

y_solution = np.exp(x) # Uniquement quand on connaît la solution analytique !
plt.plot(x, y_solution, label="Vraie solution")

# Les lignes importantes !
y_euler = euler(F, 1, x)
plt.plot(x, y_euler, '+', label="Méthode d'Euler")

plt.legend()
plt.show()

Exercice 1 : Tester le code ci-dessus et poser des questions si quelque chose n'est pas clair...

  • Méthode d'Euler
  • Vectorialisation
  • Graphiques avec matplotlib

Variantes : autres méthode explicites

$$y_{k+1} = y_k + (t_{k+1}-t_k) \times pente$$

où la $pente$ est calculable (parfois en plusieurs étapes) à partir de $y_k$ et $t_k$

Méthode de Heun

$$y_{k+1} = y_k + (t_{k+1}-t_k) \times \frac{1}{2}(F(y_k, t_k) + F(y_k + hF(y_k, t_k), t_{k+1}))$$
  • Faites un copier-coller du code ci-dessus, il "suffit" de modifier le calcul de la pente
  • Prenez votre temps. On pourra par exemple noter $p_1$ la valeur $F(y_k, t_k)$ et utiliser le fait que $t_{k+1} = t_k + h$
  • Essayez de faire simple, et procédez étape par étape
  • SyntaxError: ... : avez-vous bien fermé toutes les parenthèses à la ligne du dessus ?
  • à chaque fois que vous appelez F, écrivez np.array(F(..., ...))

Exercice 2 : Écrire la méthode de Heun en Python

In [29]:
def heun(F, y0, temps):
    ...

Trop compliqué ? à chaque étape, on calcule la pente de cette manière :

  • $p_1 = F(y_k, t_k)$
  • $p_2 = F(y_k + h \times p_1, t_k+h)$
  • $pente = (p_1+p_2)/2$
  • $y_{k+1} = y_k + h \times pente$

Méthode de Runge-Kutta : RK4

Exactement le même concept :

  • $p_1 = F(y, t)$
  • $p_2 = F(y+h/2 \times p_1, t+h/2)$
  • $p_3 = F(y+h/2 \times p_2, t+h/2)$
  • $p_4 = F(y+h \times p_3, t+h)$
  • $pente = \frac{1}{6}(p_1 + 2p_2 + 2p_3 + p_4)$

Exercice 3 : Écrire la méthode Runge-Kutta (RK4) en Python

In [18]:
def rk4(F, y0, temps):
    ...

odeint

Rien à écrire : la fonction odeint (Ordinary Differential Equation) du module scipy.integrate fonctionne exactement sur le même modèle. Donc vous pouvez utiliser odeint(F, y0, temps)

Graphiques

Exercice 4 : en vous inspirant des graphiques ci-dessus, toujours avec l'équation différentielle $y'=y$

In [22]:
# à vous

def F(y, t):
    return ...

N = 30 # nombre de points
x = np.linspace(0, 1, N)

# Tracer les courbes correspondant aux méthodes d'Euler, Heun, RK4, odeint.

plt.legend()
plt.show()
No handles with labels found to put in legend.

Ordre supérieur

On considère l'équation différentielle $y'' = -sin(y)$

Papier-Stylo (pas très difficile) : écrire la vectorialisation et comprendre pourquoi le code ci-dessous et cohérent

In [26]:
def F_pendule(z, t):
    """ F((y, y'), t) -> (y', -sin(y))"""
    y, yp = z
    return (yp, -np.sin(y))

x = np.linspace(0, 10, 100)
y_euler = euler(F_pendule, (0, 1.5), x) # On note ici que y0 est un couple : (y(t0), y'(t0))

plt.plot(x, y_euler, label="Méthode d'Euler")
plt.legend()
Out[26]:
<matplotlib.legend.Legend at 0x7f57c33fa8d0>

On remarque que Python a tracé deux courbes (avec la même légende...), alors qu'a priori on a demandé un seul tracé... C'est parce que euler renvoie ici un tableau à $n$ lignes et 2 colonnes

In [31]:
# On extrait chacune des colonnes.
# Syntaxe possible uniquement parce qu'on a un tableau np.array. Comment ferait-on avec des listes Python ?

y = y_euler[:, 0]
yp = y_euler[:, 1]

plt.plot(x, y, label="Euler : solution")
plt.plot(x, yp, label="Euler : dérivée")
plt.legend()
plt.show()

plt.plot(y, yp, label="Portrait de phase")
plt.legend()
plt.show()

Comparer Euler / Heun / RK4 / odeint

Exercice 5 : en utilisant vos fonctions ci-dessus, tracez les courbes correspondant à la résolution avec Euler, Heun, RK4 et odeint. Vous fonctions sont-elles bien adaptées à l'ordre supérieur (utilisation de np.array) ?

Testez, essayez, comparez...

In [ ]:
# à vous

Autres applications

Exercice :

  • Reprendre l'énoncé original : partie 4 et suivantes
  • Attention aux conventions qui sont parfois différentes !
In [ ]:
# à vous

Aller (beaucoup) plus loin

  • Il existe de nombreuses variantes de type Runge-Kutta. On peut les classifier facilement avec les tableaux de Butcher. En général, les coefficients utilisés sont obtenus avec les formules de Taylor.

  • Il existe des méthodes dites implicites (BDF: Backward Differentiation Formula, ou méthodes de Gear). Une des plus simples/connues est la méthode d'Euler implicite. Il est nécessaire de résoudre une équation à chaque étape ; par exemple avec la méthode de Newton-Raphson

  • Il existe d'autres méthodes plus efficaces (conservation de l'énergie, réversibilité) pour les applications physiques : méthode de Verlet, Leap Frog...

  • Dans odeint, on utilise soit une méthode BDF, soit une méthode de Adams selon que l'équation différentielle est "raide" (la dérivée de la solution est très grande) ou pas.

  • Il existe de nombreuses bibliothèques de résolution numériques des équations différentielles ordinaires. À ma connnaissance, la plus complète est écrite dans la langage Julia : OrdinaryDiffEq.jl

In [ ]: