# 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)$$

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

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

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

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$

# à 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

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()

<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

# 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...

# à vous

Autres applications¶

Exercice :

Reprendre l'énoncé original : partie 4 et suivantes
Attention aux conventions qui sont parfois différentes !

# à 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