# Capacité numérique n°6 : système d'ordre 2 non linéaire

## Script 1 : méthode odeint pour l'étude du pendule simple non amorti

import numpy as np
import matplotlib.pyplot as plt
from math import *
from scipy.integrate import odeint

# Données du problème
p_propre = 1   #  pulsation propre de 1 rad/s
theta0 = 0.25  # CI sur l'angle initial successivement 0.25 rd, 2.5 rd, 0 rd
omega0 = 0    # CI sur la vitesse angulaire initiale successivement 0 rad/s, 0 rad/s , 2.2 rad/s

# Données pour l'observation
T = 20 # durée de l'étude
n = 200

# Calculs
def pendule(X,t):
    theta = X[0] # première ligne de X(t)
    omega = X[1] # seconde ligne de X(t)
    X_point = [omega,-p_propre**2*np.sin(theta)] # dû à l'équa diff
    return X_point

t = np.linspace(0,T, n+1)
X = odeint(pendule, [theta0, omega0], t)
theta = X[:,0]
omega = X[:,1]
print(theta) # renvoie les valeurs successives de theta
print(omega) # renvoie les valeurs successives de omega

# Représentations graphiques
plt.figure(1)
plt.title('Elongation angulaire')
plt.plot(t, theta, 'b-')
plt.ylabel('theta en rad')
plt.xlabel('temps en s')
plt.grid()

plt.figure(2)
plt.title('Vitesse angulaire')
plt.plot(t, omega, 'g-')
plt.ylabel('omega en rad/s')
plt.xlabel('temps en s')
plt.grid()

plt.show()


## Script 2 : prise en compte de l'amortissement Q = 5 avec les mêmes conditions initiales précédentes

import numpy as np
import matplotlib.pyplot as plt
from math import *
from scipy.integrate import odeint

# Données du problème
p_propre = 1   #  pulsation propre de 1 rad/s
theta0 = 0  # CI sur l'angle initial successivement 0.25 rd, 2.5 rd, 0 rd
omega0 = 2.2  # CI sur la vitesse angulaire initiale successivement 0 rad/s, 0 rad/s , 2.2 rad/s

# Données pour l'observation
T = 20 # durée de l'étude
n = 200

# Calculs
def pendule(X,t):
    theta = X[0] # première ligne de X(t)
    omega = X[1] # seconde ligne de X(t)
    X_point = [omega,-p_propre**2*np.sin(theta)-p_propre*omega/5] # dû à l'équa diff
    return X_point

t = np.linspace(0,T, n+1)
X = odeint(pendule, [theta0, omega0], t)
theta = X[:,0]
omega = X[:,1]
print(theta) # renvoie les valeurs successives de theta
print(omega) # renvoie les valeurs successives de omega

# Représentations graphiques
plt.figure(1)
plt.title('Elongation angulaire')
plt.plot(t, theta, 'b-')
plt.ylabel('theta en rad')
plt.xlabel('temps en s')
plt.grid()

plt.figure(2)
plt.title('Vitesse angulaire')
plt.plot(t, omega, 'g-')
plt.ylabel('omega en rad/s')
plt.xlabel('temps en s')
plt.grid()

plt.show()