Code
import sympy as sp
from IPython.display import display, Markdown # For pretty printing
import matplotlib.pyplot as plt
import numpy as np
# plt.rcdefaults()
# plt.rcParams.update({'font.size': 16})
Gloria Faccanoni
06 mars 2026
En L1 vous avez appris comment calculer la solution générale de quelques équations différentielles ordinaires (EDO). Nous allons maintenant voir comment automatiser ces calculs avec le module sympy conçu pour du calcul symbolique/formel/exact.
Cf. https://docs.sympy.org/latest/modules/solvers/ode.html
Bien lire et comprendre les deux exemples pour résoudre des EDO du premier ordre et des systèmes d’EDO du premier ordre, puis résoudre les deux exercices proposés.
import sympy as sp
from IPython.display import display, Markdown # For pretty printing
import matplotlib.pyplot as plt
import numpy as np
# plt.rcdefaults()
# plt.rcParams.update({'font.size': 16})
Considérons le problème de Cauchy
\(
\begin{cases}
u'(x)=-3x^2u(x)+6x^2,\\
u(0)=4.
\end{cases}
\)
On peut résoudre ce problème de Cauchy de manière analytique en utilisant le module sympy de la manière suivante.
Définition du problème de Cauchy :
# Variables
x = sp.Symbol('x')
u = sp.Function('u')
# Equation différentielle
left = sp.diff(u(x),x) # la dérivée de u par rapport à x
right = 6*x**2-3*x**2*u(x) # l'expression de l'équation différentielle
edo = sp.Eq( left , right ) # l'équation différentielle
display(edo)
# Conditions initiales
x0 = 0
u0 = 4
display(sp.Eq(u(x0),u0))
\(\displaystyle \frac{d}{d x} u{\left(x \right)} = - 3 x^{2} u{\left(x \right)} + 6 x^{2}\)
\(\displaystyle u{\left(0 \right)} = 4\)
Calcul de la solution exacte (2 méthodes : une méthode directe et une méthode pas à pas) :
# Méthode 1 : méthode directe
# ============================
display(Markdown("- Méthode 1 : méthode directe"))
solpar = sp.dsolve(edo,u(x), ics={u(x0):u0}) # l'edo, la fonction u(x) inconnue, les conditions initiales comme dictionnaire
display(solpar)
\(\displaystyle u{\left(x \right)} = 2 + 2 e^{- x^{3}}\)
# Méthode 2 : méthode pas-à-pas
# ==============================
# On va calculer la solution générale, puis la solution particulière en prenant en compte
# la condition initiale
display(Markdown("- Méthode 2 : méthode pas-à-pas"))
# Étape 1 : calcul de la solution générale (qui contient la constante d'intégration) :
solgen = sp.dsolve(edo,u(x))
display(solgen)
# Étape 2 : calcul de la constante d'intégration en prenant en compte la condition initiale :
cond_init = sp.Eq( u0, solgen.rhs.subs(x,x0) )
display(cond_init)
consts = sp.solve( cond_init , dict=True)[0]
display(consts)
# Étape 3 : calcul de la solution particulière :
solpar = solgen.subs(consts)
display(solpar)
\(\displaystyle u{\left(x \right)} = C_{1} e^{- x^{3}} + 2\)
\(\displaystyle 4 = C_{1} + 2\)
{C1: 2}
\(\displaystyle u{\left(x \right)} = 2 + 2 e^{- x^{3}}\)
Conversion en une fonction de numpy et affichage :
# Transformation en fonction numpy
func = sp.lambdify(x,solpar.rhs,'numpy')
# Tracé de la solution
plt.figure(figsize=(8,5))
xx = np.linspace(0,3,101)
yy = func(xx)
plt.plot(xx,yy)
plt.xlabel('x')
plt.ylabel('u')
plt.grid();

\( \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)\).
# Variables
t = sp.Symbol('t')
# Fonctions inconnues
y1 = sp.Function(r'y_1')
y2 = sp.Function(r'y_2')
# Termes de droite
right1 = y1(t)-y2(t)
right2 = y2(t)-y1(t)
# Équations différentielles
edo1 = sp.Eq( sp.diff(y1(t),t) , right1 )
edo2 = sp.Eq( sp.diff(y2(t),t) , right2 )
display(edo1)
display(edo2)
# Conditions initiales
t0 = 0
y1_0 = 2
y2_0 = 1
display(sp.Eq(y1(t0),y1_0))
display(sp.Eq(y2(t0),y2_0))
\(\displaystyle \frac{d}{d t} y_{1}{\left(t \right)} = y_{1}{\left(t \right)} - y_{2}{\left(t \right)}\)
\(\displaystyle \frac{d}{d t} y_{2}{\left(t \right)} = - y_{1}{\left(t \right)} + y_{2}{\left(t \right)}\)
\(\displaystyle y_{1}{\left(0 \right)} = 2\)
\(\displaystyle y_{2}{\left(0 \right)} = 1\)
# Méthode 1 : méthode directe
# ============================
solpar_1, solpar_2 = sp.dsolve( [ edo1, edo2 ] ,
ics = { y1(t0):y1_0, y2(t0):y2_0 } )
display(solpar_1)
display(solpar_2)
# Méthode 2 : méthode pas-à-pas
# ==============================
# solgen = sp.dsolve([edo1,edo2],[y1(t),y2(t)])
# display(solgen)
# consts = sp.solve( [sp.Eq( y1_0, solgen[0].rhs.subs(t,t0)) ,
# sp.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)
\(\displaystyle y_{1}{\left(t \right)} = \frac{e^{2 t}}{2} + \frac{3}{2}\)
\(\displaystyle y_{2}{\left(t \right)} = \frac{3}{2} - \frac{e^{2 t}}{2}\)
# Affichage
# =========
# Transformation des solutions en fonctions numpy
func_1 = sp.lambdify(t,solpar_1.rhs,'numpy')
func_2 = sp.lambdify(t,solpar_2.rhs,'numpy')
plt.figure(figsize=(16,5))
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.title("Solutions")
plt.subplot(1,2,2)
plt.plot(yy_1,yy_2)
plt.xlabel(r'$y_1$')
plt.ylabel(r'$y_2$')
plt.title("Portrait de phase");

Calculons l’unique solution du problème de Cauchy
\( \begin{cases} y''(t) + 2y'(t) + 2y(t) = 2\sin(t), &t>0\\ y(0)=1,\\ y'(0)=0. \end{cases} \)
Méthode directe sur l’équation différentielle d’ordre 2 :
# Variables
t = sp.Symbol('t')
# Fonction inconnue
y = sp.Function('y')
# Équation différentielle
left = sp.diff(y(t),t,2) + 2*sp.diff(y(t),t) + 2*y(t)
right = 2*sp.sin(t)
edo = sp.Eq( left , right ) # l'équation différentielle
display(edo)
# Conditions initiales
t0 = 0
y0 = 1
yp0 = 0
display(sp.Eq(y(t0),y0))
display(sp.Eq(sp.diff(y(t),t).subs(t,t0),yp0))
# Méthode 1 : méthode directe
# ============================
solpar = sp.dsolve(edo,y(t), ics={y(t0):y0, sp.diff(y(t),t).subs(t,t0):yp0})
display(solpar)
\(\displaystyle 2 y{\left(t \right)} + 2 \frac{d}{d t} y{\left(t \right)} + \frac{d^{2}}{d t^{2}} y{\left(t \right)} = 2 \sin{\left(t \right)}\)
\(\displaystyle y{\left(0 \right)} = 1\)
\(\displaystyle \left. \frac{d}{d t} y{\left(t \right)} \right|_{\substack{ t=0 }} = 0\)
\(\displaystyle y{\left(t \right)} = \left(\frac{7 \sin{\left(t \right)}}{5} + \frac{9 \cos{\left(t \right)}}{5}\right) e^{- t} + \frac{2 \sin{\left(t \right)}}{5} - \frac{4 \cos{\left(t \right)}}{5}\)
Vérifions si elle satisfait bien notre problème de Cauchy :
# Qua donne la solution trouvée lorsqu'on la remplace dans l'équation différentielle ?
display( (sp.diff(solpar.rhs,t,2) + 2*sp.diff(solpar.rhs,t) + 2*solpar.rhs -2*sp.sin(t) ).simplify() )
# Que vaut la solution trouvée lorsqu'on la remplace dans les conditions initiales ?
display( solpar.rhs.subs(t,t0) )
display( sp.diff(solpar.rhs,t).subs(t,t0) )
\(\displaystyle 0\)
\(\displaystyle 1\)
\(\displaystyle 0\)
# Affichage de la solution
func = sp.lambdify(t,solpar.rhs,'numpy')
tt = np.linspace(0,10,1001)
yy = func(tt)
plt.figure(figsize=(16,5))
# Solution
# ========
plt.subplot(1,2,1)
plt.plot(tt,yy)
plt.title(r"$t\mapsto y(t)$")
plt.grid();
# Portraits de phase
# ==================
plt.subplot(1,2,2)
plt.plot(yy,np.gradient(yy,tt))
plt.title(r"$y\mapsto y'(y)$")
plt.grid();

Méthode pas à pas sur l’équation différentielle d’ordre 2 :
# Bonus : résolution de l'équation différentielle avec la méthode pas-à-pas
# ========================================================================
# Étape 1 : calcul de la solution générale (qui contient la constante d'intégration) :
solgen = sp.dsolve(edo,y(t))
display(solgen)
# Étape 2 : calcul de la constante d'intégration en prenant en compte les conditions initiales :
cond_init1 = sp.Eq( y0, solgen.rhs.subs(t,t0) )
cond_init2 = sp.Eq( yp0, sp.diff(solgen.rhs,t).subs(t,t0) )
display(cond_init1)
display(cond_init2)
consts = sp.solve( [cond_init1,cond_init2] , dict=True)[0]
display(consts)
# Étape 3 : calcul de la solution particulière :
solpar = solgen.subs(consts)
display(solpar)
\(\displaystyle y{\left(t \right)} = \left(C_{1} \sin{\left(t \right)} + C_{2} \cos{\left(t \right)}\right) e^{- t} + \frac{2 \sin{\left(t \right)}}{5} - \frac{4 \cos{\left(t \right)}}{5}\)
\(\displaystyle 1 = C_{2} - \frac{4}{5}\)
\(\displaystyle 0 = C_{1} - C_{2} + \frac{2}{5}\)
{C1: 7/5, C2: 9/5}
\(\displaystyle y{\left(t \right)} = \left(\frac{7 \sin{\left(t \right)}}{5} + \frac{9 \cos{\left(t \right)}}{5}\right) e^{- t} + \frac{2 \sin{\left(t \right)}}{5} - \frac{4 \cos{\left(t \right)}}{5}\)
Méthode directe ou pas à pas après transformation en un système d’équations différentielles d’ordre 1 :
# Bonus : transformation en un système d'équations différentielles du premier ordre
# ================================================================================
# Variables
t = sp.Symbol('t')
y = sp.Function('y')
z = sp.Function('z')
# Équations
edo1 = sp.Eq( sp.diff(y(t),t) , z(t) )
edo2 = sp.Eq( sp.diff(z(t),t) , 2*sp.sin(t) - 2*z(t) - 2*y(t) )
display(edo1)
display(edo2)
# Conditions initiales
t0 = 0
y0 = 1
yp0 = 0
display(sp.Eq(y(t0),y0))
display(sp.Eq(z(t0),yp0))
# Méthode 1 : méthode directe
# ============================
display(Markdown("- Méthode 1 : méthode directe"))
solpar = sp.dsolve( [edo1,edo2], [y(t),z(t)], ics={y(t0):y0, z(t0):yp0} )
display(solpar)
# Méthode 2 : méthode pas-à-pas
# ==============================
display(Markdown("- Méthode 2 : méthode pas-à-pas"))
# On va calculer la solution générale, puis la solution particulière en prenant en compte les conditions initiales
solgen = sp.dsolve( [edo1,edo2], [y(t),z(t)] )
display(solgen)
# Étape 2 : calcul de la constante d'intégration en prenant en compte les conditions initiales :
cond_init1 = sp.Eq( y0, solgen[0].rhs.subs(t,t0) )
cond_init2 = sp.Eq( yp0, solgen[1].rhs.subs(t,t0) )
display(cond_init1)
display(cond_init2)
consts = sp.solve( [cond_init1,cond_init2] , dict=True)[0]
display(consts)
# Étape 3 : calcul de la solution particulière :
solpar = [sol.subs(consts) for sol in solgen]
display(solpar)
\(\displaystyle \frac{d}{d t} y{\left(t \right)} = z{\left(t \right)}\)
\(\displaystyle \frac{d}{d t} z{\left(t \right)} = - 2 y{\left(t \right)} - 2 z{\left(t \right)} + 2 \sin{\left(t \right)}\)
\(\displaystyle y{\left(0 \right)} = 1\)
\(\displaystyle z{\left(0 \right)} = 0\)
[Eq(y(t), 2*sin(t)**3/5 - 4*sin(t)**2*cos(t)/5 + 2*sin(t)*cos(t)**2/5 - 4*cos(t)**3/5 + 7*exp(-t)*sin(t)/5 + 9*exp(-t)*cos(t)/5),
Eq(z(t), 4*sin(t)**3/5 + 2*sin(t)**2*cos(t)/5 + 4*sin(t)*cos(t)**2/5 + 2*cos(t)**3/5 - 16*exp(-t)*sin(t)/5 - 2*exp(-t)*cos(t)/5)]
[Eq(y(t), -(C1/2 - C2/2)*exp(-t)*cos(t) + (C1/2 + C2/2)*exp(-t)*sin(t) + 2*sin(t)**3/5 - 4*sin(t)**2*cos(t)/5 + 2*sin(t)*cos(t)**2/5 - 4*cos(t)**3/5),
Eq(z(t), C1*exp(-t)*cos(t) - C2*exp(-t)*sin(t) + 4*sin(t)**3/5 + 2*sin(t)**2*cos(t)/5 + 4*sin(t)*cos(t)**2/5 + 2*cos(t)**3/5)]
\(\displaystyle 1 = - \frac{C_{1}}{2} + \frac{C_{2}}{2} - \frac{4}{5}\)
\(\displaystyle 0 = C_{1} + \frac{2}{5}\)
{C1: -2/5, C2: 16/5}
[Eq(y(t), 2*sin(t)**3/5 - 4*sin(t)**2*cos(t)/5 + 2*sin(t)*cos(t)**2/5 - 4*cos(t)**3/5 + 7*exp(-t)*sin(t)/5 + 9*exp(-t)*cos(t)/5),
Eq(z(t), 4*sin(t)**3/5 + 2*sin(t)**2*cos(t)/5 + 4*sin(t)*cos(t)**2/5 + 2*cos(t)**3/5 - 16*exp(-t)*sin(t)/5 - 2*exp(-t)*cos(t)/5)]
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}\)
Remarque : ne pas écrire 3/2 mais sympy.Rational(3,2) pour éviter que ce terme soit évalué comme un flottant.
t = sp.Symbol('t')
y = sp.Function('y')
right = sp.Rational(3,2)*y(t)*(1-y(t)/6) # Attention à ne pas écrire 3/2 qui donne 1.5
edo = sp.Eq( sp.diff(y(t),t) , right )
display(edo)
t0 = 0
y0 = 1
display(sp.Eq(y(t0),y0))
\(\displaystyle \frac{d}{d t} y{\left(t \right)} = \frac{3 \cdot \left(1 - \frac{y{\left(t \right)}}{6}\right) y{\left(t \right)}}{2}\)
\(\displaystyle y{\left(0 \right)} = 1\)
# Méthode directe
# ===============
solpar = sp.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}}\)
# Affichage
# =========
# Étape 1 : conversion de la solution en fonction numpy
func = sp.lambdify(t,solpar.rhs,'numpy')
# Étape 2 : plot
plt.figure(figsize=(8,5))
tt = np.linspace(0,10,101)
yy = func(tt)
plt.plot(tt,yy)
plt.grid();

# Calcul de la limite
# ===================
limite = sp.limit(solpar.rhs,t,sp.oo)
display( sp.Eq( sp.Limit(y(t),t,sp.oo) , limite ) )
\(\displaystyle \lim_{t \to \infty} y{\left(t \right)} = 6\)
Bonus : dans cet exemple, si on veut d’abord trouver la solution générale, on s’aperçoit qu’il y a deux expressions possibles. On doit alors choisir celle qui convient à notre condition initiale. Comme la solution du problème de Cauchy est unique, une seule des deux expressions est correcte.
# Méthode pas-à-pas
# =================
# solgenLIST = sp.dsolve(edo,y(t))
# display(solgenLIST)
# On teste avec la première expression, si ça ne marche pas, on essaie avec la seconde
# =====================================================================================
# solgen = solgenLIST[0]
# display(solgen)
# equation = sp.Eq( y0, solgen.rhs.subs(t,t0))
# display(equation)
# consts = sp.solve( equation , dict=True)[0]
# display(consts)
# solpar = solgen.subs(consts)
# display(solpar)
Calculer la solution 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} \)
pour \(x(0)=-1\) et \(y(0)=1\).
Afficher \(t\mapsto x\) et \(t\mapsto y\) dans un même graphe.
Afficher ensuite \(x\mapsto y\).
Calculer analytiquement \(E'(t)\) où \(t\mapsto E(t)\) est l’énergie mécanique du système définie par \( E(t)=\frac{x^2(t)}{2}+\frac{y^2(t)}{2}. \)
Bonus : transformer le système en une équation différentielle d’ordre 2 et résoudre le problème de Cauchy associé.
La première équation s’écrit
\( x'(t) = y(t) \)
et la deuxième
\( y'(t) = -x(t). \)
t = sp.Symbol('t')
x = sp.Function('x')
y = sp.Function('y')
edo1 = sp.Eq( sp.diff(x(t),t) , y(t) )
edo2 = sp.Eq( sp.diff(y(t),t) , -x(t) )
display(edo1)
display(edo2)
t_0 = 0
x_0 = -1
y_0 = 1
display(sp.Eq(x(t_0),x_0))
display(sp.Eq(y(t_0),y_0))
\(\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)}\)
\(\displaystyle x{\left(0 \right)} = -1\)
\(\displaystyle y{\left(0 \right)} = 1\)
# Méthode directe
# ===============
solpar_1, solpar_2 = sp.dsolve( [edo1,edo2],
[x(t),y(t)],
ics={x(t_0):x_0, y(t_0):y_0} )
display(solpar_1)
display(solpar_2)
# Méthode pas-à-pas
# =================
# solgen = sp.dsolve([edo1,edo2],[x(t),y(t)])
# display(solgen)
# consts = sp.solve( [ sp.Eq( x_0, solgen[0].rhs.subs(t,t_0)) , sp.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)
\(\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)}\)
# Conversion en fonctions numpy
# =============================
func_1 = sp.lambdify(t,solpar_1.rhs,'numpy')
func_2 = sp.lambdify(t,solpar_2.rhs,'numpy')
# Affichage
# =========
plt.figure(figsize=(16,5))
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.title('Solutions')
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()
plt.title('Portrait de phase');

La fonction \(t\mapsto E(t) = \frac{x^2(t)}{2} + \frac{y^2(t)}{2}\) vérifie
\( \begin{aligned} E'(t) &= x(t)x'(t)+y(t)y'(t) \\ &= x(t)\varphi_1(t,x(t),y(t))+y(t)(t)\varphi_2(t,x(t),y(t)) \\ &= x(t) y(t) + y(t)(-x(t)) \\ &= 0. \end{aligned} \)
Vérifions qu’effectivement \(E\) est une constante en remplaçant \(x\) et \(y\) par leurs expressions analytiques.
# On calcule E et E' analytiquement
# ================================
# La fonction E(t) = x(t)^2/2 + y(t)^2/2
E = (solpar_1.rhs)**2/2+(solpar_2.rhs)**2/2
E = E.simplify()
display(sp.Eq( sp.Symbol('E(t)') , E ))
# La dérivée de E par rapport à t
Ep = sp.diff(E,t).simplify()
display(sp.Eq( sp.Symbol("E'(t)") , Ep ))
\(\displaystyle E(t) = 1\)
\(\displaystyle E'(t) = 0\)
Transformation et résolution du système en une seule équation différentielle du deuxième ordre
On a \(x'(t)=y(t)\) et \(y'(t)=-x(t)\), on obtient alors \(x''(t)=y'(t)=-x(t)\).
De plus, \(x(0)=-1\) et \(x'(0)=y(0)=1\).
# Variables
t = sp.Symbol('t')
# Fonction inconnue
x = sp.Function('x')
# Équation différentielle
left = sp.diff(x(t),t,2)
right = -x(t)
edo = sp.Eq( left , right ) # l'équation différentielle
display(edo)
# Conditions initiales
t0 = 0
x0 = -1
xp0 = 1
display(sp.Eq(x(t0),x0))
display(sp.Eq(sp.diff(x(t),t).subs(t,t0),xp0))
# Méthode 1 : méthode directe
# ============================
solpar_x = sp.dsolve(edo,x(t), ics={x(t0):x0, sp.diff(x(t),t).subs(t,t0):xp0})
display(solpar_x)
solpar_y = sp.diff(solpar_x.rhs,t)
display( sp.Eq( y(t) , solpar_y ) )
\(\displaystyle \frac{d^{2}}{d t^{2}} x{\left(t \right)} = - x{\left(t \right)}\)
\(\displaystyle x{\left(0 \right)} = -1\)
\(\displaystyle \left. \frac{d}{d t} x{\left(t \right)} \right|_{\substack{ t=0 }} = 1\)
\(\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)}\)