In [1]:
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
Out[1]:
In [2]:
from datetime import datetime
print('Last Updated On: ', datetime.now())
Last Updated On:  2022-03-08 10:51:41.128079

M62_TP1 EDO : calcul approché VS formel.

✨ 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, avec $t$ comme deuxième paramètre), la condition initiale $\mathbf{y}_0$ et le domaine de temps qui nous intéresse (qui commence à $t_0$). Elle retourne un tableau Numpy (même si $t$ était une liste).

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

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$

In [3]:
%reset -f
%matplotlib inline
from matplotlib.pylab import * # importe aussi numpy sans alias
from scipy.integrate import odeint
In [4]:
# def phi_odeint(y,t):
#     return 1.5*y*(1-y/6)
# ou 

# Attention à l'ordre des variables dans phi
phi_odeint = lambda y,t : 1.5*y*(1-y/6)

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

  • la fonction $\varphi$,
  • la valeur $y(t_0)$,
  • le vecteur $t$ (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$.

In [5]:
t0  = 0
y0  = 1
tt  = linspace(t0,5,201)
sol = odeint(phi_odeint,y0,tt)

La solution peut alors être tracée simplement

In [6]:
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

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. Utilisons-la pour obtenir celui associé à notre équation différentielle. La courbe rouge est la solution déterminée par odeint.

In [7]:
figure(figsize=(10,7))

# quiverplot
# define a grid and compute direction at each point
g1  = linspace(0,5,21)
g2  = linspace(0,8,21)
T,Y = meshgrid(g1,g2)         # create a grid
DT, DY = 1, phi_odeint(Y,T)   # compute growth rate on the grid
M = sqrt(DT**2+DY**2)         # norm growth rate 
M[ M==0 ] = 1                 # avoid zero division errors 
quiver(T,Y, DT/M, DY/M, M, pivot='mid')

plot(tt,sol,'r')
grid()
xlabel('t')
ylabel('y')
title('Champ des pentes');

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

In [8]:
%reset -f
%matplotlib inline
from matplotlib.pylab import * # importe aussi numpy sans alias
from scipy.integrate import odeint
In [9]:
# yy est une liste a deux composantes
pphi_odeint = lambda yy,t : [ 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

In [10]:
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_odeint,yy0,tt)
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_odeint,yy0,tt).T

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

In [11]:
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

In [12]:
figure(figsize=(10,7))
plot(sol_1,sol_2)
grid()
xlabel(r'$y_1$')
ylabel(r'$y_2$');
In [13]:
figure(figsize=(10,7))
Y1,Y2 = meshgrid(linspace(min(sol_1),max(sol_1),21),linspace(min(sol_2),max(sol_2),21))
V1,V2 = pphi_odeint([Y1,Y2],tt) 
r1=sqrt(1+V1**2)
r2=sqrt(1+V2**2)
quiver(Y1, Y2, V1/r1, V2/r2) 
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:

In [14]:
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

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

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

In [15]:
%reset -f
%matplotlib inline
from matplotlib.pylab import * # importe aussi numpy sans alias
from scipy.integrate import odeint
In [17]:
pphi_odeint = lambda yy,t : [ yy[1], -sin(yy[0])  ]

t0  = 0
yy0  = [0,1]
tt  = linspace(t0,10,1001)
sol_1, sol_2 = odeint(pphi_odeint,yy0,tt).T

figure(figsize=(18,7))

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

subplot(1,2,2)
plot(sol_1,sol_2)
xlabel(r'$y_1$')
ylabel(r'$y_2$')
grid()
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{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:

In [18]:
figure(figsize=(10,7))
# valeur pour t=0
E0 = -cos(yy0[0])+yy0[1]**2/2
plot([tt[0],tt[-1]],[E0,E0])

EE  = -cos(sol_1)+sol_2**2/2
print(max(abs(EE-E0)))
plot(tt,EE)
xlabel('t')
ylabel('E')
grid(); 
9.448590565508397e-08

🐍 Calcul analytique avec le module sympy


https://docs.sympy.org/latest/modules/solvers/ode.html

In [19]:
%reset -f
%matplotlib inline

import sympy as sym
sym.init_printing()

Exemple de résolution analytique d’une équation différentielle d'ordre 1

Considérons le problème de Cauchy $$\begin{cases} u'(x)=-3x^2u(x)+6x^2,\\ u(0)=4. \end{cases}$$

In [20]:
x   = sym.Symbol('x')
u   = sym.Function('u')
left  = sym.diff(u(x),x)
right = 6*x**2-3*x**2*u(x)
edo = sym.Eq( left , right )
display(edo)
$\displaystyle \frac{d}{d x} u{\left(x \right)} = - 3 x^{2} u{\left(x \right)} + 6 x^{2}$
In [21]:
solgen = sym.dsolve(edo,u(x))
display(solgen)
$\displaystyle u{\left(x \right)} = C_{1} e^{- x^{3}} + 2$

Prise en compte des conditions initiales:

In [22]:
x0=0
u0=4
consts = sym.solve( sym.Eq( u0, solgen.rhs.subs(x,x0)) , dict=True)[0]
consts
Out[22]:
$\displaystyle \left\{ C_{1} : 2\right\}$
In [23]:
solpar = solgen.subs(consts)
solpar
Out[23]:
$\displaystyle u{\left(x \right)} = 2 + 2 e^{- x^{3}}$

On transforme le résultat en une fonction vectorisée :

In [24]:
func = sym.lambdify(x,solpar.rhs,'numpy')
In [25]:
from matplotlib.pylab import *
figure(figsize=(10,7))
xx=linspace(0,3,101)
yy=func(xx)
plot(xx,yy)
xlabel('x')
ylabel('u')
grid();

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

$$ \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} = \begin{pmatrix} y_1(t)-y_2(t)\\ y_2(t)-y_1(t) \end{pmatrix}. $$

avec $t\in[0;10]$ et en prenant comme condition initiale le vecteur $\mathbf{y}_0=(2,1)$.

In [26]:
%reset -f
%matplotlib inline

import sympy as sym
sym.init_printing()

t  = sym.Symbol('t')

y1 = sym.Function('y1')
y2 = sym.Function('y2')

right1 = y1(t)-y2(t)
right2 = y2(t)-y1(t)

edo1 = sym.Eq( sym.diff(y1(t),t) , right1 )
edo2 = sym.Eq( sym.diff(y2(t),t) , right2 )
display(edo1)
display(edo2)

solgen = sym.dsolve([edo1,edo2],[y1(t),y2(t)])
display(solgen)
$\displaystyle \frac{d}{d t} \operatorname{y_{1}}{\left(t \right)} = \operatorname{y_{1}}{\left(t \right)} - \operatorname{y_{2}}{\left(t \right)}$
$\displaystyle \frac{d}{d t} \operatorname{y_{2}}{\left(t \right)} = - \operatorname{y_{1}}{\left(t \right)} + \operatorname{y_{2}}{\left(t \right)}$
$\displaystyle \left[ \operatorname{y_{1}}{\left(t \right)} = - C_{1} e^{2 t} - C_{2}, \ \operatorname{y_{2}}{\left(t \right)} = C_{1} e^{2 t} - C_{2}\right]$
In [27]:
t0  =0
y1_0=2
y2_0=1
consts = sym.solve( [sym.Eq( y1_0, solgen[0].rhs.subs(t,t0)) , 
                     sym.Eq( y2_0, solgen[1].rhs.subs(t,t0)) ] , dict=True)[0]
display(consts)

solpar_1=solgen[0].subs(consts)
solpar_2=solgen[1].subs(consts)
display(solpar_1)
display(solpar_2)

func_1 = sym.lambdify(t,solpar_1.rhs,'numpy')
func_2 = sym.lambdify(t,solpar_2.rhs,'numpy')

from matplotlib.pylab import *
figure(figsize=(18,7))

tt=linspace(0,3,101)
yy_1=func_1(tt)
yy_2=func_2(tt)

subplot(1,2,1)
plot(tt,yy_1,tt,yy_2)
legend([r'$y_1$',r'$y_2$'])
xlabel(r'$t$')
ylabel(r'$y$')

subplot(1,2,2)
plot(yy_1,yy_2)
xlabel(r'$y_1$')
ylabel(r'$y_2$');
$\displaystyle \left\{ C_{1} : - \frac{1}{2}, \ C_{2} : - \frac{3}{2}\right\}$
$\displaystyle \operatorname{y_{1}}{\left(t \right)} = \frac{e^{2 t}}{2} + \frac{3}{2}$
$\displaystyle \operatorname{y_{2}}{\left(t \right)} = \frac{3}{2} - \frac{e^{2 t}}{2}$

✍ Exercices

🐍 Exercice (sympy)

Calculer la solution exacte du problème de Cauchy $$\begin{cases} y'(t)=\dfrac{3}{2} y(t) \left( 1-\dfrac{y(t)}{6} \right),& t>0\\ y(0)=1. \end{cases}$$

In [28]:
%reset -f
%matplotlib inline

import sympy as sym
sym.init_printing()

t  = sym.Symbol('t')
y  = sym.Function('y')
right = sym.Rational(3,2)*y(t)*(1-y(t)/6)
edo = sym.Eq( sym.diff(y(t),t) , right )
display(edo)
solgenLIST = sym.dsolve(edo,y(t))
display(solgenLIST)
$\displaystyle \frac{d}{d t} y{\left(t \right)} = \frac{3 \left(1 - \frac{y{\left(t \right)}}{6}\right) y{\left(t \right)}}{2}$
$\displaystyle \left[ y{\left(t \right)} = \frac{6 \left(- e^{3 t} + \sqrt{e^{C_{1} + 3 t}}\right)}{e^{C_{1}} - e^{3 t}}, \ y{\left(t \right)} = - \frac{6 \left(e^{3 t} + \sqrt{e^{C_{1} + 3 t}}\right)}{e^{C_{1}} - e^{3 t}}\right]$
In [29]:
solgen = solgenLIST[0]
display(solgen)
$\displaystyle y{\left(t \right)} = \frac{6 \left(- e^{3 t} + \sqrt{e^{C_{1} + 3 t}}\right)}{e^{C_{1}} - e^{3 t}}$
In [30]:
t0 = 0 
y0 = 1
consts = sym.solve( sym.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
display(consts)

solpar = solgen.subs(consts)
display(solpar)

func = sym.lambdify(t,solpar.rhs,'numpy')

from matplotlib.pylab import *
figure(figsize=(10,7))
tt = linspace(0,3,101)
yy = func(tt)
plot(tt,yy)
grid();
$\displaystyle \left\{ C_{1} : \log{\left(25 \right)}\right\}$
$\displaystyle y{\left(t \right)} = \frac{6 \left(- e^{3 t} + 5 \sqrt{e^{3 t}}\right)}{25 - e^{3 t}}$

✨ Exercice (scipy)

Calculer la solution approchée du problème de Cauchy $$\begin{cases} y'(t)=\sin\left(t y(t)\right), & t>0\\ y(0)=a \end{cases}$$ pour $a=-2,-1,-\frac{1}{2},-\frac{1}{10},0,\frac{1}{10},\frac{1}{2},1,2$ et afficher le champ de vecteurs.

In [31]:
%reset -f
%matplotlib inline
from matplotlib.pylab import * # importe aussi numpy sans alias
from scipy.integrate import odeint

t_final = 10

phi_odeint = lambda y,t : sin(t*y)
tt  = linspace(0,t_final,201)

figure(figsize=(10,7)) # facultatif

T,Y = meshgrid(linspace(0,t_final,31),linspace(-3,3,31))
DT, DY = 1, phi_odeint(Y,T) 
M = sqrt(DT**2+DY**2)   # norm growth rate 
M[ M == 0] = 1          # avoid zero division errors 
quiver(T,Y, DT/M, DY/M, M, pivot='mid')
grid()

for y0 in [-2, -1, -0.5, -0.1, 0, 0.1, 0.5, 1, 2]:
    sol=odeint(phi_odeint,y0,tt)
    plot(tt,sol,lw=2,label='$y_0=$'+str(y0));
legend(bbox_to_anchor=(1.04,1),loc='upper left');