Bien lire et comprendre les trois exemples qui montrent comment approcher avec la fonction odeint la solution des EDO du premier ordre, des systèmes d’EDO du premier ordre ou d’une EDO d’ordre supérieur à 1 après transformation en un système d’EDO du premier ordre.
Apprendre à tracer le champ de vecteurs associé à une EDO du premier ordre ou le portrait de phase associé à un système d’EDO du premier ordre avec la fonction quiver.
Puis résoudre les deux premiers exercices proposés.
Exemples :
EDO du premier ordre + champ de vecteurs
Système d’EDO du premier ordre + portrait de phase
EDO d’ordre 2 \(\leadsto\) système d’EDO du premier ordre + portrait de phase
Exercices :
EDO du premier ordre
Système d’EDO du premier ordre
Bonus : étude de la stabilité d’un système d’EDO du premier ordre
Code
import sys #only needed to determine Python version numberprint(f'Python version: {</span>sys<span class="sc">.</span>version<span class="sc">}')import numpy as np # only needed to determine Numpy version numberprint(f"Numpy version: {</span>np<span class="sc">.</span>__version__<span class="sc">}")import scipy # only needed to determine Numpy version numberprint(f"Scipy version: {</span>scipy<span class="sc">.</span>__version__<span class="sc">}")from datetime import datetimeprint('Last Updated On: ', datetime.now())
Python version: 3.12.3 (main, Jan 17 2025, 18:03:48) [GCC 13.3.0]
Numpy version: 1.26.4
Scipy version: 1.11.4
Last Updated On: 2025-02-26 09:45:00.043319
Ce cours est une introduction aux méthodes numériques d’approximation d’Équations Différentielles Ordinaires (EDO) en utilisant Python comme langage de programmation. Au cours de ce module, nous explorerons, implémenterons et testerons plusieurs méthodes numériques.
Le module SciPy (et d’autres modules plus spécifiques) offre déjà des méthodes numériques pour résoudre des EDO, y compris celles que nous allons étudier dans ce module. Bien que ces méthodes soient déjà implémentées dans des fonctions prêtes à l’emploi, il est important de comprendre les principes de base de ces méthodes pour les utiliser correctement.
Pour illustrer l’utilisation pratique de ce module, nous allons voir sur trois exemples comment utiliser la fonction odeint du module SciPy pour approcher la solution d’abord d’une EDO ensuite d’un système d’EDO (ce qui inclut les équations différentielles d’ordre 2 ou plus).
Voici la signature de la fonction odeint pour référence :
Cette fonction résout un problème de Cauchy d’ordre 1, c’est-à-dire une ou plusieurs EDO d’ordre 1, définies par la fonction func et une condition initiale y0. La fonction func doit être définie par l’utilisateur et doit renvoyer la dérivée de y par rapport à t à un point donné t.
La fonction odeint renvoie un tableau de valeurs de y pour chaque valeur de t spécifiée dans le tableau t: dy/dt = func(y, t, ...) [or func(t, y, ...)]
Les arguments principaux sont :
func la fonction qui définit le problème de Cauchy,
y0 la condition initiale,
t le tableau des temps où on veut approcher la solution.
La méthode d’intégration utilisée est LSODA (méthode Adam/BDF adaptative).
Le principe d’utilisation de odeint (pour intégrer numériquement des équations différentielles) est le suivant: pour avoir une estimation numérique de la solution du problème de Cauchy
avec \(\mathbf{y}(t)=(y_1(t),y_2(t),\dots,y_n(t))\) le vecteur des fonctions recherchées dépendant de la variable \(t\) et \(\boldsymbol{\varphi}=(\varphi_1,\varphi_2,\dots,\varphi_n)\) une fonction de \(\mathbb{R}\times\mathbb{R}^n\to\mathbb{R}^n\), on donne comme argument la fonction \(\boldsymbol{\varphi}\) (qui doit avoir deux paramètres, même dans le cas autonome), la condition initiale \(\mathbf{y}_0\) et les points de l’interval de temps où approcher la solution (qui commence à \(t_0\)). Elle retourne un tableau Numpy (même si \(t\) était une liste).
Par défaut, l’ordre requis des deux premiers arguments de func est dans l’ordre inverse de celui classique. Pour utiliser une fonction avec la signature func(t, y, ...), c’est-à-dire d’abord \(t\) puis \(y\), l’argument tfirst doit être défini à True.
Remarque: la résolution de \(y''(t)=F(t,y(t),y'(t))\) passera par celle du système différentiel \(
\begin{cases}
y_1'(t)=y_2(t),\\
y_2'(t)=F(t,y_1(t),y_2(t))
\end{cases}
\) avec \(y(t)=y_1(t)\), \(y'(t)=y_1'(t)=y_2(t)\) et \(\varphi_1(t,y_1(t),y_2(t))=y_2(t)\), \(\varphi_2(t,y_1(t),y_2(t))=F(t,y_1(t),y_2(t))\).
Et la fonction solve_ivp ?
SciPy propose une interface plus récente et plus personnalisable pour résoudre un problème de Cauchy, il s’agit de solve_ivp. L’avantage principal est que solve_ivp offre plusieurs méthodes pour résoudre les équations différentielles alors qu’odeint est limité à une seule. Cependant, la prise en main de odeint est plus simple et suffisante pour les objectifs de ce cours.
fun la fonction qui définit le problème de Cauchy,
t_span un tuple qui contient l’intervalle de temps,
y0 la condition initiale,
method la méthode d’intégration,
t_eval le tableau des temps où on veut approcher la solution.
Parmi les méthodes d’intégration, on trouve RK45 (méthode explicite de Runge-Kutta d’ordre 5(4)), RK23 (méthode explicite de Runge-Kutta d’ordre 3(2)), DOP853 (méthode de Dormand-Prince d’ordre 8(5)), Radau (méthode implicite de Runge-Kutta d’ordre 5), BDF (méthode implicite multipas d’ordre variable) et LSODA (méthode Adam/BDF adaptative).
Les méthodes explicites de Runge-Kutta (RK23, RK45, DOP853) doivent être utilisées pour les problèmes non stiff et les méthodes implicites (‘Radau’, ‘BDF’) pour les problèmes stiff (on verra ce concept au CM-5). Parmi les méthodes de Runge-Kutta, DOP853 est recommandée pour la résolution avec une grande précision.
En cas de doute, essayez d’abord d’exécuter RK45. S’il effectue un nombre anormalement élevé d’itérations, diverge ou échoue, il est probable que le problème soit stiff et que l’on doit utiliser Radau ou BDF. LSODA peut également être un bon choix universel, mais il peut être un peu moins pratique à utiliser car il intègre un ancien code Fortran.
1.1 ✨ Champ de vecteurs et lignes de courant pour une EDO du premier ordre
Bien qu’il soit trés rare que l’on puisse résoudre explicitement une équation différentielle donnée, on peut souvent avoir une idée de l’allure des graphes des solutions en observant le champs de vecteurs associé. En effet le graphe d’une solution de l’équation \(y'(t)=\varphi(t,y(t))\) est par définition tangent à son vecteur vitesse \((1,y'(t))\) et donc au vecteur \((1, \varphi(t,y(t))\). La connaissance de la fonction \(\varphi\) en chaque point \((t,y)\) permet donc de représenter facilement ces vecteurs tangents même si l’on ne connait pas les solutions. Et si l’on en trace un grand nombre, uniformément répartis dans le plan \((t, y)\), on obtient une représentation du champ de vecteurs associé à l’équation différentielle qui permet souvent de deviner les graphes des solutions puisqu’il s’agit des courbes qui sont tangentes en tous leurs points aux vecteurs de ce champs de vecteurs.
Il est alors intéressant d’y superposer la solution numérique obtenue avec odeint.
La fonction quiver du module matplt.plotlib permet de tracer un champ de vecteurs: quiver(x_pos, y_pos, x_dir, y_dir, color)
La fonction streamplt.plot permet de tracer des lignes de courant représentant un champ vectoriel.
Utilisons-les pour obtenir les tracés associés à notre équation différentielle. La courbe rouge est la solution déterminée par odeint.
# Pour aller plus loin, décommenter un des deux appels suivants# =============================================================================# help(plt.quiver)# help(plt.streamplot)
1.2 ✨ Exemple de résolution approchée d’une équation différentielle d’ordre 1 : équation logistique simple
En guise d’exemple, considérons le problème de Cauchy \(y(0)=1\) et une équation logistique simple de la forme
NB: le nombre de points en lesquels les résultats sont évalués n’est pas (du moins directement) relié à la précision des calculs internes (ne pas imaginer que cela fixe le pas de la méthode, en particulier).
Superposons enfin le champ de vecteurs associé à cette équation différentielle à la solution approchée.
Code
# define a plt.grid g1 = np.linspace(0,5,21)g2 = np.linspace(0,8,21)T,Y = np.meshgrid(g1,g2) # create a plt.grid# compute direction at each pointDT, DY =1, phi(T,Y) # compute growth rate on the plt.grid# M = np.sqrt(DT**2+DY**2) # norm growth rate M = np.hypot(DT,DY)# Deux types de représentations : champ des pentes et lignes de courant# =============================================================================plt.figure(figsize=(20,7))# Champ des pentes# =============================================================================plt.subplot(1,2,1)# plt.quiver(T,Y, DT/M, DY/M, M, pivot='mid')plt.quiver(T,Y, DT/M, DY/M, # options de personnalisation M, cmap = plt.cm.viridis, scale_units='xy',scale=3. )plt.plot(tt,sol,'r', lw=4)plt.grid()plt.xlabel('t')plt.ylabel('y')plt.title('Champ des pentes');# Lignes de courant# =============================================================================plt.subplot(1,2,2)plt.streamplot(T,Y, DT/M, DY/M,# options de personnalisation# color=M, linewidth=1, cmap=plt.cm.viridis, density=2, arrowstyle='->', arrowsize=1.5 )plt.plot(tt,sol,'r', lw=4)plt.grid()plt.xlabel('t')plt.ylabel('y')plt.title('Lignes de courant');
1.3 ✨ Exemple de résolution approchée d’un système d’équations différentielles d’ordre 1 : modèle proie-prédateur
Considérons deux espèces : une proie (des lièvres par exemple) et un prédateur (des lynx par exemple). Ces deux populations sont représentées par \(y_1(t)\) et \(y_2(t)\) des fonctions continues du temps \(t\). Si on suppose qu’il n’y a aucune autre intervention extérieure, une modélisation possible pour ce genre de système a été proposée indépendamment par Alfred James Lotka en 1925 et Vito Volterra en 1926 :
\(
\begin{cases}
y_1'(t) = y_1(t)(a-by_2(t)) &[\stackrel{\text{déf}}{=}\varphi_1(t,y_1(t),y_2(t))] \quad\text{équation des proies}\\
y_2'(t) =-y_2(t)(c-dy_1(t)) &[\stackrel{\text{déf}}{=}\varphi_2(t,y_1(t),y_2(t))] \quad\text{équation des prédateurs}
\end{cases}
\)
\(y_1(0)=2\) unités de proies (une unités = \(1000\) animaux) et
\(y_2(0)=1\) unités de prédateurs.
On se demande comment vont évoluer les populations de ces deux espèces au cours du temps. Est-ce que les deux espèces vont coexister ou bien l’une des deux va disparaître ? Y a-t-il des cycles de croissance et de décroissance des populations ?
Pour les simulations on prendra \(a=2\), \(b=1\), \(c=1\) et \(d=0.3\).
%reset -f%matplotlib inlineimport matplotlib.pyplot as pltimport numpy as np# plt.rcdefaults()plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})from scipy.integrate import odeint
Code
# yy est une liste a deux composantespphi =lambda t,yy : [ yy[0]*(2-yy[1]) , -yy[1]*(1-0.3*yy[0]) ]
En faisant varier le temps sur l’intervalle \([0;20]\) et en prenant comme condition initiale le vecteur \(\mathbf{y}_0=(2,1)\) on écrit
Code
t0 =0yy0 = [2,1]tt = np.linspace(t0,20,201)# V1# odeint revoie une matrice, la colonne j contient l'approximation de la j-ème composante de la fonction inconnue# sol = odeint(pphi,yy0,tt,tfirst=True)# sol_1 = sol[:,0] # [y_1(t) for t in tt]# sol_2 = sol[:,1] # [y_2(t) for t in tt]# V2# Pour que odeint renvoit séparément les valeurs de la solution, il faut rajouter .T à la finsol_1,sol_2 = odeint(pphi,yy0,tt,tfirst=True).T
Le tracé des évolutions de \(y_1\) et \(y_2\) en fonction du temps \(t\) peut être obtenu par
Code
plt.figure(figsize=(10,7))plt.plot(tt,sol_1, label=r"$y_1(t)$ proies")plt.plot(tt,sol_2, label=r"$y_2(t)$ prédateurs")plt.grid()plt.xlabel('t')plt.ylabel('y')plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.);# astuce pour déplacer la légende à l'extérieur du plot
Le tracé de l’évolution de \(y_2\) en fonction de \(y_1\) peut être obtenu par
Considérons la fonction \(
E(t)=dy_1(t)+by_2(t)-c\ln(y_1(t))-a\ln(y_2(t))
\) Vérifions analytiquement que la fonction \(E\) est une intégrale première du système, c’est-à-dire que si \((y_1, y_2)\) est une solution alors l’application \(E\) est constante: \(\begin{align}
E'(t)
&=dy_1'(t)+by_2'(t)-\frac{c}{y_1(t)}y_1'(t)-\frac{a}{y_2(t)}y_2'(t)
\\
&=\left(d-\frac{c}{y_1(t)}\right)y_1'(t)+\left(b-\frac{a}{y_2(t)}\right)y_2'(t)
\\
&=\left(d-\frac{c}{y_1(t)}\right)y_1(t)(a-by_2(t))-\left(b-\frac{a}{y_2(t)}\right)y_2(t)(c-dy_1(t))
\\
&=(dy_1(t)- c)(a-by_2(t))-(by_2(t)-a)(c-dy_1(t))
\\
&=(dy_1(t)- c)(a-by_2(t))-(a-by_2(t))(dy_1(t)- c)
=0.
\end{align}\) Vérifions cette propriété numériquement:
Code
E =lambda yy : 0.3*yy[0]+yy[1]-np.log(yy[0])-2*np.log(yy[1])plt.figure(figsize=(10,7))E0 = E(yy0)EE = E([sol_1,sol_2])print(max(abs(EE-E0)))plt.hlines(E0,tt[0],tt[-1],colors='r',linestyles='dashed') # idem que plt.plot([tt[0],tt[-1]],[E0,E0],colors='r',linestyles='dashed') # valeur pour t=0plt.plot(tt,EE)plt.xlabel('t')plt.ylabel('E')plt.grid();
7.942278004158254e-07
1.4 ✨ Exemple de résolution approchée d’une équation différentielle d’ordre 2 après réduction à un système de 2 équations différentielles d’ordre 1 : Pendule simple
Considérons l’EDO \(y''(t)=-\sin(y(t))\) qui décrit le mouvement d’un pendule non amorti, équivalente au système différentiel
Considérons l’énergie suivante (somme de l’énergie potentielle et de l’énergie cinétique) : \(
E(t)=-\cos(y_1(t))+\frac{1}{2}(y_2(t))^2
\)
Vérifions analytiquement que la fonction \(E\) est conservée au cours du temps: \(\begin{align}
E'(t)
&=\sin(y_1(t))y_1'(t)+y_2(t)y_2'(t)
\\
&=-y_2'(t)y_1'(t)+y_1'(t)y_2'(t)
=0.
\end{align}\)
Vérifions cette propriété numériquement:
Code
E =lambda yy : -np.cos(yy[0])+yy[1]**2/2plt.figure(figsize=(10,7))E0 = E(yy0)EE = E([sol_1,sol_2])print(max(abs(EE-E0)))plt.hlines(E0,tt[0],tt[-1],colors='r',linestyles='dashed') # idem que plt.plot([tt[0],tt[-1]],[E0,E0],colors='r',linestyles='dashed') # valeur pour t=0plt.plot(tt,EE)plt.xlabel('t')plt.ylabel('E')plt.grid();
9.448590565508397e-08
2 ✍ Exercices
2.1 ✨ Exercice (EDO d’ordre 1)
Pour \(a=-2,-1,-\frac{1}{2},-\frac{1}{10},0,\frac{1}{10},\frac{1}{2},1,2\), afficher la solution approchée et le champ de vecteurs associé pour le problème de Cauchy
Afficher \(t\mapsto x\) et \(t\mapsto y\) dans un même graphe pour \(t\in[0,14]\).
Afficher dans un autre répère \(x\mapsto y\).
Soit une fonction \(t\mapsto E(t)\). Calculer analytiquement \(E'(t)\). Afficher ensuite \(E(t)\) calculé à partir des solutions approchées. La fonction \(E\) est définie par \(t\mapsto E(t)=\frac{x^2(t)}{2}+\frac{y^2(t)}{2}.\)
On vérifie ensuite numériquement que \(E\) est conservée (à une précision près) au cours du temps. Pour cela, on calcule la valeur exacte de \(E\) en \(t_0\), puis on trace l’évolution de \(E\) en fonction du temps. On constate que \(E\) est conservée à une \(10^{-7}\) près.
Code
# Formule qui définie l'énergie E =lambda yy : (yy[0]**2+yy[1]**2)/2# E initial exactE0 = E(yy0)# On evalue l'énergie en chaque point de la discrétisation en utlisant la solution approchéeEE = E([solx,soly])# On affiche le plus grand écart entre l'énergie initiale et l'énergie évaluée en chaque point# Si le calcul était exact, on devrait avoir max(abs(EE-E0)) = 0, norme infinie de l'erreur absolueprint(f"{</span><span class="bu">min</span>(EE<span class="op">-</span>E0) <span class="op">=</span> <span class="sc">:e}, {</span><span class="bu">max</span>(EE<span class="op">-</span>E0) <span class="op">=</span> <span class="sc">:e}") # Affichage# =============================================================================plt.figure(figsize=(17,7))# On affiche une droite qui correspond à l'énergie initialeplt.plot([tt[0],tt[-1]],[E0,E0]) # pour tracer le segment reliant (x0,y0) et (x1,y1) on peut écrire plt.plot([x0,x1],[y0,y1])# On affiche l'énergie évaluée en chaque point de la discrétisationplt.plot(tt,EE)plt.grid();
Afficher le champ des vecteurs et les solutions stationnaires.
Pour chaque point stationnaire, considérer 4 données initiales proches de ce point et calculer la solution particulière approchée associée à telle donnée. Que peut-on déduire sur la stabilité de ces points ?
L’expression du second membre est une fonction polynômiale en \(x\), \(y\), le théormème de Cauchy-Lipschitz s’applique donc et il existe une unique solution maximale pour chaque donnée de Cauchy.
Pour les données de Cauchy telles que \(
\begin{cases}
x^2+y^2=25\\
xy=12
\end{cases}
\)
la solution est globale et stationnaire. On cherche les couples \((x,y)\) qui satisfont ce système.
On peut par exemple retravailler le système comme suit: \(
\begin{cases}
x^2+y^2=25\\
xy=12
\end{cases}
\leadsto
\begin{cases}
(x+y)^2-2xy=25\\
xy=12
\end{cases}
\leadsto
\begin{cases}
(x+y)^2=49\\
xy=12
\end{cases}
\leadsto
\begin{cases}
x+y=\pm7\\
xy=12
\end{cases}
\leadsto
\begin{array}{l}
\begin{cases}
x+y=7\\
xy=12
\end{cases}
\leadsto
(3,4),(4,3)
\\
\begin{cases}
x+y=-7\\
xy=12
\end{cases}
\leadsto
(-3,-4),(-4,-3)
\end{array}
\)
soit résoludre graphiquemt en affichant les deux courbes du système. De telles données de Cauchy sont l’intersection du cercle centré en \(0\) et de rayon \(5\) et de l’hyperbole d’équation \(y =
12/x\). On obtient alors \(4\) points: \((3, 4)\), \((4, 3)\), \((-3, -4)\) et \((-4, -3)\).
Code
%reset -f%matplotlib inlineimport matplotlib.pyplot as pltimport numpy as npplt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})plt.figure(figsize=(10,7))# x^2+y^2=25 circle de centre (0,0) et rayon 5 donc d'équation x(t)=5cos(t) et y(t)=5sin(t)t = np.linspace(0,2*np.pi,101)plt.plot(5*np.cos(t),5*np.sin(t)) # pour afficher xy=12 on doit sauter 0 donc on utilise deux listes separéest = np.linspace(-5, -2, 101)plt.plot(t,12/t);t = np.linspace(2,5,101)plt.plot(t,12/t)plt.grid()plt.plot([3,4,-3,-4],[4,3,-4,-3],'o')plt.axis('equal');
