Méthode d'Euler et thermalisation en Python

Objectif de la fiche : utiliser un langage de programmation pour résoudre une équation différentielle grâce à la méthode d'Euler, à travers l'exemple de la thermalisation d'un système en contact avec un thermostat.

Cette fiche est autonome : tu n'as besoin d'aucun autre document. Chaque étape contient un script incomplet à compléter, puis à copier-coller dans Thonny (ou dans la console Basthon tout en bas de la page si tu n'as pas Thonny sous la main).

＋

Contexte : simuler sans résoudre à la main

À savoir

Un système à la température T, en contact avec un thermostat à la température T_S, voit sa température évoluer selon l'équation différentielle dT/dt = -a/c·(T - T_S), où a et c sont des constantes liées au système. Cette loi traduit un échange thermique d'autant plus rapide que l'écart de température est grand.
Cette équation possède une solution mathématique exacte, mais on va ici la retrouver numériquement, sans la résoudre à la main : c'est l'intérêt de la méthode d'Euler.
Le principe : on écrit l'équation différentielle sous la forme dT/dt = f(T, t), puis on avance pas à pas depuis la condition initiale, en approchant la courbe solution par une suite de petits segments de droite. Plus le pas est petit, plus l'approximation est fine — mais plus le calcul est long.

La fonction f(T,t) et les paramètres

À savoir

On simule un système initialement à T0 = 300 K, plongé dans un thermostat à TS = 400 K, avec des constantes a = 1 et c = 1.
La fonction f(T, t) renvoie la valeur de dT/dt à un instant donné : c'est elle qui encode l'équation différentielle.
Le script ci-dessous est déjà complet : copie-le tel quel dans Thonny (ou dans Basthon) pour commencer.

etape0_equation.py

import numpy as np
import matplotlib.pyplot as plt

# bornes de l'intervalle de temps et nombre de pas pour la résolution
t0, tf = 0., 10.
N = 1000
t = np.linspace(t0, tf, N+1)  # liste des dates pour le calcul des solutions approchées

# condition initiale et constantes du système
T0 = 300  # température initiale, en K
a = 1
c = 1
TS = 400  # température du thermostat, en K

def f(T, t):
    return -a*(T-TS)/c

Indice

Le paramètre t de f(T, t) n'est pas utilisé dans le return ici (l'équation ne dépend pas explicitement du temps), mais il faut tout de même le prévoir dans la signature de la fonction : la méthode d'Euler, plus générale, doit pouvoir fonctionner aussi avec des équations où t intervient.

Écrire la fonction euler

À savoir

La méthode d'Euler calcule T(i+1) à partir de T(i) grâce à la relation T(i+1) = T(i) + h·f(T(i), t(i)), où h est le pas (la durée entre deux instants successifs).
Comme les instants de t sont régulièrement espacés (grâce à np.linspace), h est constant : il suffit de le calculer une seule fois, comme la différence entre les deux premières valeurs de t.
La fonction doit fonctionner pour n'importe quelle équation différentielle : c'est pour cela qu'on lui passe la fonction g (qui jouera le rôle de f), la condition initiale x0, et le tableau des instants t, sans jamais écrire f, T0 ou TS à l'intérieur.

Consigne

Complète la fonction euler ci-dessous, qui renvoie un tableau numpy x contenant la solution approchée à chaque instant de t.
N'oublie pas d'initialiser x avec np.zeros(N) avant la boucle : un tableau numpy a une taille fixée dès sa création, contrairement à une liste, on ne peut pas lui ajouter des éléments avec .append(...) au fur et à mesure.

etape1_fonction_euler.py

def euler(g, x0, t):  # méthode d'Euler
    h = t[_] - t[_]  # TODO : pas constant, distance entre les deux premières valeurs de t
    N = len(_)  # on recalcule N en local pour le cas général
    x = np.zeros(_)  # TODO : initialisation obligatoire du tableau numpy
    x[_] = __  # TODO : condition initiale
    for i in range(N-_) :
        x[i+_] = x[_] + h*g(x[_], t[_])  # TODO : définit x[i+1] par Euler
    return x

Indice

range(N-1) permet à i de parcourir 0, 1, ..., N-2 : à chaque tour, on calcule x[i+1], ce qui remplit progressivement le tableau jusqu'au dernier indice N-1, sans jamais chercher à écrire au-delà.

Résultat attendu

Cette fonction ne produit aucun affichage à elle seule (elle se termine par return x) : le résultat se vérifie à l'étape suivante, en l'appliquant au problème de la thermalisation.

Résoudre et tracer avec N = 1000 pas

Consigne

Appelle ta fonction euler avec la fonction f, la condition initiale T0 et le tableau t définis à l'étape 0, et stocke le résultat dans solE.
Trace solE en fonction de t, en bleu et en trait plein ('b-'), avec une légende qui indique la valeur de N utilisée.

etape2_resolution_n1000.py

# TODO : résous l'équation avec ta fonction euler
solE = euler(_, __, _)

plt.figure()
# TODO : trace solE en fonction de t, en bleu trait plein, avec une légende indiquant N
plt.plot(_, ____, '__', label='________________'+str(N))
plt.legend(loc='best')
plt.xlabel('t') ; plt.ylabel('$T$')
plt.grid()
plt.show()

Résultat attendu

Une courbe qui part de 300 K et remonte rapidement vers 400 K (la température du thermostat), en s'aplatissant progressivement. Pour vérifier ton calcul : solE[1] doit valoir 301 (le premier petit pas, de taille h = 0,01).

Comparer avec seulement 15 pas

À savoir

Ta fonction euler est générique : elle fonctionne avec n'importe quel tableau d'instants t, quelle que soit sa taille. Il suffit de lui fournir un nouveau tableau, construit avec un nombre de pas beaucoup plus petit.

Consigne

Crée un second tableau d'instants t2, sur le même intervalle [0, 10], mais avec seulement M = 15 pas.
Appelle euler une seconde fois avec ce nouveau tableau, et stocke le résultat dans solE2.
Superpose les deux courbes sur le même graphe : solE en bleu trait plein, solE2 en rouge avec des ronds ('ro'), chacune avec une légende indiquant son nombre de pas.

etape3_comparaison_pas.py

# TODO : construis un tableau de 15 pas sur le même intervalle de temps
M = __
t2 = np.linspace(__, __, ____)

# TODO : résous l'équation avec ce nouveau tableau d'instants
solE2 = euler(_, __, __)

plt.figure()
plt.plot(t, solE, 'b-', label='Euler avec N = '+str(N))
# TODO : superpose solE2 en fonction de t2, en rouge avec des ronds
plt.plot(__, _____, '__', label='________________'+str(M))
plt.legend(loc='best')
plt.xlabel('t') ; plt.ylabel('$T$')
plt.grid()
plt.show()

Indice

Il n'est pas nécessaire de réécrire toute la fonction euler : elle a été conçue pour accepter n'importe quel tableau t en argument, c'est tout l'intérêt de l'avoir rendue générale à l'étape 1.

Résultat attendu

Les points rouges (15 pas) suivent grossièrement la courbe bleue (1000 pas) mais s'en écartent nettement au début : avec un pas h beaucoup plus grand (h ≈ 0,67 s au lieu de 0,01 s), chaque segment de droite est une moins bonne approximation de la vraie courbe. Pour vérifier : solE2[1] doit valoir environ 367.

＋

Conclusion

À retenir

La méthode d'Euler permet de résoudre numériquement une équation différentielle, sans avoir besoin d'en connaître la solution mathématique exacte, en avançant pas à pas depuis la condition initiale.
Le choix du pas (ou, de façon équivalente, du nombre de pas N) est déterminant : un pas trop grand donne une solution approchée grossière, qui peut s'écarter significativement de la réalité, en particulier là où la fonction varie vite.
Écrire une fonction générique comme euler(g, x0, t), qui ne mentionne jamais explicitement le problème traité, permet de la réutiliser telle quelle pour n'importe quelle équation différentielle du premier ordre, et pour n'importe quel découpage temporel.

Pense-bête — méthode d'Euler

dT/dt=f(T,t)	forme générale d'une équation différentielle du 1er ordre
h = t[1]-t[0]	pas constant, déduit du tableau des instants
T(i+1)=T(i)+h·f(T(i),t(i))	relation d'Euler : avancer d'un pas h

np.zeros(N)	crée un tableau numpy de N zéros, à remplir ensuite
np.linspace(a,b,N+1)	engendre N+1 instants régulièrement espacés entre a et b
'b-' / 'ro'	ligne bleue continue / ronds rouges (couleur + marqueur)

Pas de Thonny sous la main ? Utilise Basthon

Ouvrir dans un nouvel onglet ↗

Basthon est une console Python qui tourne directement dans le navigateur, sans rien installer. Colle ici (avec Ctrl+V ou Cmd+V) le code copié avec les boutons « Copier le code » ci-dessus, puis appuie sur Entrée pour l'exécuter.