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.