Chute libre : Mouvement d'un objet soumis uniquement à la force de gravité. Équation : y''(t) = -g.
Utiliser la méthode d'Euler pour résoudre le système d'équations différentielles : y'(t) = v(t), v'(t) = -g.
Soit y(t) la position, v(t) la vitesse, g = 9.8 m/s² l'accélération gravitationnelle
\(\frac{dy}{dt} = v(t)\)
\(\frac{dv}{dt} = -g\)
\(y_{n+1} = y_n + h \cdot v_n\)
\(v_{n+1} = v_n + h \cdot (-g)\)
Où h est le pas de temps
y₀ = hauteur initiale (ex: 100 m)
v₀ = vitesse initiale (ex: 0 m/s)
Initialiser y, v, t, dt, g
Tant que y > 0:
Afficher t, y, v
y = y + v * dt
v = v - g * dt
t = t + dtimport matplotlib.pyplot as plt
def chute_libre():
# Paramètres
g = 9.8
y = 100.0 # hauteur initiale
v = 0.0 # vitesse initiale
t = 0.0
dt = 0.01 # pas de temps
# Listes pour stocker les valeurs
temps = [t]
hauteurs = [y]
while y > 0:
# Méthode d'Euler
y = y + v * dt
v = v - g * dt
t = t + dt
# Stocker les valeurs
temps.append(t)
hauteurs.append(y)
return temps, hauteursLa simulation montre une décélération quadratique de la hauteur
La vitesse augmente linéairement (négligeant la résistance de l'air)
La chute libre est simulée avec succès en utilisant la méthode d'Euler. Le code permet de tracer la hauteur en fonction du temps.
• Méthode d'Euler : Approximation de la dérivée par un quotient
• Stabilité : Le pas de temps doit être suffisamment petit
• Validation : Comparer avec la solution analytique y(t) = y₀ - ½gt²
Croissance exponentielle : Phénomène où la vitesse de croissance est proportionnelle à la quantité présente. Équation : y'(t) = ky(t).
\(\frac{dy}{dt} = ky(t)\) où k est le taux de croissance
Solution analytique : \(y(t) = y_0 e^{kt}\)
\(y_{n+1} = y_n + h \cdot k \cdot y_n = y_n(1 + hk)\)
y₀ = 1 (quantité initiale)
k = 0.1 (taux de croissance)
Initialiser y, t, dt, k
Répéter n fois:
Afficher t, y
y = y * (1 + k*dt)
t = t + dtimport numpy as np
import matplotlib.pyplot as plt
def croissance_exponentielle():
k = 0.1 # taux de croissance
y = 1.0 # valeur initiale
t = 0.0
dt = 0.01 # pas de temps
n = 1000 # nombre d'itérations
temps = []
valeurs_numerique = []
valeurs_analytique = []
for i in range(n):
# Méthode d'Euler
y = y * (1 + k * dt)
t = t + dt
# Solutions
temps.append(t)
valeurs_numerique.append(y)
valeurs_analytique.append(np.exp(k * t))
return temps, valeurs_numerique, valeurs_analytiqueComparaison entre la solution numérique et la solution analytique y(t) = e^(kt)
Évaluation de l'erreur relative
Pour réduire l'erreur, on peut utiliser la méthode d'Euler améliorée ou Runge-Kutta
La simulation montre une croissance exponentielle. La méthode d'Euler approxime bien la solution analytique y(t) = e^(0.1t).
• Exponentielle : y'(t) = ky(t) ⇒ y(t) = y₀e^(kt)
• Validation : Comparer avec la solution analytique
• Précision : Réduire le pas de temps pour améliorer la précision
Modèle SIR : Modèle épidémiologique divisant la population en Susceptibles (S), Infectés (I) et Retirés (R).
\(\frac{dS}{dt} = -\beta SI\)
\(\frac{dI}{dt} = \beta SI - \gamma I\)
\(\frac{dR}{dt} = \gamma I\)
Où β est le taux de transmission et γ le taux de guérison
\(S_{n+1} = S_n - h \cdot \beta S_n I_n\)
\(I_{n+1} = I_n + h(\beta S_n I_n - \gamma I_n)\)
\(R_{n+1} = R_n + h \cdot \gamma I_n\)
S₀ = 999 (susceptibles)
I₀ = 1 (infectés)
R₀ = 0 (retirés)
β = 0.002, γ = 0.1
Initialiser S, I, R, t, dt, beta, gamma
Tant que t < Tmax:
Afficher t, S, I, R
S_new = S - dt * beta * S * I
I_new = I + dt * (beta * S * I - gamma * I)
R_new = R + dt * gamma * I
S = S_new
I = I_new
R = R_new
t = t + dtdef modele_sir():
# Paramètres
beta = 0.002 # taux de transmission
gamma = 0.1 # taux de guérison
# Conditions initiales
S = 999.0 # susceptibles
I = 1.0 # infectés
R = 0.0 # retirés
t = 0.0
dt = 0.1 # pas de temps
# Listes pour stocker les résultats
temps = [t]
susceptibles = [S]
infectes = [I]
retires = [R]
# Simulation
for i in range(500): # 500 itérations
# Méthode d'Euler
S_new = S - dt * beta * S * I
I_new = I + dt * (beta * S * I - gamma * I)
R_new = R + dt * gamma * I
# Mise à jour
S, I, R = S_new, I_new, R_new
t += dt
# Stockage
temps.append(t)
susceptibles.append(S)
infectes.append(I)
retires.append(R)
return temps, susceptibles, infectes, retiresObservation de la propagation de l'épidémie
Identification du pic d'infection
Calcul du nombre total de personnes touchées
Vérification que S + I + R reste constant (conservation de la population)
Comparaison avec des épidémies réelles
Le modèle SIR simule la propagation d'une épidémie. La simulation montre l'évolution des trois catégories de population au cours du temps.
• Conservation : S + I + R = constante
• Équations différentielles : Système couplé d'équations
• Validation : Conservation de la population totale
Modèle logistique : Modèle de croissance bornée : y'(t) = ry(t)(1 - y(t)/K), où K est la capacité maximale.
\(\frac{dy}{dt} = ry(t)\left(1 - \frac{y(t)}{K}\right)\)
Où r est le taux de croissance intrinsèque et K la capacité du milieu
\(y_{n+1} = y_n + h \cdot r \cdot y_n \left(1 - \frac{y_n}{K}\right)\)
y₀ = 10 (population initiale)
r = 0.5 (taux de croissance)
K = 1000 (capacité maximale)
\(y(t) = \frac{K}{1 + \left(\frac{K}{y_0} - 1\right)e^{-rt}}\)
import numpy as np
import matplotlib.pyplot as plt
def modele_logistique():
r = 0.5 # taux de croissance
K = 1000 # capacité maximale
y = 10.0 # population initiale
t = 0.0
dt = 0.1 # pas de temps
n = 1000 # nombre d'itérations
temps = []
valeurs_numerique = []
valeurs_analytique = []
for i in range(n):
# Méthode d'Euler
dy_dt = r * y * (1 - y/K)
y = y + dt * dy_dt
t = t + dt
# Solutions
temps.append(t)
valeurs_numerique.append(y)
# Solution analytique
y_analytique = K / (1 + (K/y0 - 1) * np.exp(-r*t))
valeurs_analytique.append(y_analytique)
return temps, valeurs_numerique, valeurs_analytiqueComparaison entre la solution numérique et analytique
Observation de la saturation à la capacité maximale
Effet du taux de croissance r
Effet de la capacité K
La simulation montre une croissance initialement exponentielle qui tend vers la capacité maximale K. La méthode d'Euler reproduit bien la solution analytique.
• Logistique : Croissance bornée par la capacité du milieu
• Saturation : La population tend vers K
• Validation : Comparaison avec la solution analytique
Loi de refroidissement de Newton : T'(t) = -k(T(t) - T∞), où T∞ est la température ambiante.
\(\frac{dT}{dt} = -k(T(t) - T_\infty)\)
Où k est le coefficient de refroidissement et T∞ la température ambiante
\(T(t) = T_\infty + (T_0 - T_\infty)e^{-kt}\)
\(T_{n+1} = T_n + h \cdot (-k)(T_n - T_\infty) = T_n - hk(T_n - T_\infty)\)
T₀ = 80°C (température initiale)
T∞ = 20°C (température ambiante)
k = 0.05 (coefficient de refroidissement)
import numpy as np
import matplotlib.pyplot as plt
def refroidissement_newton():
k = 0.05 # coefficient de refroidissement
T_inf = 20 # température ambiante
T = 80.0 # température initiale
t = 0.0
dt = 0.1 # pas de temps
n = 1000 # nombre d'itérations
temps = []
temperatures_numerique = []
temperatures_analytique = []
for i in range(n):
# Méthode d'Euler
dT_dt = -k * (T - T_inf)
T = T + dt * dT_dt
t = t + dt
# Solutions
temps.append(t)
temperatures_numerique.append(T)
T_analytique = T_inf + (80 - T_inf) * np.exp(-k * t)
temperatures_analytique.append(T_analytique)
return temps, temperatures_numerique, temperatures_analytiqueLa température diminue exponentiellement vers la température ambiante
Comparaison avec la solution analytique
Vérification que la température tend vers T∞
Calcul de la constante de temps τ = 1/k
La simulation reproduit la loi de refroidissement de Newton. La température tend exponentiellement vers la température ambiante de 20°C.
• Newton : Refroidissement proportionnel à l'écart de température
• Convergence : La température tend vers T∞
• Constante de temps : τ = 1/k
Désintégration radioactive : Processus stochastique modélisé par N'(t) = -λN(t), mais aussi simulé avec Monte Carlo.
\(\frac{dN}{dt} = -\lambda N(t)\)
Solution : \(N(t) = N_0 e^{-\lambda t}\)
Chaque noyau a une probabilité p = λΔt de se désintégrer pendant un petit intervalle Δt
Utilisation de nombres aléatoires pour décider si un noyau se désintègre
Initialiser N (nombre de noyaux), lambda, t, dt
Tant que N > 0:
Nouveau_N = N
Pour chaque noyau:
Si random() < lambda * dt:
Nouveau_N = Nouveau_N - 1
N = Nouveau_N
t = t + dtimport random
import numpy as np
def desintegration_radioactive_monte_carlo():
N = 10000 # nombre initial de noyaux
lambda_decay = 0.001 # constante de désintégration
t = 0.0
dt = 1.0 # pas de temps
n_iterations = 1000
temps = [t]
nombres_noyaux = [N]
for i in range(n_iterations):
nouveaux_noyaux = N
proba_desintegration = lambda_decay * dt
# Pour chaque noyau encore présent
for j in range(N):
if random.random() < proba_desintegration:
nouveaux_noyaux -= 1
N = nouveaux_noyaux
t += dt
# Stocker les résultats
if N > 0:
temps.append(t)
nombres_noyaux.append(N)
else:
break
return temps, nombres_noyauxSimulation multiple pour obtenir une moyenne
Comparaison avec la solution analytique N(t) = N₀e^(-λt)
Calcul de la moyenne sur plusieurs simulations
Évaluation de la variance
Concordance avec le modèle déterministe pour de grands N
Fluctuations statistiques pour de petits N
La simulation Monte Carlo reproduit le processus stochastique de désintégration radioactive. Pour de grands N, les résultats convergent vers le modèle déterministe.
• Stochastique : Processus probabiliste
• Monte Carlo : Utilisation de nombres aléatoires
• Validation : Convergence vers le modèle déterministe
Équation de diffusion : ∂C/∂t = D∇²C, où C est la concentration et D le coefficient de diffusion.
\(\frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2}\)
Discrétisation : \(\frac{C_i^{n+1} - C_i^n}{\Delta t} = D \frac{C_{i+1}^n - 2C_i^n + C_{i-1}^n}{(\Delta x)^2}\)
\(C_i^{n+1} = C_i^n + \frac{D \Delta t}{(\Delta x)^2}(C_{i+1}^n - 2C_i^n + C_{i-1}^n)\)
Concentration initiale : gaussienne ou impulsionnelle
Conditions aux limites : nulles ou périodiques
Condition de stabilité : \(\frac{D \Delta t}{(\Delta x)^2} \leq \frac{1}{2}\)
import numpy as np
import matplotlib.pyplot as plt
def diffusion_1d():
# Paramètres
D = 0.1 # coefficient de diffusion
L = 10 # longueur du domaine
nx = 100 # nombre de points spatiaux
dx = L / (nx - 1)
dt = 0.1 # pas de temps
nt = 1000 # nombre de pas temporels
# Condition initiale (gaussienne)
x = np.linspace(0, L, nx)
C = np.exp(-(x - L/2)**2 / 0.5)
# Coefficient pour le schéma
r = D * dt / (dx**2)
# Vérification de la condition de stabilité
assert r <= 0.5, "Schéma instable!"
# Simulation
for n in range(nt):
C_new = C.copy()
# Mise à jour des points internes
for i in range(1, nx-1):
C_new[i] = C[i] + r * (C[i+1] - 2*C[i] + C[i-1])
C = C_new
return x, CObservation de la propagation de la concentration
Vérification de la conservation de la masse
Modèle 2D ou 3D
Conditions aux limites différentes
La simulation montre la diffusion d'une substance dans un milieu 1D. La concentration se répartit uniformément au fil du temps.
• Stabilité : Respecter la condition CFL
• Différence finie : Approximation des dérivées partielles
• Conservation : Vérifier la conservation de la masse
Circuit RC : Circuit composé d'une résistance R et d'un condensateur C. Équation : RC·dq/dt + q = E.
\(RC\frac{dq}{dt} + q = CE\)
Soit \(\frac{dq}{dt} = \frac{E - q/C}{R} = \frac{1}{RC}(CE - q)\)
\(q(t) = CE(1 - e^{-t/RC})\)
La constante de temps est τ = RC
\(q_{n+1} = q_n + h \cdot \frac{1}{RC}(CE - q_n)\)
q₀ = 0 (condensateur déchargé)
E = 12V (tension d'alimentation)
R = 1000Ω, C = 0.001F
import numpy as np
def charge_condensateur_rc():
R = 1000.0 # résistance en ohms
C = 0.001 # capacité en farads
E = 12.0 # tension d'alimentation en volts
tau = R * C # constante de temps
q = 0.0 # charge initiale
t = 0.0
dt = 0.001 # pas de temps
n = 10000 # nombre d'itérations
temps = []
charges_numerique = []
charges_analytique = []
for i in range(n):
# Méthode d'Euler
dq_dt = (E - q/C) / R
q = q + dt * dq_dt
t = t + dt
# Solutions
temps.append(t)
charges_numerique.append(q)
charges_analytique.append(C * E * (1 - np.exp(-t/tau)))
return temps, charges_numerique, charges_analytiqueLa tension aux bornes du condensateur est V = q/C
V(t) = E(1 - e^(-t/RC))
Le condensateur se charge exponentiellement
Temps de charge ≈ 5τ pour atteindre 99% de la tension finale
La simulation montre la charge exponentielle du condensateur. La tension tend vers la tension d'alimentation E.
• Constante de temps : τ = RC
• Chargement : Exponentiel vers la tension finale
• Validation : Comparaison avec la solution analytique
Lois de Kepler : Les orbites planétaires sont des ellipses avec le Soleil à un foyer. Équations : F = -GMm/r².
\(\frac{d^2\vec{r}}{dt^2} = -\frac{GM}{r^3}\vec{r}\)
Où G est la constante gravitationnelle, M la masse du Soleil, r la distance planète-Soleil
\(\frac{d^2x}{dt^2} = -\frac{GMx}{(x^2 + y^2)^{3/2}}\)
\(\frac{d^2y}{dt^2} = -\frac{GMy}{(x^2 + y^2)^{3/2}}\)
Convertir en système d'équations du premier ordre :
\(v_x = \frac{dx}{dt}\), \(v_y = \frac{dy}{dt}\)
\(\frac{dv_x}{dt} = -\frac{GMx}{(x^2 + y^2)^{3/2}}\)
\(\frac{dv_y}{dt} = -\frac{GMy}{(x^2 + y^2)^{3/2}}\)
Position : (x₀, y₀) = (1 AU, 0)
Vitesse : (vₓ₀, vᵧ₀) = (0, 29.8 km/s)
GM = 1.327 × 10¹¹ km³/s² (paramètre gravitationnel solaire)
import numpy as np
def orbite_planetaire():
GM = 1.327e11 # paramètre gravitationnel solaire (km³/s²)
# Conditions initiales (Terre)
x, y = 1.496e8, 0.0 # position initiale (km)
vx, vy = 0.0, 29.8 # vitesse initiale (km/s)
t = 0.0
dt = 86400.0 # pas de temps (1 jour en secondes)
n = 365 # nombre d'itérations (1 an)
# Listes pour stocker les positions
x_list = [x]
y_list = [y]
for i in range(n):
# Calcul de la distance
r = np.sqrt(x**2 + y**2)
r3 = r**3
# Calcul des forces
ax = -GM * x / r3
ay = -GM * y / r3
# Mise à jour des vitesses
vx = vx + ax * dt
vy = vy + ay * dt
# Mise à jour des positions
x = x + vx * dt
y = y + vy * dt
t = t + dt
# Stockage des positions
x_list.append(x)
y_list.append(y)
return x_list, y_listUtilisation de la méthode de Verlet pour une meilleure stabilité
Conservation de l'énergie mécanique
Vérification de la conservation de l'énergie
Calcul de l'excentricité de l'orbite
La simulation reproduit l'orbite elliptique d'une planète autour du Soleil selon les lois de la gravitation universelle.
• Gravitation : Force inversement proportionnelle au carré de la distance
• Stabilité : Choisir un pas de temps suffisamment petit
• Validation : Conservation de l'énergie mécanique
Équations de Navier-Stokes : Modélisent l'écoulement des fluides visqueux. Simplification : équation de Poiseuille pour l'écoulement laminaire.
Écoulement laminaire dans un tuyau cylindrique
Vitesse radiale : \(v(r) = v_{max}(1 - \frac{r^2}{R^2})\)
Où R est le rayon du tuyau
Division du tuyau en sections annulaires
Calcul de la vitesse dans chaque section
Le débit est constant : \(Q = \int_0^R v(r) \cdot 2\pi r dr\)
import numpy as np
def ecoulement_fluide_tuyau():
R = 0.01 # rayon du tuyau (m)
vmax = 2.0 # vitesse maximale (m/s)
n_points = 100 # nombre de points radiaux
# Positions radiales
r_values = np.linspace(0, R, n_points)
# Vitesses correspondantes
v_values = vmax * (1 - (r_values/R)**2)
# Calcul du débit
delta_r = R / n_points
debit = 0
for i in range(n_points):
debit += v_values[i] * 2 * np.pi * r_values[i] * delta_r
# Simulation de particules dans l'écoulement
n_particules = 1000
positions = []
vitesses = []
for i in range(n_particules):
# Position radiale aléatoire
r_particule = np.random.uniform(0, R)
# Vitesse correspondante
v_particule = vmax * (1 - (r_particule/R)**2)
positions.append(r_particule)
vitesses.append(v_particule)
return r_values, v_values, debit, positions, vitessesLe profil est parabolique avec vitesse maximale au centre
Vitesse nulle au bord du tuyau (condition d'adhérence)
\(Q = \frac{\pi R^4 \Delta P}{8 \mu L}\) (loi de Poiseuille)
Vérification du profil parabolique
Calcul de la vitesse moyenne
La simulation montre le profil parabolique de vitesse dans un écoulement laminaire. La vitesse est maximale au centre du tuyau.
• Loi de Poiseuille : Écoulement laminaire dans un tube
• Profil parabolique : Vitesse maximale au centre
• Adhérence : Vitesse nulle au bord
