TP1 - Calcul approché VS calcul formel

import sys #only needed to determine Python version number
print(f'Python version: {</span>sys<span class="sc">.</span>version<span class="sc">}')

import numpy as np # only needed to determine Numpy version number
print(f"Numpy version: {</span>np<span class="sc">.</span>__version__<span class="sc">}")

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

Python version: 3.10.12 (main, Nov 20 2023, 15:14:05) [GCC 11.4.0]
Numpy version: 1.23.5
Last Updated On:  2024-03-22 18:34:05.853521

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

L’ECUE M62 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).

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

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 urilisée est LSODA (méthode Adam/BDF adaptative).

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

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

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

1.3 ✨ 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.pyplot import * # importe aussi numpy sans alias
import matplotlib.pyplot as plt
import numpy as np

# plt.rcdefaults()
plt. rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})
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  = 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).

1.3.1 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 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

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

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, 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, color=M)
plt.streamplot(T,Y, DT/M, DY/M, 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.4 ✨ 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(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’à 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}(t,\mathbf{y})=(\varphi_1(t,\mathbf{y}),\varphi_2(t,\mathbf{y})) \)

%reset -f
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
# plt.rcdefaults()
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})
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();

Le tracé des évolutions 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$');

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{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:

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

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

%reset -f
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
# plt.rcdefaults()
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})
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{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:

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

2 Calcul 🐍 analytique 🐍 avec le module sympy

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

%reset -f
%matplotlib inline

import sympy as sym
sym.init_printing()

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

Définition de l’edo :

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

Calcul de la solution exacte, Méthode pas-à-pas :

Étape 1 : calcul de la solution générale (qui contient la constante d’intégration) :

solgen = sym.dsolve(edo,u(x))
display(solgen)

\(\displaystyle u{\left(x \right)} = C_{1} e^{- x^{3}} + 2\)

Étape 2 : calcul de la constante d’intégration en prenant en compte la condition initiale :

x0 = 0
u0 = 4
cond_init = sym.Eq( u0, solgen.rhs.subs(x,x0) )
display(cond_init)
consts = sym.solve( cond_init , dict=True)[0]
display(consts)

\(\displaystyle 4 = C_{1} + 2\)

\(\displaystyle \left\{ C_{1} : 2\right\}\)

Étape 3 : calcul de la solution particulière :

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

\(\displaystyle u{\left(x \right)} = 2 + 2 e^{- x^{3}}\)

Calcul de la solution exacte, Méthode directe :

solpar = sym.dsolve(edo,u(x), ics={u(x0):u0})
display(solpar)

\(\displaystyle u{\left(x \right)} = 2 + 2 e^{- x^{3}}\)

Conersion en une fonction de numpy :

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

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

import matplotlib.pyplot as plt
import numpy as np
# plt.rcdefaults()
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})
plt.figure(figsize=(10,7))
xx = np.linspace(0,3,101)
yy = func(xx)
plt.plot(xx,yy)
plt.xlabel('x')
plt.ylabel('u')
plt.grid();

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

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

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

t0   = 0
y1_0 = 2
y2_0 = 1

# solgen = sym.dsolve([edo1,edo2],[y1(t),y2(t)])
# display(solgen)
# 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)

solpar_1, solpar_2 = sym.dsolve( [ edo1, edo2 ] , ics = { y1(t0):y1_0, y2(t0):y2_0 } )
display(solpar_1)
display(solpar_2)

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

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

import matplotlib.pyplot as plt
import numpy as np
# plt.rcdefaults()
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})
plt.figure(figsize=(18,7))

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

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

plt.subplot(1,2,2)
plt.plot(yy_1,yy_2)
plt.xlabel(r'$y_1$')
plt.ylabel(r'$y_2$');

3 ✍ Exercices

3.1 🐍 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} \)

Que vaut \(\lim_{t\to\infty}y(t)\)?

%reset -f
%matplotlib inline

import sympy as sym
sym.init_printing()

t0 = 0 
y0 = 1


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)

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

# Méthode pas-à-pas
# solgenLIST = sym.dsolve(edo,y(t))
# display(solgenLIST)

# solgen = solgenLIST[0]
# display(solgen)

# equation = sym.Eq( y0, solgen.rhs.subs(t,t0))
# display(equation)
# consts = sym.solve( equation , dict=True)[0]
# display(consts)

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

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

\(\displaystyle y{\left(t \right)} = \frac{6 \left(- e^{3 t} + 5 \sqrt{e^{3 t}}\right)}{25 - e^{3 t}}\)

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

import matplotlib.pyplot as plt 
import numpy as np
# plt.rcdefaults()
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})
plt.figure(figsize=(10,7))
tt = np.linspace(0,10,101)
yy = func(tt)
plt.plot(tt,yy)
plt.grid();

limite = sym.limit(solpar.rhs,t,sym.oo)

# Affichage amélioré
from IPython.core.display import Math
display( Math(  r"\lim_{t\to+\infty}y(t)=" + sym.latex(limite)  ) )

\(\displaystyle \lim_{t\to+\infty}y(t)=6\)

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

%reset -f
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np 
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})
from scipy.integrate import odeint

t_final = 10

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

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

# 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 
# 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={</span>y0<span class="sc">}$');
plt.legend(bbox_to_anchor=(1.04,1),loc='upper left');

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

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

3.3 🐍 Exercice (sympy - système)

1. Calculer la solution générale exacte du système
   \( \begin{pmatrix} x(t)\\ y(t) \end{pmatrix}' = \begin{pmatrix} 0&1\\ -1&0 \end{pmatrix} \begin{pmatrix} x(t)\\ y(t) \end{pmatrix} \)
   
2. Pour \(x(0)=-1\) et \(y(0)=1\), afficher \(t\mapsto x\) et \(t\mapsto y\) dans un même graphe.
3. Afficher ensuite \(x\mapsto y\).
4. Calculer analytiquement \(E'(t)\) avec \( E(t)=\frac{x^2(t)}{2}+\frac{y^2(t)}{2}. \)

%reset -f
%matplotlib inline

import sympy as sym
sym.init_printing()

t_0 = 0
x_0 = -1
y_0 = 1

t = sym.Symbol('t')
x = sym.Function('x')
y = sym.Function('y')
edo1 = sym.Eq( sym.diff(x(t),t) ,  y(t) )
edo2 = sym.Eq( sym.diff(y(t),t) , -x(t) )
display(edo1)
display(edo2)

\(\displaystyle \frac{d}{d t} x{\left(t \right)} = y{\left(t \right)}\)

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = - x{\left(t \right)}\)

# Méthode pas-à-pas
# solgen = sym.dsolve([edo1,edo2],[x(t),y(t)])
# display(solgen)

# consts = sym.solve( [ sym.Eq( x_0, solgen[0].rhs.subs(t,t_0)) , sym.Eq( y_0, solgen[1].rhs.subs(t,t_0)) ] , dict=True)[0]
# display(consts)

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

# Méthode compacte
solpar_1, solpar_2 = sym.dsolve( [edo1,edo2],[x(t),y(t)], ics={x(t_0):x_0, y(t_0):y_0} )
display(solpar_1)
display(solpar_2)

\(\displaystyle x{\left(t \right)} = \sin{\left(t \right)} - \cos{\left(t \right)}\)

\(\displaystyle y{\left(t \right)} = \sin{\left(t \right)} + \cos{\left(t \right)}\)

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

import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})
plt.figure(figsize=(17,7))
tt = np.linspace(0,30,101)
xx = func_1(tt)
yy = func_2(tt)

plt.subplot(1,2,1)
plt.plot(tt,xx,tt,yy)
plt.legend([r'$x(t)$',r'$y(t)$'])
plt.xlabel(r'$x$')
plt.grid()

plt.subplot(1,2,2)
plt.plot(xx,yy)
# plt.axis('square')
plt.axis('equal')
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.grid();

E = (solpar_1.rhs)**2/2+(solpar_2.rhs)**2/2
# display(E)
Ep = sym.diff(E,t).simplify()
# display(Ep)

# Affichage amélioré
from IPython.core.display import Math
display( Math(  "E(t)=" + sym.latex(E)  ) )
display( Math(  "E'(t)=" + sym.latex(Ep)  ) )

\(\displaystyle E(t)=\frac{\left(\sin{\left(t \right)} - \cos{\left(t \right)}\right)^{2}}{2} + \frac{\left(\sin{\left(t \right)} + \cos{\left(t \right)}\right)^{2}}{2}\)

\(\displaystyle E'(t)=0\)

3.4 ✨ Exercice (scipy - 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]\).
1. Afficher dans un autre répère \(x\mapsto y\). 1. On définit la fonction
\( t\mapsto E(t)=\frac{x^2(t)}{2}+\frac{y^2(t)}{2}. \)
Calculer analytiquement \(E'(t)\). Afficher ensuite \(E(t)\) calculé à partir des solutions approchées.

%reset -f
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})
from scipy.integrate import odeint

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');

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

E = lambda yy : (yy[0]**2+yy[1]**2)/2

E0 = E(yy0)
EE = E([solx,soly])

print(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}") # norme infinie de l'erreur absolue

plt.figure(figsize=(17,7))
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])
plt.plot(tt,EE)
plt.grid();

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

3.5 ✨ Exercice (scipy - système - 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, considerer 4 données initiales proches de ce point et calculer la solution particulière approchée associée à tel 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({<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é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 l’équation \(x^4-25x^2+12^2=0\):

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