Code
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
odeint
du module SciPy
odeint
solve_ivp
?sympy
Gloria FACCANONI
22 mars 2024
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
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).
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))\).
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.
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\)
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
.
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).
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
.
quiver
du module matplt.plotlib
permet de tracer un champ de vecteurs: quiver(x_pos, y_pos, x_dir, y_dir, color)
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');
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})) \)
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
Le tracé des évolutions de \(y_2\) en fonction de \(y_1\) peut être obtenu par
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
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\).
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
sympy
https://docs.sympy.org/latest/modules/solvers/ode.html
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 :
\(\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) :
\(\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 :
\(\displaystyle 4 = C_{1} + 2\)
\(\displaystyle \left\{ C_{1} : 2\right\}\)
Étape 3 : calcul de la solution particulière :
\(\displaystyle u{\left(x \right)} = 2 + 2 e^{- x^{3}}\)
Calcul de la solution exacte, Méthode directe :
\(\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 :
\(
\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)\).
\(\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({'font.size': 16})
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$');
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)\)?
\(\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}}\)
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({'font.size': 16})
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={y0}$');
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');
\(\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({'font.size': 16})
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();
\(\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\)
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({'font.size': 16})
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"{min(EE-E0) = :e}, {max(EE-E0) = :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
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}
\)
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');
[-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)