Calcul ✨ approché ✨ avec la fonction odeint du module SciPy

NotePlan de la séance de TP (partie 2)

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

# from matplotlib.pyplot import * # importe aussi numpy sans alias
import matplotlib.pyplot as plt
import numpy as np

# plt.rcdefaults()
plt. rcParams.update({'font.size': 16})
from scipy.integrate import odeint

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 :

scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)

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

Vous trouverez des informations détaillées dans la documentation officielle de SciPy: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

✨ Principe d’utilisation de odeint

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

\(\begin{cases} \mathbf{y}'(t)=\boldsymbol{\varphi}(t,\mathbf{y}(t)),\\ \mathbf{y}(t_0)=\mathbf{y}_0 \end{cases}\)

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))\).

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

Cf. https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Voici la signature de la fonction solve_ivp pour référence:

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)

Les arguments principaux sont :

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

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

Cf. https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.quiver.html et https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.streamplot.html

# Pour aller plus loin, décommenter un des deux appels suivants
# =============================================================================
# help(plt.quiver)
# help(plt.streamplot)

Exemples

✨ 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

\( y'(t)=\frac{3}{2} y(t) \left( 1-\frac{y(t)}{6} \right). \)

On cherche à afficher la solution approchée sur l’intervalle \([0,5]\).

On ajoutera ensuite le champ de vecteurs associé à cette équation différentielle.

On crée alors la fonction \(\varphi\)

La fonction odeint peut être appelée avec au minimum trois arguments :

• func: la fonction \(\varphi\),
• y0: la valeur \(y(t_0)\),
• t: le vecteur (qui commence à \(t_0\)) où la fonction \(y\) sera évaluée.

Elle renvoi un vecteur sol contenant l’évaluation de la solution en les points du vecteur \(t\).

Nota bene: par défaut la fonction \(\varphi\) doit avoir \(t\) comme deuxième variable, mais on retrouve l’ordre usuel avec le paramètre tfirst=True.

# def phi(t,y):
#     return 1.5*y*(1-y/6)
# ou 

phi = lambda t,y : 1.5*y*(1-y/6) # Attention, pas de calcul symbolique ici, on utlise des Float !

t0  = 0
y0  = 1
tt  = np.linspace(t0,5,21)
# sol = odeint(func=phi, y0=y0, t=tt, tfirst=True) # Appel avec arguments nommés: 4 paramètres, passés dans n'importe quel ordre
sol = odeint(phi, y0, tt, tfirst=True) # Appel classique: 4 paramètres, 3 passés selon l'ordre, le dernier nommé

# =============================================================================
# Bonus : on ajoute la comparaison avec les calculs de solve_ivp
# =============================================================================
from scipy.integrate import odeint, solve_ivp
sol_RK45   = solve_ivp(fun=phi, t_span=(t0,5), y0=[y0], t_eval=tt, method='RK45')
sol_RK23   = solve_ivp(fun=phi, t_span=(t0,5), y0=[y0], t_eval=tt, method='RK23')
sol_DOP853 = solve_ivp(fun=phi, t_span=(t0,5), y0=[y0], t_eval=tt, method='DOP853')
sol_Radau  = solve_ivp(fun=phi, t_span=(t0,5), y0=[y0], t_eval=tt, method='Radau')
sol_BDF    = solve_ivp(fun=phi, t_span=(t0,5), y0=[y0], t_eval=tt, method='BDF')
sol_LSODA  = solve_ivp(fun=phi, t_span=(t0,5), y0=[y0], t_eval=tt, method='LSODA')

La solution peut alors être tracée simplement

plt.figure(figsize=(10,7))
plt.plot(tt,sol,label='odeint')
plt.xlabel('t')
plt.ylabel('y')
plt.grid();


# =============================================================================
# Bonus
# =============================================================================
plt.plot(tt, sol_RK45.y[0],   'r--', label='RK45')
plt.plot(tt, sol_RK23.y[0],   'g--', label='RK23')
plt.plot(tt, sol_DOP853.y[0], 'c--', label='DOP853')
plt.plot(tt, sol_Radau.y[0],  'b--', label='Radau')
plt.plot(tt, sol_BDF.y[0],    'y--', label='BDF')
plt.plot(tt, sol_LSODA.y[0],  'm--', label='LSODA')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.);

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.

# 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 point
DT, 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');

✨ 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} \)

soit encore, en notation matricielle,

\(\begin{pmatrix} y_1\\ y_2 \end{pmatrix}'(t) = \begin{pmatrix} \varphi_1(t,y_1(t),y_2(t))\\ \varphi_2(t,y_1(t),y_2(t)) \end{pmatrix}.\)

On suppose qu’à l’instant initiale \(t=0\) il y a

• \(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\).

Pour une introduction à la modélisation de la dynamique de populations, voir par exemple cette page de vulgarisation Le modèle de Lotka-Volterra - Le lynx prédateur du lièvre

Lotka Volterra

Soit \(\mathbf{y}(t)\stackrel{\text{déf}}{=}(y_1(t),y_2(t))\) le vecteur des deux fonctions inconnues. On commence par créer la fonction vectorielle

\( \boldsymbol{\varphi}(t,\mathbf{y})=(\varphi_1(t,\mathbf{y}),\varphi_2(t,\mathbf{y})) \)

import matplotlib.pyplot as plt
import numpy as np
# plt.rcdefaults()
plt.rcParams.update({'font.size': 16})
from scipy.integrate import odeint

# yy est une liste a deux composantes
pphi = 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

t0  = 0
yy0 = [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 fin
sol_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

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

plt.figure(figsize=(10,7))
plt.plot(sol_1,sol_2)
plt.grid()
plt.xlabel(r'$y_1$')
plt.ylabel(r'$y_2$');

On peut superposer les deux graphes sur le plan des phases pour voir comment les solutions se comportent.

Attention, cette fois-ci les solutions sont des courbes paramétrées par \(t\) et non des fonctions de \(t\).

Y1,Y2 = np.meshgrid(np.linspace(min(sol_1),max(sol_1),21),np.linspace(min(sol_2),max(sol_2),21))
V1,V2 = pphi(tt,[Y1,Y2]) 
r1 = np.sqrt(1+V1**2)
r2 = np.sqrt(1+V2**2)

plt.figure(figsize=(10,7))
plt.quiver(Y1, Y2, V1/r1, V2/r2, cmap=plt.cm.viridis, scale_units='xy',scale=4.)
plt.plot(sol_1,sol_2,lw=2)
plt.grid()
plt.xlabel(r'$y_1$')
plt.ylabel(r'$y_2$');

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{aligned} 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{aligned} \) Vérifions cette propriété numériquement:

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=0
plt.plot(tt,EE)
plt.xlabel('t')
plt.ylabel('E')
plt.grid();

7.942278004158254e-07

✨ 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

\( \begin{cases} y_1'(t)=y_2(t),\\ y_2'(t)=-\sin(y_1(t)) \end{cases} \)

avec \(y(0)=0\) et \(y'(0)=1\).

On cherche à approcher la solution sur l’intervalle \([0,10]\).

import matplotlib.pyplot as plt
import numpy as np
# plt.rcdefaults()
plt.rcParams.update({'font.size': 16})
from scipy.integrate import odeint

pphi = lambda t,yy : [ yy[1], -np.sin(yy[0])  ]

t0  = 0
yy0  = [0,1]
tt  = np.linspace(t0,10,1001)
sol_1, sol_2 = odeint(pphi,yy0,tt,tfirst=True).T

plt.figure(figsize=(18,7))

plt.subplot(1,2,1)
plt.plot(tt,sol_1,tt,sol_2)
plt.xlabel('t')
plt.ylabel('y')
plt.grid()

plt.subplot(1,2,2)
plt.plot(sol_1,sol_2)
plt.axis('square')
plt.xlabel(r'$y_1$')
plt.ylabel(r'$y_2$')
plt.grid()
plt.axis('equal');

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{aligned} 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{aligned} \)

Vérifions cette propriété numériquement:

E = lambda yy : -np.cos(yy[0])+yy[1]**2/2

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=0
plt.plot(tt,EE)
plt.xlabel('t')
plt.ylabel('E')
plt.grid(); 

9.448590565508397e-08

✍ Exercices

✨ 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

\( \begin{cases} y'(t)=\sin\left(t y(t)\right), & t>0\\ y(0)=a \end{cases} \)

Bonus : savez-vous calculer la solution exacte en fonction de \(a\) ? Qu’en pense sympy ?

import matplotlib.pyplot as plt
import numpy as np 
plt.rcParams.update({'font.size': 14})
from scipy.integrate import odeint

phi = lambda t,y : np.sin(t*y)
t_final = 10
tt  = np.linspace(0,t_final,201)

plt.figure(figsize=(15,7)) 

# Champ des pentes / Lignes de courant
# À décommenter pour voir le champ des pentes au-dessus des solutions approchées
# =============================================================================
T,Y = np.meshgrid(np.linspace(0,t_final,61),np.linspace(-3,3,31))
DT, DY = np.ones_like(T), phi(T,Y) 
M = np.sqrt(DT**2+DY**2)   # norm growth rate 
# plt.quiver(T,Y, DT/M, DY/M, M, cmap=plt.cm.viridis, scale_units='xy',scale=4.)
plt.streamplot(T,Y, DT, DY,density=3)

plt.grid()
for y0 in [-2, -1, -0.5, -0.1, 0, 0.1, 0.5, 1, 2]:
    sol = odeint(phi, y0, tt, tfirst=True)
    plt.plot(tt, sol, lw=4, label=f'$y_0={y0}$');
plt.legend(bbox_to_anchor=(1.04,1),loc='upper left'); # Astuce pour déplacer la légende à l'extérieur du plot

# Affichage des deux graphes côte à côte
# =======================================

# T,Y = np.meshgrid(np.linspace(0,t_final,31),np.linspace(-3,3,31))
# DT, DY = np.ones_like(T), phi(T,Y) 
# M = np.sqrt(DT**2+DY**2)   # norm growth rate 

# fig, aax = plt.subplots(1,2,figsize=(20,7)) 


# aax[0].quiver(T,Y, DT/M, DY/M, M, pivot='mid')
# aax[0].grid()

# aax[1].streamplot(T,Y, DT/M, DY/M, color=M,density=2)
# aax[1].grid()

# for y0 in [-2, -1, -0.5, -0.1, 0, 0.1, 0.5, 1, 2]:
#     sol = odeint(phi,y0,tt,tfirst=True)
#     for ax in aax:
#         ax.plot(tt,sol,lw=4,label=f'$y_0={y0}$');
# plt.legend(bbox_to_anchor=(1.04,1),loc='upper left');

Malheureusement sympy ne peut pas nous fournir une expression analytique de la solution exacte en fonction de \(a\).

# import sympy as sym
# sym.init_printing()

# # Variables
# A = sym.symbols('A')
# t = sym.symbols('t', positive=True)
# y = sym.Function('y')

# # EDO
# edo = sym.Eq( sym.diff( y(t),t) , sym.sin(t*y(t)))
# display(edo)

# # Condition initiale
# t0 = 0
# y0 = A
# display(sym.Eq(y(t0), y0))

# # Méthode 1 : méthode directe
# # ============================
# solpar = sym.dsolve(edo,y(t), ics={y(t0):y0}) 
# display(solpar)

✨ Exercice (système)

Calculer la solution particulière approchée du système

\(\begin{cases} x'(t)= -y(t)+y^2(t)x(t)\\ y'(t)= x(t)-x^2(t)y(t) \end{cases} \)

avec \(x(0)=0\) et \(y(0)=1\).

1. Afficher \(t\mapsto x\) et \(t\mapsto y\) dans un même graphe pour \(t\in[0,14]\).
   
2. Afficher dans un autre répère \(x\mapsto y\).
3. 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}.\)

import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({'font.size': 16})
from scipy.integrate import odeint

# phi : R x R^2 -> R^2
pphi = lambda t,yy : [ -yy[1]+(yy[1])**2*yy[0] , yy[0]-(yy[0])**2*yy[1] ]

yy0  = [0,1]
tt   = np.linspace(0, 14, 201)
solx, soly = odeint(pphi, yy0, tt, tfirst=True).T

plt.figure(figsize=(17,7))

plt.subplot(1,2,1)
plt.plot(tt,solx,tt,soly)
plt.grid()
plt.xlabel(r'$t$')
plt.legend(["x","y"]);

plt.subplot(1,2,2) 
plt.plot(solx,soly)
plt.grid()
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.axis('equal');

On calcule la dérivée de \(E\) (analytiquement à la main ou, si on veut, avec sympy) pour vérifier que c’est une intégrale première du système :

\( E'(t)=x'(t)x(t)+y'(t)y(t)=\left(-y(t)+y^2(t)x(t)\right)x(t)+\left(x(t)-x^2(t)y(t)\right)y(t)=0. \)

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.

# Formule qui définie l'énergie 
E = lambda yy : (yy[0]**2+yy[1]**2)/2

# E initial exact
E0 = E(yy0)

# On evalue l'énergie en chaque point de la discrétisation en utlisant la solution approchée
EE = 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 absolue
print(f"{min(EE-E0) = :e}, {max(EE-E0) = :e}") 

# Affichage
# =============================================================================
plt.figure(figsize=(17,7))
# On affiche une droite qui correspond à l'énergie initiale
plt.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étisation
plt.plot(tt,EE)
plt.grid();

min(EE-E0) = -1.018167e-08, max(EE-E0) = 2.582059e-07

✨ Exercice (quiver - système )

Considérons le système

\( \begin{cases} x'(t) = x(t)(y(t) - x(t)^2) \\ y'(t) = y(t)(x(t) - y(t)^2) \end{cases} \)

Afficher le champ des vecteurs associé à ce système. Que remarquez-vous ?

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

pphi = lambda t,yy : [ yy[0]*(yy[1]-yy[0]**2) , yy[1]*(yy[0]-yy[1]**2) ]


# Preparation de la grille
# =============================================================================
t0  = 0
tt  = np.linspace(t0,20,201)
Y1,Y2 = np.meshgrid(np.linspace(-2,2,51),np.linspace(-2,2,51))
V1,V2 = pphi(tt,[Y1,Y2]) 
M = np.hypot(V1,V2)
V1 = V1/np.sqrt(V1**2+V2**2+1)
V2 = V2/np.sqrt(V1**2+V2**2+1)

# Affichage
# =============================================================================
plt.figure(figsize=(10,7))
# plt.quiver(Y1, Y2, V1, V2,
#         #    M, cmap = plt.cm.viridis, scale_units='xy',scale=10.
#         )
plt.streamplot(Y1, Y2, V1, V2,
               # options de personnalisation
               color=M, linewidth=1, cmap=plt.cm.viridis, density=2, arrowstyle='->', arrowsize=1.5
               )
plt.grid()
plt.xlabel(r'$y_1$')
plt.ylabel(r'$y_2$')

# Points stationnaires
# =============================================================================
plt.plot(0,0,'o',ms=10)
plt.plot(1,1,'o',ms=10);

# Une solution approchée
# =============================================================================
yy0 = [ 0.1,1]
# yy0 = [-0.1,1]
sol_1,sol_2 = odeint(pphi,yy0,tt,tfirst=True).T
plt.plot(sol_1,sol_2,lw=5)

# Courbes 
# =============================================================================
x_tmp = np.linspace(0,1,101)
plt.plot(x_tmp,x_tmp**2, 'r')
plt.plot(x_tmp**2,x_tmp, 'r');

✨ Exercice (système) : étude de la stabilité

Considérons le système

\( \begin{cases} x'(t)= x^2(t)+y^2(t)-25\\ y'(t)= x(t)y(t)-12 \end{cases} \)

1. Afficher le champ des vecteurs et les solutions stationnaires.
2. 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)\).

%reset -f
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({'font.size': 16})

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ées
t = 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');

• soit résoudre \( \begin{cases} x^2+y^2=25\\ xy=12 \end{cases} \stackrel{x\neq0}{\leadsto} \begin{cases} x^2+y^2=25\\ y=\frac{12}{x} \end{cases} \leadsto \begin{cases} x^2+\frac{12^2}{x^2}=25\\ y=\frac{12}{x} \end{cases} \leadsto \begin{cases} x^4-25x^2+12^2=0\\ y=\frac{12}{x} \end{cases} \)

On peut ainsi calculer les valeurs de \(x\) qui satisfont la première équation et en déduire les valeurs de \(y\) associées.

# numpy calcule les racines d'un polynome
x_sol = np.roots([1, 0, -25, 0, 12**2])
y_sol = 12/x_sol
print(x_sol,y_sol)

[-4. -3.  4.  3.] [-3. -4.  3.  4.]

En conclusion, on a 4 points stationnaires: \( \{(3,4) , (-4,-3), (4,3) , (-3,-4)\} \)

from scipy.integrate import odeint

pphi = lambda t,y : [ (y[0])**2+(y[1])**2-25 , y[0]*y[1]-12 ]

Y1,Y2 = np.meshgrid(np.linspace(-6,6,21),np.linspace(-6,6,21))
tt    = np.linspace(0,0.2,201)
V1,V2 = pphi(tt,[Y1,Y2])
r = np.sqrt(V1**2+V2**2)
# r[ r==0 ] = 1

plt.figure(figsize=(20,7))

plt.subplot(1,2,1)
# plt.quiver(Y1, Y2, V1/r, V2/r, r, cmap=plt.cm.viridis)
plt.quiver(Y1, Y2, V1, V2, r, cmap=plt.cm.viridis, scale_units='xy',scale=20.)
plt.plot(x_sol,y_sol,'o',ms=10)
plt.xticks(range(-6,6))
plt.yticks(range(-6,6))
plt.grid()
plt.xlabel(r'$y_1$')
plt.ylabel(r'$y_2$');

plt.subplot(1,2,2)
# plt.streamplot(Y1, Y2, V1/r, V2/r) 
plt.streamplot(Y1, Y2, V1/r, V2/r, color=r, linewidth=1, cmap=plt.cm.viridis, density=2, arrowstyle='->', arrowsize=1.5)
plt.plot(x_sol,y_sol,'o',ms=10)
plt.xticks(range(-6,6))
plt.yticks(range(-6,6))
plt.grid()
plt.xlabel(r'$y_1$')
plt.ylabel(r'$y_2$');

plt.figure(figsize=(17,14))

plt.subplot(2,2,1)
plt.title('Voisinage de (3,4)')
plt.plot(3,4,'rD',ms=20)  
plt.grid()
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
eps = 1
for y0,tfin in [ ((3+0.1,4+0.1),0.2), ((3-0.1,4+0.1),0.5) , ((3+0.1,4-0.1),0.4), ((3-0.1,4-0.1),0.2)]:
    tt  = np.linspace(0,tfin,201)
    solx1, soly1 = odeint(pphi,y0,tt,tfirst=True).T
    plt.plot(solx1,soly1,"r-",lw=2)
    plt.plot(solx1[0],soly1[0],'ro',ms=10) # point de départ proche de  (3,4)
    Y1,Y2 = np.meshgrid(np.linspace(3-eps,3+eps,11),np.linspace(4-eps,4+eps,11))
    V1,V2 = pphi(tt,[Y1,Y2])
    plt.streamplot(Y1, Y2, V1, V2, color='b',linewidth=0.1)


plt.subplot(2,2,2)
plt.title('Voisinage de (-4,-3)')
plt.plot(-4,-3,'D',c='purple',ms=20) 
plt.grid()
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
eps = 0.2
for y0,tfin in [ ((-4+0.1,-3+0.1),0.2), ((-4-0.1,-3+0.1),0.4) , ((-4+0.1,-3-0.1),0.4), ((-4-0.1,-3-0.1),0.2)]:
    tt  = np.linspace(0,tfin,201)
    solx1,soly1 = odeint(pphi,y0,tt,tfirst=True).T
    plt.plot(solx1,soly1,"-",c='purple',lw=2) 
    plt.plot(solx1[0],soly1[0],'o',c='purple',ms=10) # point de départ proche de  (-4,-3)
    Y1,Y2 = np.meshgrid(np.linspace(-4-eps,-4+eps,11),np.linspace(-3-eps,-3+eps,11))
    V1,V2 = pphi(tt,[Y1,Y2])
    plt.streamplot(Y1, Y2, V1, V2, color='b',linewidth=0.1)


plt.subplot(2,2,3)
plt.title('Voisinage de (4,3)')
plt.plot(4,3,'mD',ms=20) 
plt.grid()
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
eps = 1
for y0,tfin in [ ((4+0.1,3+0.1),0.2), ((4-0.1,3+0.1),0.4) , ((4+0.1,3-0.1),0.4), ((4-0.1,3-0.1),0.2)]:
    tt  = np.linspace(0,tfin,201)
    solx1, soly1 = odeint(pphi,y0,tt,tfirst=True).T
    plt.plot(solx1,soly1,"m-",lw=2) 
    plt.plot(solx1[0],soly1[0],'mo',ms=10) # point de départ proche de   (4,3) 
    Y1,Y2 = np.meshgrid(np.linspace(4-eps,4+eps,11),np.linspace(3-eps,3+eps,11))
    V1,V2 = pphi(tt,[Y1,Y2])
    plt.streamplot(Y1, Y2, V1, V2, color='b',linewidth=0.1)

    
plt.subplot(2,2,4)
plt.title('Voisinage de (-3,-4)')
plt.plot(-3,-4,'gD',ms=20) 
plt.grid()
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
eps = 0.2
for y0,tfin in [ ((-3+0.1,-4+0.1),0.2), ((-3-0.1,-4+0.1),0.4) , ((-3+0.1,-4-0.1),0.4), ((-3-0.1,-4-0.1),0.2)]:
    tt  = np.linspace(0,tfin,201)
    solx1, soly1 = odeint(pphi,y0,tt,tfirst=True).T
    plt.plot(solx1,soly1,"g-",lw=2) 
    plt.plot(solx1[0],soly1[0],'go',ms=10); # point de départ proche de (-3,-4)  
    Y1,Y2 = np.meshgrid(np.linspace(-3-eps,-3+eps,11),np.linspace(-4-eps,-4+eps,11))
    V1,V2 = pphi(tt,[Y1,Y2])
    plt.streamplot(Y1, Y2, V1, V2, color='b',linewidth=0.1)