Méthode d'Euler

Résumé en cinq points

Le problème à résoudre

Équation différentielle d'ordre 1 et problème de Cauchy

On cherche une fonction y du temps t qui vérifie y'(t) = f(t, y(t)) sur un intervalle [t₀, t_f], avec une condition initiale y(t₀) = y₀. C'est un problème de Cauchy. La fonction f est donnée, elle traduit la loi d'évolution. Exemples.

• Décroissance radioactive ou décharge d'un condensateur, y' = − y / τ, solution exacte y₀ exp(− (t − t₀) / τ).
• Modèle logistique de Verhulst, N' = r N (1 − N / K), où r est le taux de croissance et K la capacité d'accueil du milieu. La solution exacte est N(t) = K / (1 + (K / N₀ − 1) exp(− r t)) avec t₀ = 0. Elle croît de N₀ vers K en forme de S quand N₀ < K / 2.

Discrétisation

On ne calcule pas y partout mais en des instants régulièrement espacés, t_k = t₀ + k h, où h est le pas de discrétisation. Avec n pas, h = (t_f − t₀) / n. On construit des valeurs approchées y_k de y(t_k).

La méthode d'Euler explicite

Principe

Pour h petit, la dérivée est proche du taux d'accroissement, y'(t_k) ≈ (y(t_k+1) − y(t_k)) / h. En remplaçant y'(t_k) par f(t_k, y(t_k)), on obtient la relation de récurrence de la méthode d'Euler explicite.

y_k+1 = y_k + h f(t_k, y_k).

Interprétation géométrique. À partir du point (t_k, y_k), on avance d'un pas h en suivant la tangente de pente f(t_k, y_k). On remplace ainsi la courbe par une ligne brisée. La méthode est dite explicite car y_k+1 s'obtient directement à partir de y_k, sans résoudre d'équation.

Programme

def euler(f, t0, tf, y0, n):
    # f : fonction de deux variables, f(t, y). Renvoie les listes des t_k et des y_k.
    h = (tf - t0) / n
    t = [t0]
    y = [y0]
    for k in range(n):
        y.append(y[k] + h * f(t[k], y[k]))
        t.append(t[k] + h)
    return t, y

La fonction f est passée en paramètre, ce qui permet de résoudre n'importe quelle équation avec le même programme. Pour le modèle logistique, on définit f(t, N) = r N (1 − N / K), qui ne dépend pas de t mais doit en recevoir un. Les listes renvoyées se tracent directement, et contiennent n + 1 points.

Exemple d'utilisation, avec r = 0,5 par jour, K = 1 000 individus et N₀ = 10 sur 20 jours.

r = 0.5
K = 1000

def f_logistique(t, N):
    return r * N * (1 - N / K)

t, N = euler(f_logistique, 0, 20, 10, 200)

La fonction f_logistique utilise r et K définis dans l'espace global du module, conformément à la règle de portée lexicale. Avec 200 pas, h = 0,1 jour, soit vingt fois moins que le temps caractéristique 1 / r = 2 jours. La courbe obtenue a l'allure en S attendue et approche K à la fin de l'intervalle.

Coût

Chaque pas demande une évaluation de f et quelques opérations. Le temps de calcul est donc proportionnel au nombre de pas n = (t_f − t₀) / h. Diviser le pas par 10 multiplie le temps de calcul par 10.

Impact du pas de discrétisation

Qualité des résultats, erreur d'ordre 1

L'écart entre la valeur approchée y_n et la valeur exacte y(t_f) est l'erreur globale. On comprend son comportement en deux temps. Sur un seul pas, remplacer la courbe par sa tangente commet une erreur de l'ordre de h², d'après le développement de Taylor y(t + h) = y(t) + h y'(t) + (h² / 2) y''(t) + …, dont Euler ne garde que les deux premiers termes. Sur l'intervalle entier, il y a n = (t_f − t₀) / h pas, et ces erreurs locales s'accumulent, ce qui donne une erreur globale de l'ordre de n × h², c'est-à-dire de l'ordre de h. Pour la méthode d'Euler et une fonction f régulière, elle est à peu près proportionnelle au pas h, on dit que la méthode est d'ordre 1. Diviser h par 2 divise l'erreur finale par environ 2. Pour gagner un chiffre significatif, il faut dix fois plus de pas, donc dix fois plus de temps.

Exemple chiffré. Pour y' = − y avec y(0) = 1 sur [0, 1], la valeur exacte est e⁻¹ ≈ 0,3679. Avec un seul pas, h = 1, Euler donne y₁ = 1 − 1 = 0. Avec 10 pas, h = 0,1, on obtient (1 − 0,1)¹⁰ = 0,9¹⁰ ≈ 0,3487, erreur d'environ 0,019. Avec 100 pas, 0,99¹⁰⁰ ≈ 0,3660, erreur d'environ 0,0018. L'erreur est bien divisée par 10 quand le pas l'est.

On met cet effet en évidence en traçant sur un même graphique la solution exacte et plusieurs solutions approchées de pas différents, ou en traçant l'erreur finale en fonction de h en échelle logarithmique, où la pente vaut 1.

Instabilité quand le pas est trop grand

Pour y' = − a y avec a > 0, la récurrence s'écrit y_k+1 = (1 − a h) y_k. La solution exacte décroît vers 0 sans osciller. La suite approchée décroît sans osciller si a h < 1, oscille en changeant de signe si 1 < a h < 2, et diverge si a h > 2. Un pas trop grand devant le temps caractéristique 1 / a donne donc une courbe qui n'a plus rien à voir avec la réalité. Règle pratique, choisir un pas nettement inférieur au plus petit temps caractéristique du problème, par exemple au moins dix fois plus petit.

Compromis entre précision et temps de calcul

Réduire le pas améliore la précision tant que l'erreur de méthode domine. Pour des pas extrêmement petits, les erreurs d'arrondi des flottants, accumulées sur un très grand nombre de pas, finissent par compter, et le temps de calcul devient prohibitif. Le choix du pas est donc un compromis, qu'on règle en diminuant h jusqu'à ce que la courbe ne change plus visiblement.

Comparaisons avec d'autres résolutions

Comparaison avec une résolution exacte

Quand la solution exacte est connue, décroissance exponentielle, modèle logistique, on la trace à côté de la solution approchée. C'est le meilleur moyen de valider un programme et de mesurer l'erreur. Sur le modèle logistique, la solution d'Euler approche bien la sigmoïde pour un pas petit devant 1 / r, et prend de l'avance ou du retard visible pour un pas de l'ordre de 1 / r.

Méthode de Heun

La méthode de Heun corrige Euler en utilisant la pente moyenne sur le pas. On calcule d'abord une prédiction par Euler, p = y_k + h f(t_k, y_k), puis on corrige.

y_k+1 = y_k + (h / 2) (f(t_k, y_k) + f(t_k+1, p)).

Elle demande deux évaluations de f par pas au lieu d'une, mais elle est d'ordre 2. Son erreur est à peu près proportionnelle à h², diviser le pas par 10 divise l'erreur par environ 100. À précision égale, elle est donc bien plus rapide qu'Euler. Sur l'exemple y' = − y, un seul pas h = 1 donne déjà 0,5 au lieu de 0, et dix pas donnent une erreur de l'ordre de 7 × 10⁻⁴, contre 0,019 avec Euler.

Fonction de résolution fournie par un module

Les bibliothèques scientifiques fournissent des fonctions de résolution plus précises, qui adaptent elles-mêmes le pas. La fonction odeint du module scipy.integrate en est un exemple. Attention, elle attend une fonction dont les arguments sont dans l'ordre (y, t), et non (t, y). Le programme précise qu'aucune connaissance sur un module n'est exigible, sa documentation doit être fournie. On l'utilise comme référence pour juger la qualité d'une méthode programmée à la main. Concrètement, on calcule sur les mêmes instants la solution d'Euler, celle de Heun et celle du module, on trace les trois courbes, puis on trace l'écart de chaque méthode à la référence. On constate que, pour un même pas, l'écart de Heun est bien plus faible que celui d'Euler, et qu'il diminue beaucoup plus vite quand on réduit le pas.

Compléments, systèmes et équations d'ordre 2

Systèmes différentiels d'ordre 1

Ce complément fait partie des thèmes libres proposés par le programme, il n'est pas exigible. Pour un système de plusieurs équations couplées, on applique la même formule à chaque inconnue, en calculant toutes les pentes avec les valeurs de l'instant t_k avant de mettre à jour.

Modèle proies-prédateurs de Lotka-Volterra, x' = a x − b x y pour les proies, y' = − c y + d x y pour les prédateurs, avec a, b, c, d positifs. Un pas d'Euler s'écrit

def pas_lotka(x, y, h, a, b, c, d):
    nx = x + h * (a * x - b * x * y)
    ny = y + h * (-c * y + d * x * y)
    return nx, ny

Les deux nouvelles valeurs sont calculées à partir des anciennes. Mettre x à jour avant de calculer y avec la nouvelle valeur change la méthode et fausse la comparaison. Les solutions exactes sont périodiques, et la méthode d'Euler explicite les transforme en spirales qui s'écartent lentement, un défaut qui se réduit avec le pas.

Modèle épidémique SIR, S' = − β S I, I' = β S I − γ I, R' = γ I, où S, I et R sont les proportions de sensibles, d'infectés et de guéris. La somme S + I + R reste constante, et c'est un bon test de la programmation, puisque la méthode d'Euler la conserve exactement, aux arrondis près.

Équations d'ordre 2

Une équation d'ordre 2 se ramène à un système de deux équations d'ordre 1 en introduisant la dérivée comme nouvelle inconnue. Pour l'oscillateur amorti x'' + (ω₀ / Q) x' + ω₀² x = 0, on pose v = x', et on obtient x' = v, v' = − (ω₀ / Q) v − ω₀² x. On applique alors Euler au couple (x, v).

Pour un oscillateur non amorti, x'' + ω₀² x = 0, un pas d'Euler s'écrit x_k+1 = x_k + h v_k et v_k+1 = v_k − h ω₀² x_k. Un calcul direct montre que la quantité ω₀² x² + v², proportionnelle à l'énergie mécanique, est multipliée exactement par 1 + ω₀² h² à chaque pas. Elle croît donc géométriquement. La méthode d'Euler explicite fait croître l'énergie à chaque pas, la trajectoire dans le plan (x, v) spirale vers l'extérieur au lieu de rester sur une ellipse. C'est un défaut connu de la méthode, qu'on met en évidence en traçant l'énergie au cours du temps.

Système différentiel issu d'une cinétique chimique

Pour deux réactions successives A → B → C d'ordre 1, de constantes k₁ et k₂, on a [A]' = − k₁ [A], [B]' = k₁ [A] − k₂ [B], [C]' = k₂ [B]. Euler donne l'évolution des trois concentrations, et le maximum de [B] se lit sur la courbe. Le pas doit rester petit devant 1 / k₁ et 1 / k₂.

Méthodes

• Écrire la récurrence y_k+1 = y_k + h f(t_k, y_k) avant le programme, et identifier f, t₀, t_f, y₀ et n.
• Vérifier le nombre de points, n + 1, et le dernier instant, t_f.
• Choisir un pas petit devant le plus petit temps caractéristique, puis le diviser par 2 ou par 10 pour vérifier que la courbe ne change plus.
• Valider le programme sur un cas à solution exacte connue avant de l'utiliser sur un modèle sans solution explicite.
• Pour un système, calculer toutes les pentes avec les anciennes valeurs avant toute mise à jour.
• Pour une équation d'ordre 2, poser v = x' et se ramener à un système d'ordre 1.
• Commenter une courbe obtenue en reliant son allure au pas choisi, retard, oscillations parasites, divergence.

Ce qu'il faut savoir refaire en khôlle

1. Établir la relation de récurrence de la méthode d'Euler explicite à partir du taux d'accroissement, et l'interpréter géométriquement.
2. Écrire la fonction euler prenant f en paramètre et l'appliquer au modèle logistique.
3. Comparer la solution d'Euler à la solution exacte de y' = − y, et montrer que l'erreur est proportionnelle au pas.
4. Montrer que pour y' = − a y, un pas tel que a h > 2 fait diverger la méthode.
5. Présenter la méthode de Heun et expliquer pourquoi elle est d'ordre 2.
6. Résoudre par Euler le système de Lotka-Volterra ou le modèle SIR, et ramener une équation d'ordre 2 à un système d'ordre 1.

Les résultats à connaître

En informatique, les « résultats » sont la correction des méthodes, leur coût et leur comportement quand le pas varie.

• Récurrence d'Euler explicite y_k+1 = y_k + h f(t_k, y_k), avec t_k = t₀ + k h et h = (t_f − t₀)/n : attendu, à savoir établir.
• n pas produisent n + 1 valeurs, de y₀ à y_n, et le dernier instant vaut t_f aux arrondis près : attendu.
• Coût : n évaluations de f pour Euler, 2n pour Heun. Le temps de calcul est proportionnel à 1/h : attendu, à savoir justifier.
• Si f ne dépend que de t, Euler coïncide avec la méthode des rectangles à gauche appliquée à ∫ f, et Heun avec la méthode des trapèzes : à savoir établir.
• Pour y' = a y, Euler donne y_k = (1 + a h)^k y₀, une suite géométrique : à savoir établir.
• Pour y' = − a y avec a > 0 : la suite décroît vers 0 sans changer de signe si a h < 1. Elle tend vers 0 en alternant si 1 < a h < 2. Elle ne tend pas vers 0 si a h ≥ 2 : à savoir démontrer.
• Erreur globale d'Euler à peu près proportionnelle à h, méthode d'ordre 1 : admis en général, à illustrer numériquement. Le cas y' = a y se démontre (exercice 9).
• Heun est d'ordre 2 : admis, à illustrer numériquement.
• Une fonction de module comme odeint sert de référence. Aucune connaissance sur elle n'est exigible, sa documentation est fournie.
• Système d'ordre 1 : toutes les pentes se calculent avec les valeurs de l'instant t_k : thème libre, non exigible.
• Équation d'ordre 2 ramenée à un système en posant v = x' : thème libre, non exigible.
• Oscillateur non amorti : Euler multiplie ω₀² x² + v² par 1 + ω₀² h² à chaque pas : thème libre, démontrable par un calcul direct.
• Programmation vectorisée : peut être proposée, n'est pas un attendu du programme.
• Hors programme signalé : méthode d'Euler implicite, méthodes de Runge-Kutta d'ordre supérieur, pas adaptatif.

Clés de démonstration

Récurrence d'Euler. On intègre l'équation entre t_k et t_k+1, ce qui donne y(t_k+1) = y(t_k) + ∫ f(s, y(s)) ds. On approche l'intégrale par la valeur à gauche multipliée par h. L'astuce est ce point de vue intégral : il donne aussi Heun, qui remplace les rectangles par les trapèzes.

Stabilité de y' = − a y. La récurrence est géométrique de raison q = 1 − a h, donc y_k = q^k y₀. Elle tend vers 0 si et seulement si |q| < 1, c'est-à-dire 0 < a h < 2, et elle garde un signe constant si q ≥ 0. L'astuce est de tout ramener à la raison d'une suite géométrique.

Ordre 1 pour y' = a y. On écrit y_n = y₀ exp(n ln(1 + a T/n)) avec T = t_f − t₀, puis on développe le logarithme à l'ordre 2. L'astuce est le développement limité ln(1 + x) = x − x²/2 + o(x²), qui fait apparaître un terme en 1/n, donc en h.

Énergie de l'oscillateur. On développe ω₀² (x + h v)² + (v − h ω₀² x)². Les termes croisés 2 h ω₀² x v s'annulent, et il reste (1 + ω₀² h²)(ω₀² x² + v²). L'astuce est de repérer cette annulation avant de développer entièrement.

Conservation dans SIR. La somme des trois pentes est nulle en tout point. Un pas d'Euler ajoute h fois cette somme à S + I + R, donc ne la change pas. L'astuce vaut pour tout invariant linéaire : Euler le conserve exactement, aux arrondis près.

Exemples et contre-exemples à connaître

1. Pour y' = y, y(0) = 1 sur [0, 1], Euler donne (1 + 1/n)ⁿ : 2 pour n = 1, 2,5937 pour n = 10, 2,7048 pour n = 100, contre e ≈ 2,7183. La méthode sous-estime, car la solution est convexe et la tangente passe sous la courbe.
2. Pour y' = 2t, y(0) = 0, sur [0, 1], Euler donne exactement 1 − h au lieu de 1. C'est le cas où l'erreur est exactement proportionnelle au pas.
3. Heun est exacte pour y' = 2t, car la méthode des trapèzes est exacte pour une fonction affine.
4. Pour y' = − a y, le choix a h = 1 donne y₁ = 0, puis y_k = 0 pour tout k ≥ 1. Le choix a h = 2 donne y₀, − y₀, y₀… sans amortissement.
5. Pour le modèle logistique avec 1 < r h < 2, la suite d'Euler dépasse K puis oscille autour. Or la vraie solution ne franchit jamais K quand N₀ < K. Avec r h juste au-dessus de 2, elle alterne entre deux valeurs.
6. Pour y' = y², y(0) = 1, la solution exacte 1/(1 − t) explose en t = 1. Euler sur [0, 2] ne le signale pas : valeur finie absurde, de l'ordre de 10¹⁰⁴, avec h = 0,1, ou inf par dépassement de capacité avec h = 0,01, sans message si f calcule y * y. Il faut connaître le modèle avant de faire confiance au calcul.
7. Pour l'oscillateur non amorti avec ω₀ = 1 et h = 0,1, l'énergie double en 70 pas environ, soit un peu plus d'une période. La trajectoire en spirale n'est pas physique.
8. Dans SIR, S + I + R reste constant avec une mise à jour simultanée. Une dérive de cette somme signale une erreur de programmation.
9. odeint attend une fonction de (y, t). Lui passer la fonction f(t, y) écrite pour Euler échange les rôles de t et de y. En général, aucun message n'apparaît et le résultat est faux. Avec f(t, y) = − y, odeint résout en fait y' = − t.

Techniques du chapitre

Technique 1, mettre l'équation sous forme résolue. Mini-exercice : circuit RC, R C u' + u = E. On isole u' : f(t, u) = (E − u)/(R C). Le temps caractéristique est τ = R C, le pas se choisit devant τ.

Technique 2, faire deux pas à la main. Mini-exercice : y' = t − y, y(0) = 1, h = 0,5. y₁ = 1 + 0,5 × (0 − 1) = 0,5, puis y₂ = 0,5 + 0,5 × (0,5 − 0,5) = 0,5. La solution exacte t − 1 + 2 e^−t vaut 0,7358 en t = 1. Ce calcul sert de test pour le programme.

Technique 3, choisir le pas à partir des temps caractéristiques. Mini-exercice : réactions A → B → C avec k₁ = 0,2 s⁻¹ et k₂ = 5 s⁻¹, sur 30 s. Le plus petit temps est 1/k₂ = 0,2 s. On prend h ≤ 0,02 s, donc n ≥ 1 500.

Technique 4, mesurer l'ordre par le rapport des erreurs. Mini-exercice : pour y' = − y sur [0, 1], l'erreur finale vaut 0,0192 avec n = 10 et 0,0094 avec n = 20. Le rapport est proche de 2, donc l'ordre est 1. Pour Heun, on trouve un rapport proche de 4, donc l'ordre est 2.

Technique 5, compter les évaluations de f. Mini-exercice : Euler à 1 000 pas contre Heun à 100 pas. Le coût est de 1 000 évaluations contre 200. Pour y' = − y, les erreurs sont d'environ 1,8 × 10⁻⁴ et 6 × 10⁻⁶. Heun gagne sur les deux plans.

Technique 6, calculer les pentes avant de mettre à jour. Mini-exercice : pour x' = − y et y' = x, on écrit dx = -y et dy = x, puis x = x + h * dx et y = y + h * dy. Les deux pentes utilisent bien les valeurs de t_k.

Technique 7, tester sur une équation qui dépend de t. Mini-exercice : un programme qui évalue f au mauvais instant donne le bon résultat sur y' = − y, car f n'y dépend pas de t. Sur y' = 2t avec n = 2, on doit trouver 0,5. Une évaluation en t_k+1 donne 1,5.

Exercices classiques corrigés

Choisis ton niveau. Difficulté indicative ; les exercices qui se suivent restent dans l'ordre.

def euler_final(f, t0, tf, y0, n):
    # renvoie y_n, approximation de y(tf)
    h = (tf - t0) / n
    y = y0
    for k in range(n):
        t = t0 + k * h
        y = y + h * f(t, y)
    return y

from math import exp

def ecart_logistique(r, K, N0, T, n):
    h = T / n
    N = N0
    emax = 0
    for k in range(1, n + 1):
        N = N + h * r * N * (1 - N / K)
        e = N - K / (1 + (K / N0 - 1) * exp(-r * k * h))
        if e < 0:
            e = -e
        if e > emax:
            emax = e
    return emax

def heun(f, t0, tf, y0, n):
    h = (tf - t0) / n
    t = [t0]
    y = [y0]
    for k in range(n):
        p1 = f(t[k], y[k])
        p2 = f(t[k] + h, y[k] + h * p1)
        y.append(y[k] + h / 2 * (p1 + p2))
        t.append(t[k] + h)
    return t, y

def sir_pic(beta, gamma, S0, I0, T, n):
    h = T / n
    S = S0
    I = I0
    imax = I0
    tmax = 0
    for k in range(1, n + 1):
        dS = -beta * S * I
        dI = beta * S * I - gamma * I
        S = S + h * dS
        I = I + h * dI
        if I > imax:
            imax = I
            tmax = k * h
    return imax, tmax

Erreurs fréquentes