#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Dec 13 15:28:22 2024

@author: laurentdumas
"""
import numpy as np
import matplotlib.pyplot as plt

omega = 1/2  
epsilon = 0.1  
t_max = 200    
dt = 0.01      
theta0 = np.pi / 8  
dtheta0 = 0.0      


def fonction_pendule(t, Y, omega, epsilon):
    theta, dtheta = Y  
    d2theta = -omega**2 * (1 + epsilon * np.cos(t)) * theta
    return np.array([dtheta, d2theta])

def runge_kutta4(f, t0, Y0, t_max, dt, omega, epsilon):
    t = np.arange(t0, t_max, dt)
    Y = np.zeros((len(t), len(Y0)))
    Y[0] = Y0
    
    for i in range(1, len(t)) : 
        k1 = dt * f(t[i-1], Y[i-1], omega, epsilon)
        k2 = dt * f(t[i-1] + dt/2, Y[i-1] + k1/2, omega, epsilon)
        k3 = dt * f(t[i-1] + dt/2, Y[i-1] + k2/2, omega, epsilon)
        k4 = dt * f(t[i-1] + dt, Y[i-1] + k3, omega, epsilon)
        Y[i] = Y[i-1] + (k1 + 2*k2 + 2*k3 + k4) / 6
        
    return t, Y

Y0 = [theta0, dtheta0]
t, Y = runge_kutta4(fonction_pendule, 0, Y0, t_max, dt, omega, epsilon)

plt.figure(figsize=(10, 5))
plt.plot(t, Y[:, 0], label=r"$\theta(t)$")
plt.title("Évolution temporelle de l'angle θ(t)")
plt.xlabel("Temps t")
plt.ylabel(r"Angle $\theta$ (rad)")
plt.grid()
plt.legend()
plt.show()