import numpy as np
import matplotlib.pyplot as plt

# Paramètres
k = 0.05  # s^-1
T_env = 25.0  # °C
T0 = 800.0  # °C

# EDO : dT/dt = f(t, T)
def f_cooling(t, T):
    return -k * (T - T_env)

# Solution analytique
def T_exact(t):
    return T_env + (T0 - T_env) * np.exp(-k * t)

def euler_method(f, t0, tf, y0, h):
    N = int((tf - t0) / h) + 1  # Nombre de points dans la grille
    t = np.linspace(t0, tf, N)  # Grille de temps
    y = np.zeros(N)             # Tableau pour stocker les solutions numériques
    y[0] = y0                   # Condition initiale
    
    # Pour chaque point de la grille, calculer la solution numérique
    for n in range(N - 1):
        y[n + 1] = y[n] + h * f(t[n], y[n])
    
    return t, y

# Résolution avec différents pas de temps
t0, tf = 0, 100  # secondes

t_euler_1, T_euler_1 = euler_method(f_cooling, t0, tf, T0, 1)
t_euler_5, T_euler_5 = euler_method(f_cooling, t0, tf, T0, 5)
t_euler_20, T_euler_20 = euler_method(f_cooling, t0, tf, T0, 20)

# Calcul de l'erreur pour différents h
h_range = np.array([20, 10, 5, 2, 1, 0.5, 0.2])
errors = []

for h in h_range:
    t_euler, T_euler = euler_method(f_cooling, t0, tf, T0, h)
    T_exact_at_tf = T_exact(tf)
    error = abs(T_euler[-1] - T_exact_at_tf)
    errors.append(error)

errors = np.array(errors)

def rk2_method(f, t0, tf, y0, h):
    """
    Résolution d'une EDO par la méthode de Runge-Kutta d'ordre 4.
    """
    N = int((tf - t0) / h) + 1
    t = np.linspace(t0, tf, N)
    y = np.zeros(N)
    y[0] = y0
    
    for n in range(N - 1):
        k1 = f(t[n], y[n])
        k2 = f(t[n] + h/2, y[n] + h/2 * k1)
        
        y[n + 1] = y[n] + h * k2
    
    return t, y

def rk4_method(f, t0, tf, y0, h):
    """
    Résolution d'une EDO par la méthode de Runge-Kutta d'ordre 4.
    """
    N = int((tf - t0) / h) + 1
    t = np.linspace(t0, tf, N)
    y = np.zeros(N)
    y[0] = y0
    
    for n in range(N - 1):
        k1 = f(t[n], y[n])
        k2 = f(t[n] + h/2, y[n] + h/2 * k1)
        k3 = f(t[n] + h/2, y[n] + h/2 * k2)
        k4 = f(t[n] + h, y[n] + h * k3)
        
        y[n + 1] = y[n] + (h / 6) * (k1 + 2*k2 + 2*k3 + k4)
    
    return t, y

# Résolution avec le même pas de temps h = 10 s
h = 10
t_euler, T_euler = euler_method(f_cooling, t0, tf, T0, h)
t_rk2, T_rk2 = rk2_method(f_cooling, t0, tf, T0, h)
t_rk4, T_rk4 = rk4_method(f_cooling, t0, tf, T0, h)

print("Erreur Euler :", abs(T_euler[-1] - T_exact(tf)))
print("Erreur RK4   :", abs(T_rk4[-1] - T_exact(tf)))

Erreur Euler : 4.465072986791238
Erreur RK4   : 0.020714566028654957

from scipy.integrate import solve_ivp

sol = solve_ivp(f, [t0, tf], [y0], method='RK45', rtol=1e-6, atol=1e-9, dense_output=True)

from scipy.integrate import solve_ivp

# Résolution du refroidissement avec solve_ivp
sol = solve_ivp(f_cooling, [t0, tf], [T0], method='RK45', 
                rtol=1e-6, atol=1e-9, dense_output=True)

print("Nombre de pas utilisés :", len(sol.t))
print("Temps final :", sol.t[-1])
print("T finale (numérique) :", sol.y[0, -1])
print("T finale (exacte)    :", T_exact(tf))
print("Erreur :", abs(sol.y[0, -1] - T_exact(tf)))

Nombre de pas utilisés : 23
Temps final : 100.0
T finale (numérique) : 30.22192017312949
T finale (exacte)    : 30.221908924291238
Erreur : 1.1248838251987081e-05

# Évaluer la solution à des temps arbitraires
t_dense = np.linspace(t0, tf, 1000)
T_dense = sol.sol(t_dense)[0]

# EDO raide
lam = -1000

def f_stiff(t, y):
    return lam * y

# Euler avec h = 0.003
t0_stiff, tf_stiff = 0, 0.01
y0_stiff = 1.0
h_unstable = 0.003

t_unstable, y_unstable = euler_method(f_stiff, t0_stiff, tf_stiff, y0_stiff, h_unstable)

print("y(tf) numérique :", y_unstable[-1])
print("y(tf) exact     :", np.exp(lam * tf_stiff))

y(tf) numérique : -8.0
y(tf) exact     : 4.5399929762484854e-05

# Avec h = 0.001 (stable)
h_stable = 0.001
t_stable, y_stable = euler_method(f_stiff, t0_stiff, tf_stiff, y0_stiff, h_stable)

t_exact_stiff = np.linspace(t0_stiff, tf_stiff, 500)
y_exact_stiff = np.exp(lam * t_exact_stiff)

def euler_backward(f, t0, tf, y0, h, tol=1e-8, max_iter=50):
    N = int((tf - t0) / h) + 1  # Nombre de points dans la grille
    t = np.linspace(t0, tf, N)  # Grille de temps
    y = np.zeros(N)             # Tableau pour stocker les solutions numériques
    y[0] = y0                   # Condition initiale
    eps = 1e-8                  # Utilise différences finies pour approximer la dérivée du résidu
    for n in range(N - 1):
        y_new = y[n] + h * f(t[n], y[n])  # Prédiction initiale (Euler explicite)
        # Newton-Raphson : résout g(y_{n+1}) = y_{n+1} - y_n - h*f(t_{n+1}, y_{n+1}) = 0
        for iteration in range(max_iter):  # Fait des itérations de Newton-Raphson
            g = y_new - y[n] - h * f(t[n + 1], y_new)  # Résidu
            # Dérivée de g par rapport à y_new g'(y) = 1 - h*f'(y), où f'(y) est la derivée de f par rapport à y_new
            f_y = (f(t[n + 1], y_new + eps) - f(t[n + 1], y_new)) / eps  # Approximation avec différences finies
            g_prime = 1 - h * f_y
            # Mise à jour de y_new avec la formule de Newton-Raphson
            y_new = y_new - g / g_prime
            # Si le résidu est suffisamment petit, on peut arrêter les itérations
            if abs(g) < tol:
                break
        # Enregistre la solution et passe au point suivant
        y[n + 1] = y_new
    return t, y

# Résolution avec Euler implicite, h = 0.005 (plus grand que la limite de stabilité explicite)
h_backward = 0.005
t_backward, y_backward = euler_backward(f_stiff, t0_stiff, tf_stiff, y0_stiff, h_backward)

print("y(tf) numérique :", y_backward[-1])
print("y(tf) exact     :", np.exp(lam * tf_stiff))

y(tf) numérique : 0.02777777777777778
y(tf) exact     : 4.5399929762484854e-05

# Résolution avec RK45 (explicite) — observe la difficulté
sol_explicit = solve_ivp(f_stiff, [t0_stiff, tf_stiff], [y0_stiff], method='RK45', rtol=1e-6, atol=1e-9)

# Résolution avec BDF (implicite)
sol_implicit = solve_ivp(f_stiff, [t0_stiff, tf_stiff], [y0_stiff], method='BDF', rtol=1e-6, atol=1e-9)

print("RK45 (explicite) — nombre de pas :", len(sol_explicit.t))
print("BDF (implicite)  — nombre de pas :", len(sol_implicit.t))
print("RK45 (explicite) — nombre d'évaluations de f :", sol_explicit.nfev)
print("BDF (implicite)  — nombre d'évaluations de f :", sol_implicit.nfev)

RK45 (explicite) — nombre de pas : 41
BDF (implicite)  — nombre de pas : 107
RK45 (explicite) — nombre d'évaluations de f : 248
BDF (implicite)  — nombre d'évaluations de f : 214

Jour	Sujet
03 mars	Fondamentaux de la programmation en Python
04 mars (matin)	Fonctions, listes approfondies, fichiers et graphiques
04 mars (après-midi)	Calcul numérique avec NumPy
11 mars (matin)	Tableaux en mémoire, matrices & précision numérique
11 mars (après-midi)	Résolution des systèmes linéaires
18 mars (matin)	Interpolation & recherche des racines
18 mars (après-midi)	Systèmes d'équations non linéaires
25 mars (matin)	Intégration numérique I
25 mars (après-midi)	Intégration numérique II (Monte Carlo)
01 avril (matin)	Équations différentielles ordinaires I
01 avril (après-midi)	Équations différentielles ordinaires II
08 avril	Examen

Horaire	Sujet
0:00 – 0:55	Méthode d'Euler (explicite)
0:55 – 1:50	Méthodes de Runge-Kutta (RK2, RK4)
1:50 – 2:00	☕ Pause
2:00 – 2:30	Pas adaptatif & `scipy.integrate.solve_ivp`
2:30 – 3:40	Stabilité : méthodes explicites vs implicites
3:40 – 3:45	Récapitulatif

Phénomène	Équation différentielle
Refroidissement (loi de Newton)	\(\frac{\mathrm{d}T}{\mathrm{d}t} = -k(T - T_{\text{env}})\)
Cinétique de transformation (Avrami)	\(\frac{\mathrm{d}\xi}{\mathrm{d}t} = n k t^{n-1}(1-\xi)\)
Fluage (loi puissance)	\(\frac{\mathrm{d}\varepsilon}{\mathrm{d}t} = A\sigma^n e^{-Q/(RT)}\)
Oscillateur amorti	\(m\frac{\mathrm{d}^2x}{\mathrm{d}t^2} + c\frac{\mathrm{d}x}{\mathrm{d}t} + kx = 0\)

Méthode	Ordre	Stabilité	Usage
Euler explicite	$\mathcal{O}(h)$	Limitée	Problèmes simples où avec un petit pas de temps inhérent
RK4	$\mathcal{O}(h^4)$	Bonne pour problèmes lisses	Cheval de bataille général
RK45 adaptatif	$\mathcal{O}(h^5)$	Automatique	`solve_ivp` par défaut
Euler implicite	$\mathcal{O}(h)$	A-stable	Problèmes raides, plasticity
BDF, Radau	Ordre élevé	A-stable	Problèmes raides complexes

USMA1Q — Méthodes Numériques¶

Séance 10 — Équations différentielles ordinaires I : fondamentaux¶

Plan du cours¶

Plan de la séance (indicatif)¶

Motivation — pourquoi résoudre des EDO numériquement ?¶

Forme générale d'une EDO¶

1 · Méthode d'Euler (explicite)¶

Exemple 1 : loi de refroidissement de Newton¶

Implémentation de la méthode d'Euler¶

Observations¶

L'erreur locale — l'erreur fait à chaque pas¶

Développement de Taylor de la vraie solution¶

Ce que fait Euler¶

L'erreur de troncature locale¶

Interprétation géométrique¶

L'erreur globale — l'erreur cumulée sur tout l'intervalle $[t_0, t_f]$¶

Convergence empirique¶

À vous de faire les exercices 1.1 et 1.2 !¶

2 · Méthodes de Runge-Kutta¶

Méthode du point milieu (RK2)¶

Analyse d'erreur de RK2 (méthode du point milieu)¶

Pourquoi RK2 est-il plus précis qu'Euler ?¶

L'intuition géométrique¶

Erreur locale vs globale¶

Méthode de Runge-Kutta d'ordre 4 (RK4)¶

Implémentation de RK4¶

Comparaison Euler vs RK2 vs RK4¶

Observation¶

À vous de faire les exercices 2.1 et 2.2 !¶

☕ Pause — 10 minutes¶

3 · Pas adaptatif et scipy.integrate.solve_ivp¶

Méthodes à pas adaptatif¶

scipy.integrate.solve_ivp¶

Exemple avec solve_ivp¶

Observation¶

Interpolation dense de la solution¶

À vous de faire les exercices 3.1 et 3.2 !¶

4 · Stabilité : méthodes explicites vs implicites¶

Exemple d'instabilité¶

Que se passe-t-il ?¶

Problèmes raides (stiff)¶

Méthode d'Euler implicite (backward)¶

Implémentation d'Euler implicite¶

Méthodes implicites dans solve_ivp¶

Observation¶

À vous de faire les exercices 4.1 et 4.2 !¶

Points clés de la séance¶

3 · Pas adaptatif et `scipy.integrate.solve_ivp`¶

`scipy.integrate.solve_ivp`¶

Exemple avec `solve_ivp`¶

Méthodes implicites dans `solve_ivp`¶