import matplotlib.pyplot as plt
import math

# =============================================
# Exercice 1 : Méthode d'Euler explicite
# =============================================

def euler(f, t0, y0, h, N):
    """Méthode d'Euler explicite : y' = f(t,y), y(t0)=y0"""
    t = [t0]
    y = [y0]
    
    for i in range(N):
        y_suiv = y[-1] + h * f(t[-1], y[-1])
        t_suiv = t[-1] + h
        t.append(t_suiv)
        y.append(y_suiv)
    
    return t, y


# Application à y' = -y, y(0)=1
def f1(t, y):
    return -y


t1, y1 = euler(f1, 0, 1, 0.1, 50)

# Tracé
plt.plot(t1, y1, 'b-o', label="Euler h=0.1", markersize=3)
plt.plot(t1, [math.exp(-t) for t in t1], 'r-', label="Solution exacte e^{-t}")
plt.grid()
plt.xlabel("t")
plt.ylabel("y(t)")
plt.title("Méthode d'Euler - y' = -y")
plt.legend()
plt.show()


# Avec h plus petit
t2, y2 = euler(f1, 0, 1, 0.01, 500)
plt.plot(t2, y2, 'b-', label="Euler h=0.01")
plt.plot(t2, [math.exp(-t) for t in t2], 'r-', label="Exacte")
plt.grid()
plt.legend()
plt.title("Euler avec h=0.01")
plt.show()


# Erreur maximale
erreur01 = max(abs(y1[i] - math.exp(-t1[i])) for i in range(len(y1)))
erreur001 = max(abs(y2[i] - math.exp(-t2[i])) for i in range(len(y2)))
print(f"Erreur max avec h=0.1  : {erreur01:.5f}")
print(f"Erreur max avec h=0.01 : {erreur001:.5f}")


# =============================================
# Exercice 2 : Oscillateur harmonique (système)
# =============================================

def euler2(t0, u0, v0, h, N):
    """Euler pour le système : u' = v, v' = -u"""
    t = [t0]
    u = [u0]
    v = [v0]
    
    for _ in range(N):
        u_suiv = u[-1] + h * v[-1]
        v_suiv = v[-1] + h * (-u[-1])
        t_suiv = t[-1] + h
        
        t.append(t_suiv)
        u.append(u_suiv)
        v.append(v_suiv)
    
    return t, u, v


t_osc, u_osc, v_osc = euler2(0, 1, 0, 0.1, 200)

plt.figure(figsize=(10,5))
plt.plot(t_osc, u_osc, 'b-', label="u(t)")
plt.plot(t_osc, v_osc, 'r-', label="v(t)")
plt.grid()
plt.legend()
plt.title("Oscillateur harmonique - u et v")
plt.show()

# Plan de phase (u,v)
plt.plot(u_osc, v_osc, 'g-')
plt.grid()
plt.xlabel("u")
plt.ylabel("v")
plt.title("Plan de phase (u,v)")
plt.axis("equal")
plt.show()


# Avec h=0.5 (instable)
t_gros, u_gros, v_gros = euler2(0, 1, 0, 0.5, 40)
plt.plot(t_gros, u_gros, label="u avec h=0.5")
plt.plot(t_gros, v_gros, label="v avec h=0.5")
plt.grid()
plt.legend()
plt.title("h=0.5 → instabilité")
plt.show()


# =============================================
# Exercice 3 : Croissance logistique + Point-milieu + Heun
# =============================================

def logistique(t, y):
    return 0.5 * y * (1 - y/100)


def euler3(f, t0, y0, h, N):          # Point-milieu
    t = [t0]
    y = [y0]
    for _ in range(N):
        k = y[-1] + (h/2) * f(t[-1], y[-1])
        y_suiv = y[-1] + h * f(t[-1] + h/2, k)
        t.append(t[-1] + h)
        y.append(y_suiv)
    return t, y


def euler4(f, t0, y0, h, N):          # Heun
    t = [t0]
    y = [y0]
    for _ in range(N):
        k1 = f(t[-1], y[-1])
        k2 = f(t[-1] + h, y[-1] + h*k1)
        y_suiv = y[-1] + (h/2)*(k1 + k2)
        t.append(t[-1] + h)
        y.append(y_suiv)
    return t, y


# Tracés
t_log, y_log = euler(logistique, 0, 10, 0.1, 200)
t_log3, y_log3 = euler3(logistique, 0, 10, 0.1, 200)
t_log4, y_log4 = euler4(logistique, 0, 10, 0.1, 200)

exact = [1000 / (10 + 90 * math.exp(-0.5*t)) for t in t_log]

plt.plot(t_log, y_log, 'b-', label="Euler explicite")
plt.plot(t_log3, y_log3, 'g-', label="Point-milieu")
plt.plot(t_log4, y_log4, 'm-', label="Heun")
plt.plot(t_log, exact, 'r--', label="Solution exacte")
plt.grid()
plt.legend()
plt.title("Croissance logistique")
plt.show()


# Différentes conditions initiales
for y0 in [5, 10, 100, 150]:
    t, y = euler(logistique, 0, y0, 0.1, 200)
    plt.plot(t, y, label=f"y0 = {y0}")
plt.grid()
plt.legend()
plt.title("Logistique - Différentes conditions initiales")
plt.show()


# =============================================
# Exercice 4 : Méthode de Newton
# =============================================

def newton(f, df, x0, N):
    """Méthode de Newton"""
    x = [x0]
    for _ in range(N):
        if abs(df(x[-1])) < 1e-12:      # éviter division par zéro
            print("Dérivée nulle !")
            break
        x_suiv = x[-1] - f(x[-1]) / df(x[-1])
        x.append(x_suiv)
    return x


# Racine de sqrt(2)
def f_sqrt2(x):
    return x**2 - 2

def df_sqrt2(x):
    return 2*x

approx = newton(f_sqrt2, df_sqrt2, 1.0, 10)
print("Approximation de √2 :", approx[-1])


# Racines du polynôme x³ - 3x + 1
def f_poly(x):
    return x**3 - 3*x + 1

def df_poly(x):
    return 3*x**2 - 3

print("\nRacines avec différents x0 :")
print("x0=0   →", newton(f_poly, df_poly, 0, 8)[-1])
print("x0=0.7 →", newton(f_poly, df_poly, 0.7, 8)[-1])
print("x0=2   →", newton(f_poly, df_poly, 2, 8)[-1])


# Tracé des itérations pour x0=0
x_iter = newton(f_poly, df_poly, 0, 8)
X = [x for x in x_iter]
Y = [f_poly(x) for x in X]

plt.plot(X, Y, 'bo-', label="Itérations de Newton")
plt.axhline(0, color='black', linestyle='--')
plt.grid()
plt.legend()
plt.title("Méthode de Newton - Convergence")
plt.show()