5. Rappel sur les EDO#
Marc BUFFAT, dpt mécanique, Université Claude Bernard Lyon 1
%matplotlib inline
import numpy as np
import sympy as sp
from sympy.plotting import plot
import matplotlib.pyplot as plt
from IPython.core.display import HTML
from IPython.display import display
5.1. Equation différentielle (EDO) d’ordre 1#
5.1.1. EDO d’ordre 1: forme générale explicite#
5.1.2. EDO d’ordre 1 linéaire homogène#
5.1.3. EDO d’ordre 1 linéaire non homogène#
solution = solution générale homogène + solution particulière
solution particulière : méthode variation de la constante
solution particulière évidente
from sympy import Derivative, Eq, dsolve
t = sp.symbols('t')
y = sp.symbols('y',cls=sp.Function)
eq = Eq(Derivative(y(t),t),y(t)/2+1)
CI = {y(0):2}
display("EDO",eq,"CI",CI)
sol = dsolve(eq,ics=CI)
display("solution ",sol)
plot(sol.rhs,(t,0,4),title="solution",ylabel="y(t)");
'EDO'
'CI'
{y(0): 2}
'solution '
5.1.4. Système d’équations différentielles d’ordre 1#
5.1.4.1. système EDO d’ordre 1: forme générale#
5.1.4.2. Exemple 1#
soit
\(\textbf{2 valeurs propres de A }\): \(\pm i\)
changement de variable
les coefficients sont les vecteurs propres de A:
[i,1] vecteur propre associé à \(-i\)
[-i,1] vecteur propre associé à \(i\)
équations découplées
solutions découplées
solutions générales
pour \(C_1=C_2=1/2\)
from sympy import Matrix
y1, y2 = sp.symbols("y_1 y_2", cls=sp.Function)
Y = Matrix([[y1(t)],[y2(t)]])
A = Matrix([[0,-1],[1,0]])
EQ = Eq(Derivative(Y,t).doit(),A*Y)
Y0 = Matrix([[2],[-1]])
display("systeme EDO",EQ, "CI",Y0)
VP = A.eigenvects()
display("VP:",VP)
l1 = VP[0][0]
l2 = VP[1][0]
D = Matrix([[l1,0],[0,l2]])
display("Matrice diagonale D v.p. :",D)
L1 = VP[0][2][0]
L2 = VP[1][2][0]
Q = Matrix([[L1[0],L2[0]],[L1[1],L2[1]]])
display("Matrice de passage Q et Q^-1=",Q,Q.inv())
display("A = Q D Q^-1",Q * D * Q.inv())
z1, z2 = sp.symbols("z_1 z_2", cls=sp.Function)
Z = Matrix([[z1(t)],[z2(t)]])
display("changement de variable ",Eq(Z,Q.inv()*Y))
EQ1 = Eq(Derivative(Z,t).doit(),D*Z)
Z0 = Q.inv()*Y0
display("EDO diagonalisé",EQ1,"C.I",Z0)
sol1 = dsolve(Derivative(z1(t),t)-l1*z1(t),ics={z1(0):Z0[0]})
sol2 = dsolve(Derivative(z2(t),t)-l2*z2(t),ics={z2(0):Z0[1]})
display("solution en z:",sol1,sol2)
Ys = Q*Matrix([[sol1.rhs],[sol2.rhs]])
display("solution en Y",Eq(Y,Ys))
'systeme EDO'
'CI'
'VP:'
[(-I,
1,
[Matrix([
[-I],
[ 1]])]),
(I,
1,
[Matrix([
[I],
[1]])])]
'Matrice diagonale D v.p. :'
'Matrice de passage Q et Q^-1='
'A = Q D Q^-1'
'changement de variable '
'EDO diagonalisé'
'C.I'
'solution en z:'
'solution en Y'
5.1.4.3. EDO d’ordre 1 linéaire homogène#
si \(\mathbf{A}\) possède n valeurs propres réelles \(\lambda_k\), la solution générale s’écrit
les \(\alpha_k\) sont solution du système linéaire
On a donc diagonalisée la matrice \(\mathbf{A} = Q*\mathbf{D}*Q^{-1}\) pour obtenir un système EDO linéaire découplée $\( \dot{Z} = \mathbf{D} Z \mbox{ avec } Z = \mathbf{Q}^{-1} {Y} \)$
5.1.4.4. Exemple 2#
from sympy import Matrix
y1, y2 = sp.symbols("y_1 y_2", cls=sp.Function)
Y = Matrix([[y1(t)],[y2(t)]])
A = Matrix([[2,-1],[-1,2]])
EQ = Eq(Derivative(Y,t).doit(),A*Y)
Y0 = Matrix([[2],[-1]])
display("systeme EDO",EQ, "CI",Y0)
VP = A.eigenvects()
l1 = VP[0][0]
l2 = VP[1][0]
D = Matrix([[l1,0],[0,l2]])
display("Matrice diagonale D v.p. :",D)
L1 = VP[0][2][0]
L2 = VP[1][2][0]
Q = Matrix([[L1[0],L2[0]],[L1[1],L2[1]]])
display("Matrice de passage Q=",Q)
display("Q D Q^-1",Q * D * Q.inv())
z1, z2 = sp.symbols("z_1 z_2", cls=sp.Function)
Z = Matrix([[z1(t)],[z2(t)]])
display("changement de variable ",Eq(Z,Q.inv()*Y))
EQ1 = Eq(Derivative(Z,t).doit(),D*Z)
Z0 = Q.inv()*Y0
display("EDO diagonalisé",EQ1,"C.I",Z0)
sol1 = dsolve(Derivative(z1(t),t)-l1*z1(t),ics={z1(0):Z0[0]})
sol2 = dsolve(Derivative(z2(t),t)-l2*z2(t),ics={z2(0):Z0[1]})
display("solution en z:",sol1,sol2)
Ys = Q*Matrix([[sol1.rhs],[sol2.rhs]])
display("solution en Y",Eq(Y,Ys))
'systeme EDO'
'CI'
'Matrice diagonale D v.p. :'
'Matrice de passage Q='
'Q D Q^-1'
'changement de variable '
'EDO diagonalisé'
'C.I'
'solution en z:'
'solution en Y'
5.2. Equation différentielle (EDO) d’ordre 2#
5.2.1. EDO d’ordre 2: forme générale explicite#
5.2.2. EDO d’ordre 2 linéaire homogène#
Cas général:
racines du polynôme caractéristique \(a\lambda^2 + b\lambda + c\)
Equation en mécanique: système dynamique avec force de rappel (fonction de y), et un amortissement (fonction de \( \dot{y}\)
solution fonction du signe de \(\Delta = \lambda_0^2 -\omega_0^2\)
\(\Delta>0\) régime apériodique
\(\Delta=0\) régime critique
\(\Delta<0\) régime pseudo-péroidique
lambda0, omega0 = sp.symbols("lambda_0 omega_0")
eq = Derivative(y(t),t,t) + 2*lambda0*Derivative(y(t),t) + omega0**2*y(t)
display("EDO ",eq)
# amortie
Valnum={lambda0:2,omega0:1}
# critique
#Valnum={lambda0:1,omega0:1}
# pseudo per.
#Valnum={lambda0:1/sp.S(4),omega0:1}
# periodique
Valnum={lambda0:0,omega0:1}
Delta = float((lambda0*2-omega0**2).subs(Valnum))
display(Valnum,"Delta=",Delta)
CI={y(0):1,Derivative(y(t),t).subs(t,0):0}
display("CI:",CI)
sol = dsolve(eq.subs(Valnum),ics=CI)
display("solution :",sol)
plot(sol.rhs,(t,0,20),title="solution EDO pour $\Delta ={}$".format(Delta),ylabel="y(t)");
'EDO '
{lambda_0: 0, omega_0: 1}
'Delta='
-1.0
'CI:'
{y(0): 1, Subs(Derivative(y(t), t), t, 0): 0}
'solution :'
5.2.3. EDO d’ordre 2 linéaire avec forçage \(f(t)\)#
solution = solution équation homogène + solution particulière
si forçage harmonique \(f(t)= A \cos(\omega_1 t)\)
la solution dépend de la valeur de \(\omega_1\) par rapport à \(\omega_0\)
si \(\omega_1 \neq \omega_0\) solution particulière \( y(t) = \alpha \cos{\omega_1 t}\)
si \(\omega_1 = \omega_0\) résonnance: solution particulière \( y(t) = t \sin{\omega_0 t}\)
omega1, A = sp.symbols("omega_1 A")
eq = Eq(Derivative(y(t),t,t) + omega0**2*y(t),A*sp.cos(omega1*t))
display("EDO ",eq)
# sans forcage
Valnum={omega0:1,omega1:1/sp.S(2),A:0/sp.S(3)}
# avec forcage
#Valnum={omega0:1,omega1:1/sp.S(2),A:1/sp.S(3)}
# avec forcage +
#Valnum={omega0:1,omega1:6/sp.S(5),A:1/sp.S(3)}
# resonnance
#Valnum={omega0:1,omega1:1,A:1/sp.S(3)}
#
sol = dsolve(eq.subs(Valnum),ics=CI)
display("solution: ",sol)
plot(sol.rhs,(t,0,40),title="solution EDO pour ${}$".format(Valnum),ylabel="y(t)");
'EDO '
'solution: '
5.2.4. Transformation en EDO d’ordre 1#
Tous système d’équations différentielles peut se mettre sous la forme générique d’un système différentielle d’ordre 1
La transformation utilise un changement de variable
pour une EDO d’ordre n, on prend comme inconnues la fonction et toutes les dérivées jusqu’à l’ordre n-1
Important cette forme d’EDO d’ordre 1 est la forme générique pour la résolution numérique des EDO
5.2.4.1. transformation EDO d’ordre 2#
équation générique d’ordre 2
système d’ordre 1 équivalent
5.2.4.2. exemple#
équation d’un système dynamique sans amortissement:
transformation en EDO d’ordre 1
valeurs propres:
équation caractéristique
\[\det(\lambda Id -A) = \lambda^2 + \omega_0^2 = 0\]2 racines complexes conjuguées \(\lambda = \pm i \omega_0 \)
d’où la solution générale
\[ y(t) = Z_1 e^{i \omega_0 t} + Z_2 e^{- i \omega_0 t} \]qui correspond à la solution périodique
\[ y(t) = A \cos{\omega_0 t} + B \sin{\omega_0 t} \]
5.3. Résolution numérique des EDO#
transformation en un système d’EDO d’ordre 1
résolution d’EDO d’ordre 1 :
5.3.1. méthode d’Euler#
recherche de la solution \(y=F(t)\) de l’EDO \(\frac{dy}{dt} = f(y(t),t) \)$
choix pas de discrétisation \(h=dt\)
calcul d’une solution approchée aux instants \(t= 0, h, 2h, .. kh\)
notation \(t_k = k h\), solution approchée \(y_k \approx F(t_k)\)
calcul itératif de \(y_{k+1}\) en fonction de \(y_k\) par discrétisation de l’EDO
5.3.2. algorithme d’Euler#
données: fonction second membre
f(y,t)parametres temps total
Td’intégration et nbre de pas d’intégrationnh= T/n calcul du pas d’intégrationY= tableau dimension n+1 solution approchée:Y[0] = y0 condition initiale t = 0 pour k de 1 a n Y[k] = Y[k-1] + h*f(Y[k_1],t) t = t + h fin boucle
5.3.3. méthode de Runge Kutta 2#
pour améliorer la précision de la solution approchée (ordre 1 avec Euler), on va utiliser une discrétisation plus précise de la dérivée entre \(t_k\) et \(t_k+h\), en calculant la valeur au milieu de l’intervalle à \(t_{k+1/2} = t_k + h/2\) à l’aide de la formule de RungeKutta 2 suivante:
dans laquelle la valeur inconnue \(y_{k+1/2}\) est calculée avec la méthode d’Euler
5.3.4. algorithme de Runge Kutta 2#
la boucle d’itération de Runge Kutta 2 s’écrit
Y[0] = y0 condition initiale
t = 0
pour k de 1 a n
ym = Y[k-1] + h/2*f(Y[k_1],t)
Y[k]= Y[k-1] + h*f(ym,t+h/2)
t = t + h
fin boucle