# -*- coding: utf-8 -*-
"""
Created on Tue Jul 19 11:42:04 2022

@author: starons
"""

"""
Ce programme résoud numériquement l'équation différentielle de l'oscillateur 
harmonique non amorti en adaptant la méthode d'Euler explicite et compare à
la solution analytique. 
"""

# On importe les bibliothèques utiles
# -----------------------------------
import numpy as np
from matplotlib import pyplot as plt

# Valeurs numériques caractéristiques du problème
# -----------------------------------------------
omega0 = 2                              # pulsation propre
x0 = 4                                  # conditions initiales 
xp0 = 0

# Paramètres de la résolution numérique
# ------------------------------------- 
t_fin = 15/omega0                       # durée
N = 100                                 # nombre de points
h = t_fin/(N-1)                         # pas de temps

# Implémentation de la méthode d'Euler explicite
# ---------------------------------------------- 
t = np.linspace(0,t_fin,N)              # découpage du temps 
x = np.zeros(N)                         # création de tableaux vides à N éléments (remplis de zéros)
xp = np.zeros(N)
x[0] = x0                               # initialisation 
xp[0] = xp0
 
for n in range(N-1):                    # remplissage par récurrence
    x[n+1] = x[n]+xp[n]*h               # schéma numérique 
    xp[n+1] = xp[n]-omega0**2*x[n]*h    

# On connait une solution analytique
# ---------------------------------- 
t_analytique = np.linspace(0,t_fin,1000)# pour disposer d'un tableau plus dense de valeurs de temps   
x_analytique = x0*np.cos(omega0*t_analytique)+xp0/omega0*np.sin(omega0*t_analytique)
xp_analytique = -x0*omega0*np.sin(omega0*t_analytique)+xp0*np.cos(omega0*t_analytique)

# Représentations graphiques
# --------------------------
plt.subplot(1,2,1)
plt.plot(t,x,color='black',marker='.',label="solution numérique")
plt.plot(t_analytique,x_analytique,color='red',linestyle='dashed',label="solution analytique")
plt.xlabel("temps"),plt.ylabel("position")
plt.title("Oscillateur harmonique non amorti")
plt.legend(loc="lower left"),plt.grid()

plt.subplot(1,2,2)
plt.plot(x,xp/omega0,color='black',marker='.',label="solution numérique")
plt.plot(x_analytique,xp_analytique/omega0,color='red',linestyle='dashed',label="solution analytique")
plt.axis("equal")                         # pour avoir un repère orthonormé
plt.xlabel("position"),plt.ylabel("vitesse")
plt.title("Trajectoire de phase")
plt.legend(),plt.grid()