from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())

import sys #only needed to determine Python version number
print(f'Python version: {sys.version}')

import numpy as np #only needed to determine Numpy version number
print(f"Numpy version: {np.__version__}")

from datetime import datetime
print('Last Updated On: ', datetime.now())

Python version: 3.10.6 (main, Nov 14 2022, 16:10:14) [GCC 11.3.0]
Numpy version: 1.21.5
Last Updated On:  2023-02-21 13:20:58.761298

M62_TP1 EDO : calcul approché VS formel.¶

1 ✨ Calcul approché avec la fonction odeint du module SciPy
2 🐍 Calcul analytique avec le module sympy
- 2.1 Exemple de résolution analytique d’une équation différentielle d'ordre 1
- 2.2 Exemple de résolution analytique d’un système d'équations différentielles d'ordre 1
3 ✍ Exercices

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

L'ECUE M62 n'est qu'une introduction aux méthodes numériques d'approximation d'EDO en utilisant Python comme langage commun de programmation. Nous allons étudier, programmer et tester plusieurs méthodes numériques.

Cependant toutes ces méthodes (et bien plus) se trouvent déjà dans le module SciPy.
Bien sûr vous devez toujours faire un peu de programmation avec Python pour les utiliser et une compréhension des fondements de la méthode numérique que vous utilisez est toujours indispensables.

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

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 $$\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 forme quelconque, 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).

Cf. https://docs.scipy.org/doc/scipy-1.2.1/reference/generated/scipy.integrate.odeint.html

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

Remarque: 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-1.2.1/reference/generated/scipy.integrate.solve_ivp.html

Exemple de résolution approchée d’une équation différentielle d'ordre 1¶

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 crée alors la fonction $\varphi$

%reset -f
%matplotlib inline
from matplotlib.pylab import * # importe aussi numpy sans alias
# rcdefaults()
rcParams.update({'font.size': 16})
from scipy.integrate import odeint

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)

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

La solution peut alors être tracée simplement

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

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

Champ de vecteurs et lignes de courant¶

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 matplotlib permet de tracer un champ de vecteurs: quiver(x_pos, y_pos, x_dir, y_dir, color)
La fonction streamplot 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.

# define a grid 
g1  = linspace(0,5,21)
g2  = linspace(0,8,21)
T,Y = meshgrid(g1,g2)         # create a grid

# compute direction at each point
DT, DY = 1, phi(T,Y)          # compute growth rate on the grid
# M = sqrt(DT**2+DY**2)         # norm growth rate 
M = hypot(DT,DY)

figure(figsize=(20,7))

# Champ des pentes
subplot(1,2,1)
# quiver(T,Y, DT/M, DY/M, M, pivot='mid')
quiver(T,Y, DT/M, DY/M, M, cmap=cm.viridis, scale_units='xy',scale=3.)
# https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.quiver.html
plot(tt,sol,'r', lw=4)
grid()
xlabel('t')
ylabel('y')
title('Champ des pentes');

# Lignes de courant
subplot(1,2,2)
# streamplot(T,Y, DT/M, DY/M, color=M)
streamplot(T,Y, DT/M, DY/M, color=M, linewidth=1, cmap=cm.viridis, density=2, arrowstyle='->', arrowsize=1.5)
# cf. https://matplotlib.org/stable/tutorials/colors/colormaps.html
plot(tt,sol,'r', lw=4)
grid()
xlabel('t')
ylabel('y')
title('Lignes de courant');

Exemple de résolution approchée d’un système d'équations différentielles d'ordre 1¶

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érieur, 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(y_1(t),y_2(t),t)] \quad\text{équation équation des proies}\\ y_2'(t) =-y_2(t)(c-dy_1(t)) &[\stackrel{\text{déf}}{=}\varphi_2(y_1(t),y_2(t),t)] \quad\text{équation équation des prédateurs} \end{cases} $$ soit encore, en notation matricielle,

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

On suppose qu’à ce jour 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 et on se demande comment vont évoluer les populations de ces deux espèces. Pour les simulations on prendra $a=2$, $b=1$, $c=1$ et $d=0.3$.

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}(\mathbf{y},t)=(\varphi_1(\mathbf{y},t),\varphi_2(\mathbf{y},t)) $$

%reset -f
%matplotlib inline
from matplotlib.pylab import * # importe aussi numpy sans alias
# rcdefaults()
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  = 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

figure(figsize=(10,7))
plot(tt,sol_1,tt,sol_2)
grid()
xlabel('t')
ylabel('y')
legend([r"$y_1(t)$ proies","$y_2(t)$ prédateurs"]);

Le tracé des évolutions de $y_2$ en fonction de $y_1$ peut être obtenu par

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

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

figure(figsize=(10,7))
quiver(Y1, Y2, V1/r1, V2/r2, cmap=cm.viridis, scale_units='xy',scale=4.)
plot(sol_1,sol_2)
grid()
xlabel(r'$y_1$')
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{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:

figure(figsize=(10,7))
# valeur pour t=0
E0 = 0.3*yy0[0]+yy0[1]-log(yy0[0])-2*log(yy0[1])
plot([tt[0],tt[-1]],[E0,E0])

EE  = 0.3*sol_1+sol_2-log(sol_1)-2*log(sol_2)
print(max(abs(EE-E0)))
plot(tt,EE)
xlabel('t')
ylabel('E')
grid();

7.942278004158254e-07