Etude de l'équation différentielle du pendule non linéaire:¶
$y'' = -\sin{y}$ avec la condition initiale $y(0) = \frac{17\pi}{18},\quad y'(0) = 0.$
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
Question 1. Résolution avec solveur intégré de Python¶
On réecrit l'équation sous forme d'un système du premier ordre en posant $y[0] = y$ et $y[1] = y'$
#fonction definissant l'équation différentielle
def pend(t, y):
dydt = [y[1], - np.sin(y[0])]
return dydt
# pendule initial quasiment vertical (170°) à vitesse nulle
y0 = [np.pi*17/18, 0.0]
#duree de simulation
T=64
# resolution avec schéma RK45
sol = solve_ivp(pend, [0, T], y0, dense_output=True)
# on veut suffisamment de points pour le graphe
t = np.linspace(0, T, 100)
y = sol.sol(t)
#tracé des courbes solutions
plt.plot(t, y[0], 'b', label='y')
plt.plot(t, y[1], 'g', label=" y' ")
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()
#plan de phases
plt.plot(y[0], y[1], 'b')
plt.axis('equal')
plt.xlabel('y')
plt.ylabel("y'")
plt.title('schéma RK45 de Python')
plt.grid()
plt.show()
Question 2. Schéma implicite-explicite¶
# duree de simulation
T = 64
# pas
h = 1e-2
# nombre de points de la subdivision
N = int(T/h)
#initialisation des variables
q = np.zeros(N)
p = np.zeros(N)
# conditions initiales
q[0], p[0] = y0
# schéma Euler symplectique
# p_{n+1} = p_n - h sin q_n
# q_{n+1} = q_n + h p_{n+1}
for k in range(N-1):
p[k+1] = p[k] - h*np.sin(q[k])
q[k+1] = q[k] + h*p[k+1]
#trace des courbes solutions
t=np.linspace(0,T,N)
plt.plot(t, q, 'b', label=r'$y$')
plt.plot(t, p, 'g', label=r"$y'$")
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()
# plan de phases
plt.plot(q, p, '-')
plt.axis('equal')
plt.xlabel('y')
plt.ylabel("y'")
plt.title('schéma Euler implicite-explicite')
plt.grid()
plt.show()
Question 3. Conservation de l'énergie.¶
#test de conservation de l'énergie
t=np.linspace(0,T,N)
y = sol.sol(t)
E= -np.cos(q[0])
print('Pour le schéma implicite-explicite:')
print('max |E(t) - E(0)|=',np.max(np.abs(0.5*p**2 - np.cos(q) - E)))
print('qmplitude de variation =',np.max(0.5*p**2 - np.cos(q)) - np.min(0.5*p**2 - np.cos(q)))
print('Pour le schéma RK45:')
print('max |E(t) - E(0)|=',np.max(np.abs(0.5*y[1]**2 - np.cos(y[0])-E)))
print('amplitude de variation=',np.max(0.5*y[1]**2 - np.cos(y[0]))-\
np.min(0.5*y[1]**2 - np.cos(y[0])))
Pour le schéma implicite-explicite: max |E(t) - E(0)|= 0.0076617528911484545 qmplitude de variation = 0.015308219870332818 Pour le schéma RK45: max |E(t) - E(0)|= 0.04505415112335387 amplitude de variation= 0.05104294341917559
t=np.linspace(0,T,N)
y = sol.sol(t)
plt.plot(t,0.5*p**2 - np.cos(q))
plt.plot(t,0.5*y[1]**2 - np.cos(y[0]))
plt.plot(t,-np.cos(y0[0])*np.ones(len(t)))
plt.xlabel('t')
plt.ylabel('energie')
plt.legend(['E(t) implicite-explicite','E(t) RK45','E(0)'])
plt.show()
Le schéma implicite-explicite conserve mieux l'énergie que le schéma standard de Python. Pourtant le schéma standard de Python est d'ordre 5 alors que le schéma implicite-explicite est d'ordre 1. On voit aussi sur que la trajectoire dans le plan de phase se referme mieux avec le schéma implicite-explicite. Avec le schéma de Python, les oscillations du pendule semblent perdre peu à peu de l'amplitude. Le pendule semble perdre de l'énergie avec Python.
Question 4. Calcul des périodes numériques.¶
# Obtenir les temps où la vitesse angulaire change de signe (passage par un maximum)
t_max = t[np.where(np.diff(np.sign(p)) < 0)]
# Calculer la période comme la différence entre deux temps successifs de passage par un maximum
periode = np.mean(np.diff(t_max))
print('la periode numérique est',round(periode,2))
la periode numérique est 15.33
# detecte quand la vitesse y' change de signe positif
# en devenant ensuite negative
# cela correspond aux maximums de y
def retour(t, y): return y[1]
retour.terminal = False
retour.direction = -1
sol = solve_ivp(pend, [0, T], y0,
dense_output=True, events=retour)
sol.t_events
[array([ 0. , 14.433816 , 28.04291236, 41.20155337, 54.22983391])]
# Obtenir les temps où les événements se produisent (passage par un maximum)
t_max = sol.t_events[0]
# Calculer la période comme la différence entre deux temps successifs de passage par un maximum
periode = np.mean(np.diff(t_max))
print('la periode numérique est',round(periode,2))
la periode numérique est 13.56
On trouve sur wikipedia https://fr.wikipedia.org/wiki/Pendule_simple#Cas_pleinement_non_lin%C3%A9aire" que la période théorique du pendule est 2.44 x 2 $\pi$
print('la periode théorique est',np.round(2*np.pi*2.44,2))
la periode théorique est 15.33
Le schéma implicite-explicite donne un période correcte, ce n'est pas le cas de solve_ivp qui sous-estime la période.
Question 5. Changement de condition initiale.¶
# conditions initiales
q[0], p[0] = 19*np.pi/20, 0
# schéma Euler symplectique
# p_{n+1} = p_n - h sin q_n
# q_{n+1} = q_n + h p_{n+1}
for k in range(N-1):
p[k+1] = p[k] - h*np.sin(q[k])
q[k+1] = q[k] + h*p[k+1]
# plan de phases
plt.plot(q, p, '-')
plt.axis('equal')
plt.xlabel('y')
plt.ylabel("y'")
plt.title('schéma Euler implicite-explicite')
plt.grid()
plt.show()
#trace des courbes solutions
t=np.linspace(0,T,N)
plt.plot(t, q, 'b', label=r'$y$')
plt.plot(t, p, 'g', label=r"$y'$")
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()
# pendule initial quasiment vertical (171°) à vitesse nulle
y0 = [np.pi*19/20, 0.0]
#duree de simulation
T=64
# resolution avec schéma RK45
sol = solve_ivp(pend, [0, T], y0, dense_output=True)
# on veut suffisamment de points pour le graphe
t = np.linspace(0, T, 100)
y = sol.sol(t)
#tracé des courbes solutions
plt.plot(t, y[0], 'b', label='y')
plt.plot(t, y[1], 'g', label=" y' ")
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()
#plan de phases
plt.plot(y[0], y[1], 'b')
plt.axis('equal')
plt.xlabel('y')
plt.ylabel("y'")
plt.title('schéma RK45 de Python')
#plt.grid()
plt.show()
Avec solve_ivp, on obtient un résultat non physique, après 3 périodes, le pendule fait un tour complet au lieu de revenir! Cela semble dû à un défaut de conservation de l'énergie. Le schéma implicite-explicite se comporte bien.
t=np.linspace(0,T,N)
y = sol.sol(t)
plt.plot(t,0.5*p**2 - np.cos(q))
plt.plot(t,0.5*y[1]**2 - np.cos(y[0]))
plt.plot(t,-np.cos(y0[0])*np.ones(len(t)))
plt.xlabel('t')
plt.ylabel('energie')
plt.legend(['E(t) implicite-explicite','E(t) RK45','E(0)'])
plt.show()
On voit un excès d'energie non physique à t= 50 environ là où le pendule fait un tour complet.