Soit l'équation de Mathieu :
x¨+[delta+epsilon*cos(t)]x=0
*Ecrire un programme python qui resoud numriquement cette equation pour plusieurs couples (delta,epsilon) judicieusement choisi. On intègrera sur une seule période de forcage qui sera à déterminer (seule 1 valeur est possible et elle est déterminer)
*Toujours sur le programme faire une systématisation (diagramme de phase des oscillations pour chaque couple (delta,epsilon )
*Toujours sur le programme à l'aide de la méthode de dichotomie, trouver le point à partir duquel on quitte d'une zone stable à une zone instable
*faire une discussion et une interprétation
*Conclure
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# Définition de la fonction f(t, x)
def f(t, x):
return np.array([-x[1], -δ*x[0] - ε*np.cos(t)*x[0]])
# Définition des paramètres
delta_values = [0.1, 0.5, 1.0]
epsilon_values = [0.1, 0.5, 1.0]
# Résolution numérique de l'équation pour chaque couple (δ, ε)
for delta in delta_values:
for epsilon in epsilon_values:
# Définition des conditions initiales
x0 = np.array([1.0, 0.0])
# Définition de la période de forçage
T = 2*np.pi
# Résolution numérique de l'équation
t_span = [0.0, T]
sol = solve_ivp(f, t_span, x0, rtol=1e-5, atol=1e-6)
# Extraction des solutions
t = sol.t
x = sol.y[0]
# Diagramme de phase
plt.figure()
plt.plot(x, sol.y[1])
plt.xlabel('x')
plt.ylabel('dx/dt')
plt.title('Diagramme de phase (δ = {}, ε = {})'.format(delta, epsilon))
plt.show()
# Recherche du point de bifurcation
tol = 1e-6
x_bif = 0.0
while abs(x[int(T/2)] - x[int(T/2) + 1]) > tol:
x_bif = (x[int(T/2)] + x[int(T/2) + 1]) / 2
print('Point de bifurcation (δ = {}, ε = {}) : x = {}'.format(delta, epsilon, x_bif))
- Badges
- Science topic
- Similar topics
- Engineering
- Modeling