Code FFT Glissante
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import wave
import os
from pygments import pigmentize
from matplotlib.animation import PillowWriter
# ----------------------------
# PARAMÈTRE : chemin audio
# ----------------------------
# filename = "audio.wav"
filename = "no_audio.wav"
# ----------------------------
# Chargement ou génération du signal
# ----------------------------
if os.path.exists(filename):
wav = wave.open(filename, 'rb')
fs = wav.getframerate()
n_samples = wav.getnframes()
audio = wav.readframes(n_samples)
signal = np.frombuffer(audio, dtype=np.int16).astype(float)
if wav.getnchannels() == 2:
signal = signal[::2]
signal /= np.max(np.abs(signal))
wav.close()
else:
fs = 4400
duration = 3
t = np.arange(0, duration, 1/fs)
freqs_notes = [262, 2300, 392, 1230]
signal = np.zeros_like(t)
segment = len(t) // len(freqs_notes)
for i, f in enumerate(freqs_notes):
start = i * segment
end = (i + 1) * segment
signal[start:end] = np.sin(2*np.pi*f*t[start:end])
t = np.arange(len(signal)) / fs
# ----------------------------
# PARAMÈTRES FFT
# ----------------------------
window_size = 1024
hop_size = 50
window = np.hanning(window_size)
freqs = np.fft.rfftfreq(window_size, 1/fs)
# ----------------------------
# BANQUE DE 22 FILTRES RECTANGULAIRES
# ----------------------------
f_min = 20
#f_max = 20000
f_max = fs
n_filters = 22
width_min = 20
width_max = 500
centers = np.linspace(f_min, f_max, n_filters)
widths = np.linspace(width_min, width_max, n_filters)
def build_filterbank(freqs):
filters = np.zeros((n_filters, len(freqs)))
for i in range(n_filters):
f0 = centers[i]
w = widths[i] / 2
mask = (freqs >= f0 - w) & (freqs <= f0 + w)
filters[i, mask] = 1.0
return filters
filterbank = build_filterbank(freqs)
# ----------------------------
# FIGURE : 4 FENÊTRES
# ----------------------------
fig, (ax_global, ax_time, ax_fft, ax_filt) = plt.subplots(4, 1, figsize=(12, 12), dpi = 60)
# -------------------------------------------------
# 1) SIGNAL GLOBAL + FENÊTRE GLISSANTE
# -------------------------------------------------
ax_global.set_title("Signal global avec fenêtre glissante")
ax_global.plot(t, signal, color='steelblue', lw=1)
ax_global.set_xlim(0, t[-1])
ax_global.set_ylim(-1.1, 1.1)
# rectangle représentant la fenêtre
window_patch = ax_global.axvspan(
0,
window_size/fs,
color='red',
alpha=0.3
)
# -------------------------------------------------
# 2) SIGNAL LOCAL
# -------------------------------------------------
ax_time.set_title("Signal local analysé")
ax_time.set_xlim(0, window_size/fs)
ax_time.set_ylim(-1.1, 1.1)
line_time, = ax_time.plot([], [], lw=1.5)
# -------------------------------------------------
# 3) FFT GLISSANTE
# -------------------------------------------------
ax_fft.set_title("FFT glissante")
ax_fft.set_xlim(0, np.max(freqs))
ax_fft.set_ylim(0, 1)
line_fft, = ax_fft.plot([], [], lw=1.5, color='orange')
# -------------------------------------------------
# 4) ÉNERGIE FILTRÉE
# -------------------------------------------------
ax_filt.set_title("Énergie dans les 22 bandes fréquentielles")
ax_filt.set_xlim(-0.5, n_filters - 0.5)
ax_filt.set_ylim(0, 1)
bars = ax_filt.bar(np.arange(n_filters), np.zeros(n_filters))
plt.tight_layout()
# ----------------------------
# ANIMATION
# ----------------------------
def update(frame):
start = frame * hop_size
end = start + window_size
if end >= len(signal):
return line_time, line_fft, bars
# segment courant
segment = signal[start:end]
# fenêtrage
windowed = segment * window
# FFT
spectrum = np.abs(np.fft.rfft(windowed))
spectrum /= np.max(spectrum) + 1e-12
# -------------------------------------------------
# 1) MISE À JOUR FENÊTRE GLISSANTE
# -------------------------------------------------
t_start = start / fs
t_end = end / fs
window_patch.set_xy([
[t_start, -1.1],
[t_start, 1.1],
[t_end, 1.1],
[t_end, -1.1],
[t_start, -1.1]
])
# -------------------------------------------------
# 2) SIGNAL LOCAL
# -------------------------------------------------
time_axis = np.arange(window_size) / fs
line_time.set_data(time_axis, segment)
# -------------------------------------------------
# 3) FFT
# -------------------------------------------------
line_fft.set_data(freqs, spectrum)
# -------------------------------------------------
# 4) ÉNERGIES FILTRÉES
# -------------------------------------------------
energies = filterbank @ spectrum
energies /= np.max(energies) + 1e-12
for i, b in enumerate(bars):
b.set_height(energies[i])
return line_time, line_fft, bars
# animation
frames = (len(signal) - window_size) // hop_size
ani = animation.FuncAnimation(
fig,
update,
frames=frames,
interval=30,
blit=False
)
# CODE POUR SAUVEGARDER DES GIF PAS TROP VOLUMINEUX
# METTRE LE dpi = 60 DANS LES FIGURES
# from matplotlib.animation import PillowWriter
# ani = FuncAnimation(fig, update, frames=800, interval=20)
ani.save("FFT_GLISSANTE.gif",
writer=PillowWriter(fps=20)) #, save_count=80)
plt.show()
# pygmentize -f pdf -o FFT_Glissante_v3.pdf FFT_Glissante_v3.pyCode Aliasing Shannon
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter
# =========================
# PARAMÈTRES
# =========================
f_signal = 50
T = 0.2
fe_ok = 500
fe_bad = 70
t_cont = np.linspace(0, T, 5000)
t_ok = np.arange(0, T, 1/fe_ok)
t_bad = np.arange(0, T, 1/fe_bad)
frames = 120
# =========================
# FIGURE
# =========================
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 9), dpi = 60)
# =========================
# 1. TEMPS + SHANNON / NYQUIST
# =========================
ax1.plot(t_cont, np.sin(2*np.pi*f_signal*t_cont), 'k', alpha=0.3)
line_ok, = ax1.plot([], [], 'bo-', label="Sans aliasing")
line_bad, = ax1.plot([], [], 'ro-', label="Avec aliasing")
nyquist_line = ax1.axhline(0, color='green', linestyle='--', label="Nyquist (fe/2)")
ax1.set_xlim(0, T)
ax1.set_ylim(-1.2, 1.2)
ax1.set_title("Visualisation de l'aliasing (fsignal = 50 Hz, fech1 = 500 Hz, fech2 = 70 Hz")
ax1.legend()
ax1.grid()
# =========================
# 2. SPECTRE + REPLIEMENT
# =========================
def spectrum(x, fe):
N = len(x)
X = np.fft.fft(x)
f = np.fft.fftfreq(N, 1/fe)
return f[:N//2], np.abs(X[:N//2]) / N
# spectres fixes
x_ok = np.sin(2*np.pi*f_signal*t_ok)
x_bad = np.sin(2*np.pi*f_signal*t_bad)
f1, X1 = spectrum(x_ok, fe_ok)
f2, X2 = spectrum(x_bad, fe_bad)
stem_ok = ax2.stem(f1, X1, linefmt='b-', markerfmt='bo', basefmt=" ", label="OK")
stem_bad = ax2.stem(f2, X2, linefmt='r-', markerfmt='ro', basefmt=" ", label="Aliasing")
f_max_line = ax2.axvline(f_signal, color='black', linestyle='--', label="f_max")
alias_line = ax2.axvline(0, color='red', linestyle='--', label="f_alias")
ax2.set_xlim(0, 200)
ax2.set_ylim(0, 0.6)
ax2.set_title("Repliement spectral")
ax2.grid()
ax2.legend()
# =========================
# UPDATE
# =========================
def update(i):
# -------------------------
# temps
# -------------------------
idx_ok = max(1, int(len(t_ok)*i/frames))
idx_bad = max(1, int(len(t_bad)*i/frames))
line_ok.set_data(t_ok[:idx_ok], np.sin(2*np.pi*f_signal*t_ok[:idx_ok]))
line_bad.set_data(t_bad[:idx_bad], np.sin(2*np.pi*f_signal*t_bad[:idx_bad]))
# -------------------------
# SHANNON DYNAMIQUE
# on fait varier fe_bad progressivement
# -------------------------
fe_dynamic = np.linspace(40, 200, frames)[i]
nyquist = fe_dynamic / 2
shannon_limit = 2 * f_signal
# ligne Nyquist
nyquist_line.set_ydata([0, 0]) # stable visuellement
nyquist_line.set_label(f"Nyquist = {nyquist:.1f} Hz")
# -------------------------
# SPECTRE + ALIAS
# -------------------------
fe_current = np.linspace(40, 200, frames)[i]
if f_signal > fe_current / 2:
f_alias = abs(f_signal - fe_current)
else:
f_alias = f_signal
alias_line.set_xdata([f_alias, f_alias])
return line_ok, line_bad, nyquist_line, alias_line
ani = FuncAnimation(fig, update, frames=frames, interval=50, blit=False)
# CODE POUR SAUVEGARDER DES GIF PAS TROP VOLUMINEUX
# METTRE LE dpi = 60 DANS LES FIGURES
# from matplotlib.animation import PillowWriter
# ani = FuncAnimation(fig, update, frames=80, interval=20)
# ani.save("Shannon.gif",
# writer=PillowWriter(fps=20)) #, save_count=80)
plt.tight_layout()
plt.show()Code SHFT
import numpy as np
import matplotlib.pyplot as plt
import wave
import os
# ----------------------------
# PARAMÈTRE : chemin audio
# ----------------------------
filename = "audio.wav" # <-- change ici
# ----------------------------
# Essai de chargement du fichier audio
# ----------------------------
if os.path.exists(filename):
print("Fichier audio trouvé : utilisation du fichier")
wav = wave.open(filename, 'rb')
fs = wav.getframerate()
n_samples = wav.getnframes()
audio = wav.readframes(n_samples)
signal = np.frombuffer(audio, dtype=np.int16)
# Normalisation
signal = signal / np.max(np.abs(signal))
# Conversion stéréo → mono
if wav.getnchannels() == 2:
signal = signal[::2]
wav.close()
t = np.arange(len(signal)) / fs
else:
print("Fichier non trouvé → génération d'un signal musical")
fs = 1000
duration = 5
t = np.arange(0, duration, 1/fs)
freqs_notes = [262, 330, 392, 523] # Do - Mi - Sol - Do
signal = np.zeros_like(t)
segment = len(t) // len(freqs_notes)
for i, f in enumerate(freqs_notes):
start = i * segment
end = (i + 1) * segment
signal[start:end] = np.sin(2*np.pi*f*t[start:end])
# ----------------------------
# STFT PARAMÈTRES
# ----------------------------
window_size = 512
hop_size = 128
window = np.hanning(window_size)
spectrogram = []
times = []
# ----------------------------
# STFT glissante
# ----------------------------
for start in range(0, len(signal) - window_size, hop_size):
segment = signal[start:start + window_size]
windowed = segment * window
fft = np.fft.rfft(windowed)
spectrogram.append(np.abs(fft))
times.append(start / fs)
spectrogram = np.array(spectrogram).T
freqs = np.fft.rfftfreq(window_size, 1/fs)
# ----------------------------
# AFFICHAGE
# ----------------------------
# 1) Signal temporel
plt.figure(figsize=(12, 4))
plt.plot(t, signal)
plt.xlabel("Temps (s)")
plt.ylabel("Amplitude")
plt.title("Signal temporel")
plt.grid(True)
# 2) Spectrogramme
plt.figure(figsize=(12, 6))
plt.pcolormesh(times, freqs, spectrogram, shading='auto')
plt.xlabel("Temps (s)")
plt.ylabel("Fréquence (Hz)")
plt.title("Spectrogramme (STFT)")
plt.colorbar(label="Amplitude")
plt.ylim(0, fs/2)
plt.show()
Vibration membrane cochlée
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter
# ==================================================
# PARAMETRES
# ==================================================
c = 1.0
k = 8.0
omega = c * k
L = 10.0
sigma = 1.2
x01 = -5.0
delay = 8.0
x02 = x01 - c * delay
r = 1.0
x = np.linspace(-15, 15, 1200)
# ==================================================
# PAQUET INCIDENT
# ==================================================
def packet(x, t, x0):
envelope = np.exp(
-((x - x0 - c*t)**2)/(2*sigma**2)
)
carrier = np.cos(k*x - omega*t)
return envelope * carrier
# ==================================================
# PAQUET REFLECHI
# ==================================================
def reflected_packet(x, t, x0):
tc = (L - x0)/c
if t < tc:
return np.zeros_like(x)
envelope = np.exp(
-((x - (2*L - x0 - c*t))**2)/(2*sigma**2)
)
carrier = np.cos(-k*x - omega*t)
return r * envelope * carrier
# ==================================================
# MEMBRANE BASILAIRE SIMPLIFIEE
# ==================================================
x_bf = 4 # position résonante
sigma_bf = 0.8 # largeur de la zone active
gain_bf = np.exp(
-((x - x_bf)**2)/(2*sigma_bf**2)
)
# ==================================================
# FIGURE
# ==================================================
fig, (ax1, ax2) = plt.subplots(
2, 1,
figsize=(10, 8)
)
# ------------------------
# Panneau supérieur
# ------------------------
line_total, = ax1.plot(
[], [],
color='navy',
lw=2
)
ax1.axvline(
L,
color='red',
linestyle=':',
label='Interface'
)
ax1.set_xlim(-15, 15)
ax1.set_ylim(-2.5, 2.5)
ax1.set_ylabel("Amplitude")
ax1.set_title("Propagation du paquet d'onde")
ax1.legend()
# ------------------------
# Panneau inférieur
# ------------------------
line_bm, = ax2.plot(
[],
[],
color='darkgreen',
lw=2
)
ax2.axvline(
x_bf,
color='orange',
linestyle='--',
label='Position caractéristique'
)
ax2.set_xlim(-15, 15)
ax2.set_ylim(-2.5, 2.5)
ax2.set_xlabel("Position cochléaire")
ax2.set_ylabel("Déplacement")
ax2.set_title(
"Réponse localisée de type membrane basilaire"
)
ax2.legend()
# ==================================================
# ANIMATION
# ==================================================
def update(frame):
t = frame * 0.08
psi1 = packet(x, t, x01)
psi2 = packet(x, t, x02)
refl = (
reflected_packet(x, t, x01)
+ reflected_packet(x, t, x02)
)
total = psi1 + psi2 + refl
# ------------------------------------------------
# Réponse cochléaire simplifiée
# ------------------------------------------------
membrane = gain_bf * total
line_total.set_data(x, total)
line_bm.set_data(x, membrane)
ax1.set_title(
f"Propagation du paquet - t = {t:.2f}"
)
return line_total, line_bm
ani = FuncAnimation(
fig,
update,
frames=500,
interval=30,
blit=True
)
# CODE POUR SAUVEGARDER DES GIF PAS TROP VOLUMINEUX
# METTRE LE dpi = 60 DANS LES FIGURES
# from matplotlib.animation import PillowWriter
# ani = FuncAnimation(fig, update, frames=800, interval=20)
# ani.save("Cochlée_Paquet_Ondes.gif",
# writer=PillowWriter(fps=20)) #, save_count=80)
plt.tight_layout()
plt.show()
# ani.save('Vibration_Mecanique_Membrane_Cochlee_4.gif', writer='pillow', fps=30)Etrier Membrane cochlée
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# --- Paramètres ---
L_membrane = 0.5
L_water = 20.0
c_water = 148.2
A = 0.5
freq = 2
Nx = 200
Nt = 500
dx = L_water / Nx
dt = 0.02
# --- grille eau ---
x_water = np.linspace(L_membrane, L_membrane + L_water, Nx)
# --- état bille ---
ball_x = 0.2
ball_speed = 0.0
ball_v = 0.02
# --- clavier ---
keys = {"left": False, "right": False}
def on_key_press(event):
if event.key == "left":
keys["left"] = True
if event.key == "right":
keys["right"] = True
def on_key_release(event):
if event.key == "left":
keys["left"] = False
if event.key == "right":
keys["right"] = False
# --- onde eau (simple propagation) ---
def water_wave(x, t):
return A * np.sin(2 * np.pi * freq * (t - (x - L_membrane) / c_water))
# --- figure ---
fig, ax = plt.subplots(figsize=(10, 6))
line_water, = ax.plot([], [], lw=2, color='green', label="Eau")
# membrane verticale (ligne pointillée)
membrane_line = ax.axvline(L_membrane, color='k', linestyle='--')
# membrane "bleue" (petit segment déformé horizontalement)
membrane_top, = ax.plot([], [], lw=4, color='blue')
# bille
ball, = ax.plot([], [], 'ro', markersize=10)
ax.set_xlim(0, L_membrane + L_water)
ax.set_ylim(-1.5, 1.5)
ax.set_title("Membrane + onde + bille contrôlée")
ax.legend()
fig.canvas.mpl_connect("key_press_event", on_key_press)
fig.canvas.mpl_connect("key_release_event", on_key_release)
# --- init ---
def init():
line_water.set_data([], [])
membrane_top.set_data([], [])
ball.set_data([], [])
return line_water, membrane_top, ball
# --- update ---
def update(frame):
global ball_x, ball_speed
t = frame * dt
# -----------------------
# onde eau
# -----------------------
y_water = water_wave(x_water, t)
line_water.set_data(x_water, y_water)
# -----------------------
# membrane (ligne verticale + petite déformation horizontale)
# -----------------------
y_mem = np.linspace(-1, 1, 50)
x_mem = L_membrane + 0.05 * np.sin(2*np.pi*freq*t) # vibration horizontale
membrane_top.set_data(
[x_mem]*len(y_mem),
y_mem
)
# -----------------------
# bille (mouvement horizontal)
# -----------------------
if keys["left"]:
ball_x -= ball_v
if keys["right"]:
ball_x += ball_v
ball_x = np.clip(ball_x, 0, L_membrane + L_water)
ball.set_data([ball_x], [0])
return line_water, membrane_top, ball
ani = FuncAnimation(fig, update, frames=Nt, init_func=init, blit=True, interval=20)
plt.show()Potentiel d’Action
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.widgets import Slider
from matplotlib.animation import PillowWriter
# -------------------------
# Paramètres
# -------------------------
L = 30
Nx = 500
x = np.linspace(0, L, Nx)
v = 10.0
tau_spike = 0.2
dt = 0.05
lambda_rate = 3.0
tau_ref = 0.08 # réfractaire
# -------------------------
# état
# -------------------------
spikes = []
last_spike_time = -1e9
def poisson_event(lam, dt):
return np.random.rand() < lam * dt
def spike_shape(x, t, t0):
return np.exp(-((x - v*(t - t0))**2) / (2*tau_spike**2))
# -------------------------
# figure
# -------------------------
fig, ax = plt.subplots(dpi = 60)
plt.subplots_adjust(bottom=0.25)
line, = ax.plot(x, np.zeros_like(x))
ax.set_ylim(0, 1.5)
ax.set_xlim(0, L)
ax.set_title("Neurone avec période réfractaire")
ax.set_xlabel("temps(ms)")
ax.set_ylabel("Potentiel membranaire (mV)")
# -------------------------
# SLIDER λ
# -------------------------
ax_lambda = plt.axes([0.25, 0.1, 0.5, 0.03])
slider_lambda = Slider(ax_lambda, "λ (fréquence)", 0.1, 20.0, valinit=lambda_rate)
def update_lambda(val):
global lambda_rate
lambda_rate = slider_lambda.val
slider_lambda.on_changed(update_lambda)
# -------------------------
# update animation
# -------------------------
def update(frame):
global spikes, last_spike_time
t = frame * dt
# -------------------------
# réfractaire
# -------------------------
can_fire = (t - last_spike_time) > tau_ref
if can_fire and poisson_event(lambda_rate, dt):
spikes.append(t)
last_spike_time = t
# -------------------------
# champ
# -------------------------
y = np.zeros_like(x)
new_spikes = []
for t0 in spikes:
if t - t0 < L / v:
y += spike_shape(x, t, t0)
new_spikes.append(t0)
spikes = new_spikes
line.set_ydata(y)
return line,
ani = FuncAnimation(fig, update, frames=2000, interval=20, blit=True)
# CODE POUR SAUVEGARDER DES GIF PAS TROP VOLUMINEUX
# METTRE LE dpi = 60 DANS LES FIGURES
# from matplotlib.animation import PillowWriter
# ani = FuncAnimation(fig, update, frames=500, interval=20)
# ani.save("Potentiel_d_Action.gif",
# writer=PillowWriter(fps=20)) #, save_count=80)
plt.show()
Potentiel d Action sans periode réfractaire
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.widgets import Slider
# -------------------------
# Paramètres physiques
# -------------------------
L = 10.0 # longueur axone
Nx = 500
x = np.linspace(0, L, Nx)
v = 2.0 # vitesse propagation spike
tau_spike = 0.2 # largeur du potentiel d'action
dt = 0.05
# -------------------------
# Etat des spikes
# -------------------------
spikes = [] # chaque spike = (t0)
# -------------------------
# lambda initial (Poisson)
# -------------------------
lambda_rate = 2.0
# -------------------------
# fonction Poisson
# -------------------------
def poisson_event(lam, dt):
return np.random.rand() < lam * dt
# -------------------------
# potentiel d'action
# -------------------------
def spike_shape(x, t, t0):
return np.exp(-((x - v*(t - t0))**2) / (2*tau_spike**2))
# -------------------------
# figure
# -------------------------
fig, ax = plt.subplots()
plt.subplots_adjust(bottom=0.25)
line, = ax.plot(x, np.zeros_like(x), lw=2)
ax.set_ylim(0, 1.5)
ax.set_xlim(0, L)
ax.set_title("Train de potentiels d'action (processus de Poisson)")
# -------------------------
# Slider λ
# -------------------------
ax_lambda = plt.axes([0.25, 0.1, 0.5, 0.03])
slider_lambda = Slider(ax_lambda, 'λ (fréquence spikes)', 0.1, 20.0, valinit=lambda_rate)
def update_lambda(val):
global lambda_rate
lambda_rate = slider_lambda.val
slider_lambda.on_changed(update_lambda)
# -------------------------
# update animation
# -------------------------
def update(frame):
global spikes
t = frame * dt
# -------------------------
# génération Poisson
# -------------------------
if poisson_event(lambda_rate, dt):
spikes.append(t)
# -------------------------
# champ total
# -------------------------
y = np.zeros_like(x)
new_spikes = []
for t0 in spikes:
# ne garder que les spikes visibles
if t - t0 < L / v:
y += spike_shape(x, t, t0)
new_spikes.append(t0)
spikes = new_spikes
line.set_ydata(y)
return line,
ani = FuncAnimation(fig, update, frames=2000, interval=20, blit=True)
plt.show()Comparaison Euler RK4
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.signal import find_peaks
from scipy import stats
from matplotlib.animation import PillowWriter
# ==========================================================
# PARAMETRES
# ==========================================================
g = 9.81
L = 1.0
dt = 0.01
T = 25
N = int(T / dt)
t = np.linspace(0, T, N)
# ==========================================================
# CONDITIONS INITIALES
# ==========================================================
theta_rk4 = np.zeros(N)
theta_rk2 = np.zeros(N)
theta_euler = np.zeros(N)
omega_rk4 = np.zeros(N)
omega_rk2 = np.zeros(N)
omega_euler = np.zeros(N)
theta_rk4[0] = np.pi/4
theta_rk2[0] = np.pi/4
theta_euler[0] = np.pi/4
# ==========================================================
# SYSTEME
# ==========================================================
def f(theta, omega):
return omega, -(g/L)*np.sin(theta)
# ==========================================================
# RK4
# ==========================================================
for i in range(N-1):
th, om = theta_rk4[i], omega_rk4[i]
k1_t, k1_w = f(th, om)
k2_t, k2_w = f(th+0.5*dt*k1_t, om+0.5*dt*k1_w)
k3_t, k3_w = f(th+0.5*dt*k2_t, om+0.5*dt*k2_w)
k4_t, k4_w = f(th+dt*k3_t, om+dt*k3_w)
theta_rk4[i+1] = th + dt/6*(k1_t+2*k2_t+2*k3_t+k4_t)
omega_rk4[i+1] = om + dt/6*(k1_w+2*k2_w+2*k3_w+k4_w)
# ==========================================================
# RK2
# ==========================================================
for i in range(N-1):
th, om = theta_rk2[i], omega_rk2[i]
k1_t, k1_w = f(th, om)
k2_t, k2_w = f(th+0.5*dt*k1_t, om+0.5*dt*k1_w)
theta_rk2[i+1] = th + dt*k2_t
omega_rk2[i+1] = om + dt*k2_w
# ==========================================================
# EULER
# ==========================================================
for i in range(N-1):
dth, dom = f(theta_euler[i], omega_euler[i])
theta_euler[i+1] = theta_euler[i] + dt*dth
omega_euler[i+1] = omega_euler[i] + dt*dom
# ==========================================================
# PEAKS (post-process)
# ==========================================================
idx_rk2, _ = find_peaks(theta_rk2, distance=100)
idx_eu, _ = find_peaks(theta_euler, distance=100)
peaks_rk2 = theta_rk2[idx_rk2]
peaks_eu = theta_euler[idx_eu]
times_rk2 = t[idx_rk2]
times_eu = t[idx_eu]
# ==========================================================
# Z SCORE
# ==========================================================
def z(x):
if len(x) < 2:
return np.zeros(len(x))
return (x - np.mean(x)) / np.std(x)
z_rk2 = z(peaks_rk2)
z_eu = z(peaks_eu)
# ==========================================================
# FIGURE
# ==========================================================
fig = plt.figure(figsize=(10,6), dpi = 60)
gs = fig.add_gridspec(2,2)
ax_p = fig.add_subplot(gs[:,0])
ax_t = fig.add_subplot(gs[:,1])
# ax_z = fig.add_subplot(gs[1,1])
# pendule
line_rk4, = ax_p.plot([], [], lw=2, label="RK4")
line_rk2, = ax_p.plot([], [], "--", label="RK2")
line_eu, = ax_p.plot([], [], ":", label="Euler")
mass_rk4, = ax_p.plot([], [], "o")
mass_rk2, = ax_p.plot([], [], "o")
mass_eu, = ax_p.plot([], [], "o")
ax_p.set_xlim(-1.2,1.2)
ax_p.set_ylim(-1.2,1.2)
ax_p.set_aspect("equal")
ax_p.grid()
ax_p.legend()
ax_p.set_title("Pendule")
# theta
line_t_rk4, = ax_t.plot([], [], label="RK4")
line_t_rk2, = ax_t.plot([], [], label="RK2")
line_t_eu, = ax_t.plot([], [], label="Euler")
ax_t.set_xlim(0,T)
ax_t.set_ylim(-2,2)
ax_t.grid()
ax_t.legend()
ax_t.set_title("Angle")
# z
# line_z_rk2, = ax_z.plot([], [], "o-", label="RK2")
# line_z_eu, = ax_z.plot([], [], "o-", label="Euler")
# ax_z.grid()
# ax_z.legend()
# ax_z.set_title("Z-score")
# ==========================================================
# VIDEO WRITER
# ==========================================================
Writer = animation.FFMpegWriter
writer = Writer(fps=30, metadata=dict(artist='ChatGPT'), bitrate=1800)
# ==========================================================
# UPDATE
# ==========================================================
def update(i):
# pendule coords
xrk4 = L*np.sin(theta_rk4[i])
yrk4 = -L*np.cos(theta_rk4[i])
xrk2 = L*np.sin(theta_rk2[i])
yrk2 = -L*np.cos(theta_rk2[i])
xeu = L*np.sin(theta_euler[i])
yeu = -L*np.cos(theta_euler[i])
line_rk4.set_data([0,xrk4],[0,yrk4])
line_rk2.set_data([0,xrk2],[0,yrk2])
line_eu.set_data([0,xeu],[0,yeu])
mass_rk4.set_data([xrk4],[yrk4])
mass_rk2.set_data([xrk2],[yrk2])
mass_eu.set_data([xeu],[yeu])
# theta curves
line_t_rk4.set_data(t[:i], theta_rk4[:i])
line_t_rk2.set_data(t[:i], theta_rk2[:i])
line_t_eu.set_data(t[:i], theta_euler[:i])
## # z-score progressive
## nrk2 = np.sum(times_rk2 <= t[i])
## neu = np.sum(times_eu <= t[i])
##
## if nrk2 > 1:
## line_z_rk2.set_data(range(nrk2), z_rk2[:nrk2])
##
## if neu > 1:
## line_z_eu.set_data(range(neu), z_eu[:neu])
##
## return (line_rk4,line_rk2,line_eu,
## mass_rk4,mass_rk2,mass_eu,
## line_t_rk4,line_t_rk2,line_t_eu,
## line_z_rk2,line_z_eu)
# ==========================================================
# SAVE ANIMATION
# ==========================================================
ani = animation.FuncAnimation(fig, update, frames=N, interval=20)
# CODE POUR SAUVEGARDER DES GIF PAS TROP VOLUMINEUX
# METTRE LE dpi = 60 DANS LES FIGURES
# from matplotlib.animation import PillowWriter
# ani = FuncAnimation(fig, update, frames=800, interval=20)
# ani.save("Pendule_Euler_RK4.gif",
# writer=PillowWriter(fps=20)) #, save_count=80)
# ani.save("pendule_erreur.mp4", writer=writer)
plt.show()Comparaison Euler RK4 (avec z-score)
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.signal import find_peaks
# ==========================================================
# PARAMETRES PHYSIQUES
# ==========================================================
g = 9.81
L = 1.0
# ==========================================================
# PARAMETRES NUMERIQUES
# ==========================================================
dt = 0.01
T = 35
N = int(T / dt)
# ==========================================================
# CONDITIONS INITIALES
# ==========================================================
theta0 = np.pi / 4
omega0 = 0.0
# ==========================================================
# TABLEAUX
# ==========================================================
t = np.zeros(N)
theta_rk4 = np.zeros(N)
omega_rk4 = np.zeros(N)
theta_rk2 = np.zeros(N)
omega_rk2 = np.zeros(N)
theta_euler = np.zeros(N)
omega_euler = np.zeros(N)
# Initialisation
theta_rk4[0] = theta0
omega_rk4[0] = omega0
theta_rk2[0] = theta0
omega_rk2[0] = omega0
theta_euler[0] = theta0
omega_euler[0] = omega0
# ==========================================================
# SYSTEME DIFFERENTIEL
# ==========================================================
def f(theta, omega):
dtheta = omega
domega = -(g/L) * np.sin(theta)
return dtheta, domega
# ==========================================================
# SIMULATION RK4 (REFERENCE)
# ==========================================================
for n in range(N - 1):
theta = theta_rk4[n]
omega = omega_rk4[n]
k1_t, k1_w = f(theta, omega)
k2_t, k2_w = f(theta + 0.5 * dt * k1_t,
omega + 0.5 * dt * k1_w)
k3_t, k3_w = f(theta + 0.5 * dt * k2_t,
omega + 0.5 * dt * k2_w)
k4_t, k4_w = f(theta + dt * k3_t,
omega + dt * k3_w)
theta_rk4[n+1] = theta + (dt/6)*(k1_t + 2*k2_t + 2*k3_t + k4_t)
omega_rk4[n+1] = omega + (dt/6)*(k1_w + 2*k2_w + 2*k3_w + k4_w)
t[n+1] = t[n] + dt
# ==========================================================
# SIMULATION RK2
# ==========================================================
for n in range(N - 1):
theta = theta_rk2[n]
omega = omega_rk2[n]
k1_t, k1_w = f(theta, omega)
k2_t, k2_w = f(theta + 0.5 * dt * k1_t,
omega + 0.5 * dt * k1_w)
theta_rk2[n+1] = theta + dt * k2_t
omega_rk2[n+1] = omega + dt * k2_w
# ==========================================================
# SIMULATION EULER
# ==========================================================
for n in range(N - 1):
theta = theta_euler[n]
omega = omega_euler[n]
dtheta, domega = f(theta, omega)
theta_euler[n+1] = theta + dt * dtheta
omega_euler[n+1] = omega + dt * domega
# ==========================================================
# DETECTION DES MAXIMA
# ==========================================================
##def find_peaks(theta, t,
## min_peak_height=0.1,
## min_distance=100):
##
## peaks_theta = []
## peaks_time = []
## peaks_index = []
##
## for i in range(1, len(theta)-1):
##
## if (theta[i] > theta[i-1]
## and theta[i] > theta[i+1]
## and theta[i] > min_peak_height
## and (len(peaks_index) == 0
## or i - peaks_index[-1] > min_distance)):
##
## peaks_theta.append(theta[i])
## peaks_time.append(t[i])
## peaks_index.append(i)
##
## return peaks_theta, peaks_time
##
##
### Maxima RK2
##peaks_rk2, times_rk2 = find_peaks(theta_rk2, t)
##
### Maxima Euler
##peaks_euler, times_euler = find_peaks(theta_euler, t)
peaks_idx, _ = find_peaks(theta_euler,
height=0.1,
distance=100)
# ==========================================================
# DETECTION DES MAXIMA AVEC SCIPY
# ==========================================================
# RK2
idx_rk2, _ = find_peaks(theta_rk2,
height=0.1,
distance=100)
peaks_rk2 = theta_rk2[idx_rk2]
times_rk2 = t[idx_rk2]
# Euler
idx_euler, _ = find_peaks(theta_euler,
height=0.1,
distance=100)
peaks_euler = theta_euler[idx_euler]
times_euler = t[idx_euler]
# ==========================================================
# CALCUL DES Z-SCORES
# ==========================================================
zscore_rk2 = stats.zscore(peaks_rk2)
zscore_euler = stats.zscore(peaks_euler)
# ==========================================================
# FIGURE : 2 SUBPLOTS
# ==========================================================
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7, 6))
# ==========================================================
# SUBPLOT 1 : SIMULATIONS + BARRES θmax
# ==========================================================
ax1.plot(t, theta_rk4,
label='RK4 (référence)',
linewidth=2)
ax1.plot(t, theta_rk2,
label='RK2',
linestyle='--',
linewidth=2)
ax1.plot(t, theta_euler,
label='Euler',
linestyle=':',
linewidth=2)
# Barres verticales RK2
for time_peak, peak in zip(times_rk2, peaks_rk2):
ax1.vlines(time_peak,
0,
peak,
color='blue',
alpha=0.4)
# Barres verticales Euler
for time_peak, peak in zip(times_euler, peaks_euler):
ax1.vlines(time_peak,
0,
peak,
color='orange',
alpha=0.4)
# Points maxima
ax1.plot(times_rk2, peaks_rk2,
marker = 'o', markersize=5)
ax1.plot(times_euler, peaks_euler,
marker = 'o', color='orange',
markersize=5)
ax1.set_xlabel('Temps (s)')
ax1.set_ylabel('θ (rad)')
ax1.set_title('Simulation temporelle et θ_max')
ax1.legend()
ax1.grid(True)
# ==========================================================
# SUBPLOT 2 : Z-SCORES
# ==========================================================
ax2.plot(zscore_rk2,
'o-',
label='RK2',
color='blue',
markersize=5)
ax2.plot(zscore_euler,
'o-',
label='Euler',
color='orange',
markersize=5)
ax2.axhline(0,
color='gray',
linestyle='--',
linewidth=1)
ax2.set_xlabel('Nombre de périodes')
ax2.set_ylabel('Z-score')
ax2.set_title('Z-score des θ_max')
ax2.legend()
ax2.grid(True)
plt.tight_layout()
plt.rcParams.update({'font.size': 12, 'axes.titlesize': 7, 'legend.fontsize': 6})
plt.savefig("Comparaison_Euler_RK2.svg", format="svg", bbox_inches='tight')
plt.savefig("Comparaison_Euler_RK2.eps", format="eps", bbox_inches='tight')
plt.show()Regression Lineaire avec et sans incertitude
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
# ============================
# A. DÉFINITION DU MODÈLE
# ============================
def modele_lineaire(x, params):
"""Modèle linéaire : y = a x + b"""
a, b = params
return a * x + b
# ============================
# B. GÉNÉRATION DES DONNÉES
# ============================
rng = np.random.default_rng(0)
# vrais paramètres
a_vrai = 2.0
b_vrai = 5.0
x_obs = np.linspace(0, 10, 30)
y_theo = modele_lineaire(x_obs, [a_vrai, b_vrai])
# incertitudes expérimentales (hétérogènes pour montrer l'effet)
sigma = 0.5 + 0.3 * x_obs # incertitudes croissantes
# données bruitées
y_obs = y_theo + rng.normal(0, sigma)
# ============================
# C. ESTIMATION DES PARAMÈTRES
# ============================
# --- 1) Régression linéaire simple (OLS) ---
def chi2_OLS(params):
y_mod = modele_lineaire(x_obs, params)
return np.sum((y_obs - y_mod)**2)
params_init = [1.0, 1.0]
res_OLS = minimize(chi2_OLS, params_init)
params_OLS = res_OLS.x
a_OLS, b_OLS = params_OLS
# --- 2) Régression linéaire pondérée (WLS / χ²) ---
def chi2_WLS(params):
y_mod = modele_lineaire(x_obs, params)
return np.sum(((y_obs - y_mod) / sigma)**2)
res_WLS = minimize(chi2_WLS, params_init)
params_WLS = res_WLS.x
a_WLS, b_WLS = params_WLS
print("=== Paramètres vrais ===")
print(f"a = {a_vrai:.3f}, b = {b_vrai:.3f}\n")
print("=== Régression OLS (non pondérée) ===")
print(f"a = {a_OLS:.3f}, b = {b_OLS:.3f}")
print("\n=== Régression WLS (pondérée par 1/sigma^2) ===")
print(f"a = {a_WLS:.3f}, b = {b_WLS:.3f}")
# ============================
# D. VISUALISATION
# ============================
x_fine = np.linspace(0, 10, 300)
plt.figure(figsize=(8,6))
# Données expérimentales
plt.errorbar(x_obs, y_obs, yerr=sigma, fmt='o', label="Données expérimentales")
# Modèle vrai
plt.plot(x_fine, modele_lineaire(x_fine, [a_vrai, b_vrai]),
'k--', label="Modèle vrai")
# Fit OLS
plt.plot(x_fine, modele_lineaire(x_fine, params_OLS),
'b-', label="Régression OLS (non pondérée)")
# Fit WLS
plt.plot(x_fine, modele_lineaire(x_fine, params_WLS),
'r-', label="Régression WLS (pondérée)")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Régression linéaire : OLS vs WLS (avec incertitudes)")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()
Modele decroissance exponentielle
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
# ============================
# A. DÉFINITION DU MODÈLE
# ============================
def modele_exponentiel(t, params):
"""
Modèle physique : décroissance exponentielle
N(t) = N0 * exp(-t / tau)
params = [N0, tau]
"""
N0, tau = params
return N0 * np.exp(-t / tau)
# ============================
# B. GÉNÉRATION / CHARGEMENT DES DONNÉES
# ============================
# Ici on simule des données expérimentales pour l'exemple
rng = np.random.default_rng(0)
# "vrais" paramètres physiques
N0_vrai = 100.0
tau_vrai = 2.5
t_obs = np.linspace(0, 10, 30) # temps mesurés
N_theo = modele_exponentiel(t_obs, [N0_vrai, tau_vrai])
# incertitudes expérimentales (par ex. bruit de comptage)
sigma = 5.0 * np.ones_like(t_obs)
# données observées = théorie + bruit gaussien
N_obs = N_theo + rng.normal(0, sigma)
# ============================
# C. ESTIMATION DES PARAMÈTRES
# ============================
def chi2(params):
N_mod = modele_exponentiel(t_obs, params)
return np.sum(((N_obs - N_mod) / sigma)**2)
# point de départ pour l'optimisation
params_init = [80.0, 1.0]
result = minimize(chi2, params_init)
params_fit = result.x
N0_fit, tau_fit = params_fit
print("Paramètres vrais : N0 = {:.2f}, tau = {:.2f}".format(N0_vrai, tau_vrai))
print("Paramètres ajustés : N0 = {:.2f}, tau = {:.2f}".format(N0_fit, tau_fit))
# ============================
# D. VALIDATION / DIAGNOSTICS
# ============================
N_mod_fit = modele_exponentiel(t_obs, params_fit)
residus = N_obs - N_mod_fit
z_scores = residus / sigma
chi2_val = chi2(params_fit)
ddl = len(t_obs) - len(params_fit) # degrés de liberté
chi2_reduit = chi2_val / ddl
print("χ² = {:.2f}, ddl = {}, χ² réduit = {:.2f}".format(chi2_val, ddl, chi2_reduit))
# ============================
# E. VISUALISATION
# ============================
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7,8), sharex=True)
# Données + modèle
ax1.errorbar(t_obs, N_obs, yerr=sigma, fmt='o', label="Données expérimentales")
t_fine = np.linspace(0, 10, 300)
# ax1.plot(t_fine, modele_exponentiel(t_fine, [N0_vrai, tau_vrai]),'k--', label="Modèle vrai")
ax1.plot(t_fine, modele_exponentiel(t_fine, params_fit), 'r-', label="Modèle ajusté")
ax1.set_ylabel("N(t)")
ax1.set_title("Décroissance exponentielle : données vs modèle")
ax1.legend()
ax1.grid(True)
# Résidus normalisés (z-scores)
ax2.axhline(0, color='k', linewidth=1)
# ax2.errorbar(t_obs, z_scores, yerr=np.ones_like(z_scores), fmt='o', label="z-scores")
ax2.errorbar(t_obs, z_scores, fmt='o', label="z-scores")
ax2.set_xlabel("Temps t")
ax2.set_ylabel("z-score")
ax2.set_title("Résidus normalisés (z-scores)")
ax2.grid(True)
plt.tight_layout()
plt.show()
Battements 435Hz – 440 Hz
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# Paramètres du signal
frequence1 = 435 # Hz
frequence2 = 440 # Hz
amplitude = 1
duree = 10 # secondes (20 ms)
fe = 44100 # fréquence d'échantillonnage (Hz)
tau = 0.3
# Axe temporel
t = np.linspace(0, duree, int(fe * duree), endpoint=False)
# Génération du signal sinusoïdal
audio_435Hz = amplitude * np.sin(2 * np.pi * frequence1 * t)
audio_440Hz = amplitude * np.sin(2 * np.pi * frequence2 * t)
audio_Battements = audio_435Hz + audio_440Hz
audio_Battements_Decay = audio_Battements * np.exp(-t/tau)
# FFT
N = len(audio_Battements)
fft_values = np.fft.fft(audio_Battements)
fft_freqs = np.fft.fftfreq(N, 1 / fe)
# On garde uniquement les fréquences positives
mask = fft_freqs >= 0
fft_freqs = fft_freqs[mask]
fft_values = fft_values[mask]
# Amplitude (normalisée)
amplitude_spectre = np.abs(fft_values) / N
# Génération du spectrogramme
frequences, temps, Sxx = signal.spectrogram(
audio_Battements,
fe,
nperseg=4096
)
# Affichage
fig, axes = plt.subplots(2, 1, figsize=(15, 10))
##plt.title("Signal sinusoïdal de 435 Hz")
##plt.xlabel("Temps (s)")
##plt.ylabel("Amplitude")
##plt.grid(True)
##
##plt.tight_layout()
##plt.show()
# LA FONCTION title, ... S'APPLIQUE A L'OBJET plt. SI L'oN VEUT UTILISER LES OBJETS GENERES PAR LE SUBPLOT (ax[], ...) IL FAUT UTILISER LES FONCTIONS set_title(), ...
axes[0].plot(t[:25000], audio_Battements[:25000], color='green', Label = 'Sans attenuation')
axes[0].plot(t[:25000], audio_Battements_Decay[:25000], color='blue', Label = 'Avec attenuation')
axes[0].set_title("Battements (435 Hz - 440 Hz)")
axes[0].set_xlabel("Temps (s)")
axes[0].set_ylabel("Amplitude")
axes[0].legend(loc="upper left")
axes[0].grid(True)
# axes[1].plot(fft_freqs, amplitude_spectre, color='red')
axes[1].loglog(fft_freqs,amplitude_spectre, color='red')
axes[1].set_title("Spectre de fréquence (FFT)")
axes[1].set_xlabel("Fréquence (Hz)")
axes[1].set_ylabel("Amplitude")
# axes[1].set_xlim(0, 1000) # zoom utile
axes[1].grid(True)
plt.rcParams.update({'font.size': 12, 'axes.titlesize': 14, 'legend.fontsize': 10})
plt.savefig("Battements_435Hz_440Hz.svg", format="svg", bbox_inches='tight')
plt.savefig("Battements_435Hz_440Hz.eps", format="eps", bbox_inches='tight')
plt.tight_layout()
# Spectrogramme
fig, ax2 = plt.subplots()
ax2.set_title("Spectrogramme")
ax2.set_xlabel("Temps (s)")
ax2.set_ylabel("Fréquence (Hz)")
ax2.set_ylim(0, 2000)
spectre = ax2.pcolormesh( temps, frequences, 10 * np.log10(Sxx), shading='gouraud', cmap='viridis' )
# Barre de couleur
fig.colorbar(
spectre,
ax=ax2,
label="Puissance (dB)"
)
# Ajustement automatique
plt.tight_layout()
# Affichage
plt.show()
Tube de Kundt
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter
# --- Paramètres ---
L = 10.0 # Longueur de la corde (m)
A = 1.0 # Amplitude (m)
f = 400 # Fréquence (Hz)
v = 340 # Vitesse de l'onde (m/s)
modes = [1, 2, 3, 4] # Modes harmoniques à afficher
# --- Création de la figure ---
fig, axes = plt.subplots(len(modes), 1, figsize=(12, 8), sharex=True, dpi = 60)
if len(modes) == 1:
axes = [axes]
# --- Grille spatiale et temporelle ---
x = np.linspace(0, L, 1000)
t = np.linspace(0, 2, 100)
# --- Liste pour stocker les animations ---
anis = []
# --- Animation pour chaque mode ---
for i, n in enumerate(modes):
ax = axes[i]
k = n * np.pi / L
omega = 2 * np.pi * f * n
# Tracer les nœuds et ventres
nodes = np.linspace(0, L, n+1)
for node in nodes:
ax.axvline(x=node, color='red', linestyle='--', alpha=0.5, label='Nœud' if i==0 else "")
antinodes = np.linspace(L/(2*n), L, n)
for antinode in antinodes:
ax.axvline(x=antinode, color='green', linestyle=':', alpha=0.5, label='Ventre' if i==0 else "")
# Initialisation
line, = ax.plot([], [], lw=2, color='blue')
ax.set_ylim(-1.5*A, 1.5*A)
ax.set_xlim(0, L)
ax.set_ylabel(f'Mode {n}')
ax.legend(loc='upper right')
ax.grid(True)
# Fonction d'initialisation
def init():
line.set_data([], [])
return line,
# Fonction d'animation (capture de n, k, omega, x, line)
def update(frame, n=n, k=k, omega=omega, x=x, line=line):
y = A * np.sin(k * x) * np.cos(omega * frame)
line.set_data(x, y)
return line,
# Animation (blit=False pour éviter les conflits)
ani = FuncAnimation(
fig, update, frames=t,
init_func=init, blit=False, interval=50, repeat=True
)
# anis.append(ani) # Stocker l'animation
# CODE POUR SAUVEGARDER DES GIF PAS TROP VOLUMINEUX
# METTRE LE dpi = 60 DANS LES FIGURES
# from matplotlib.animation import PillowWriter
# ani = FuncAnimation(fig, update, frames=80, interval=20)
# ani.save("Tube_Kundt.gif",
# writer=PillowWriter(fps=20)) #, save_count=80)
anis.append(ani) # Stocker l'animation
# --- Affichage ---
plt.suptitle("Ondes stationnaires : Nœuds (rouge) et Ventres (vert)", y=1.02)
plt.xlabel('Position (m)')
plt.tight_layout()
plt.show()Explications Ondes Progressives Battements
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
import tkinter as tk
from tkinter import ttk
# Paramètres par défaut
L = 50.0 # Longueur du tuyau
N = 1000
x = np.linspace(0, L, N)
# Fonction pour générer une onde plane progressive
def onde_plane_progressive(x, A, k, omega, c, phi, x0, t):
return A * np.cos(k * (x - x0) - omega * t + phi)
def onde_plane_regressive(x, A, k, omega, c, phi, x0, t):
return A * np.cos(-k * (x - L + x0) - omega * t - phi)
# Classe pour l'application Tkinter
class OndeApp:
def __init__(self, root):
self.root = root
self.root.title("Simulation d'ondes planes dans un tuyau")
# Variables de contrôle
self.A1 = tk.DoubleVar(value=1.0)
self.A2 = tk.DoubleVar(value=1.0)
self.k1 = tk.DoubleVar(value=2*np.pi*100/L)
self.k2 = tk.DoubleVar(value=2*np.pi*100/L)
self.omega1 = tk.DoubleVar(value=1.0)
self.omega2 = tk.DoubleVar(value=1.0)
self.c1 = tk.DoubleVar(value=1.0)
self.c2 = tk.DoubleVar(value=1.0)
self.phi1 = tk.DoubleVar(value=0.0)
self.phi2 = tk.DoubleVar(value=0.0)
self.x01 = tk.DoubleVar(value=0.0)
self.x02 = tk.DoubleVar(value=0.0)
self.L = tk.DoubleVar(value=L)
# Création de l'interface
self.create_widgets()
self.create_plot()
def create_widgets(self):
# Frame pour les contrôles
control_frame = ttk.LabelFrame(self.root, text="Paramètres des ondes")
control_frame.pack(side=tk.LEFT, fill=tk.BOTH, expand=False, padx=5, pady=5)
# Paramètres Onde 1
ttk.Label(control_frame, text="Onde 1").grid(row=0, column=0, columnspan=2, sticky=tk.W, pady=2)
ttk.Label(control_frame, text="Amplitude:").grid(row=1, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=5.0, variable=self.A1, orient=tk.HORIZONTAL).grid(row=1, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Vecteur d'onde (k):").grid(row=2, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=20.0, variable=self.k1, orient=tk.HORIZONTAL).grid(row=2, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Pulsation (ω):").grid(row=3, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=10.0, variable=self.omega1, orient=tk.HORIZONTAL).grid(row=3, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Célérité:").grid(row=4, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=5.0, variable=self.c1, orient=tk.HORIZONTAL).grid(row=4, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Déphasage:").grid(row=5, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0, to=2*np.pi, variable=self.phi1, orient=tk.HORIZONTAL).grid(row=5, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Position départ:").grid(row=6, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0, to=self.L.get(), variable=self.x01, orient=tk.HORIZONTAL).grid(row=6, column=1, sticky=tk.EW)
# Paramètres Onde 2
ttk.Label(control_frame, text="Onde 2").grid(row=7, column=0, columnspan=2, sticky=tk.W, pady=2)
ttk.Label(control_frame, text="Amplitude:").grid(row=8, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=5.0, variable=self.A2, orient=tk.HORIZONTAL).grid(row=8, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Vecteur d'onde (k):").grid(row=9, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=20.0, variable=self.k2, orient=tk.HORIZONTAL).grid(row=9, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Pulsation (ω):").grid(row=10, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=10.0, variable=self.omega2, orient=tk.HORIZONTAL).grid(row=10, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Célérité:").grid(row=11, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=5.0, variable=self.c2, orient=tk.HORIZONTAL).grid(row=11, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Déphasage:").grid(row=12, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0, to=2*np.pi, variable=self.phi2, orient=tk.HORIZONTAL).grid(row=12, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Position départ:").grid(row=13, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0, to=self.L.get(), variable=self.x02, orient=tk.HORIZONTAL).grid(row=13, column=1, sticky=tk.EW)
# Paramètres du tuyau
ttk.Label(control_frame, text="Longueur du tuyau:").grid(row=14, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=10, to=100, variable=self.L, orient=tk.HORIZONTAL).grid(row=14, column=1, sticky=tk.EW)
# Bouton pour redémarrer l'animation
ttk.Button(control_frame, text="Redémarrer", command=self.restart_animation).grid(row=15, column=0, columnspan=2, pady=5)
def create_plot(self):
# Frame pour le graphique
plot_frame = ttk.Frame(self.root)
plot_frame.pack(side=tk.RIGHT, fill=tk.BOTH, expand=True)
# Figure Matplotlib
self.fig, self.ax = plt.subplots(figsize=(15, 8))
self.ax.set_xlim(0, self.L.get())
self.ax.set_ylim(-5, 5)
self.ax.set_xlabel("Position dans le tuyau (m)")
self.ax.set_ylabel("Amplitude")
self.ax.set_title("Propagation d'ondes planes dans un tuyau")
self.ax.grid(True)
# Canvas Tkinter
self.canvas = FigureCanvasTkAgg(self.fig, master=plot_frame)
self.canvas.draw()
self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=True)
# Animation
self.line1, = self.ax.plot([], [], label='Onde 1', color='blue', alpha=0.7)
self.line2, = self.ax.plot([], [], label='Onde 2', color='red', alpha=0.7)
self.line_total, = self.ax.plot([], [], label='Somme des ondes', color='green', linewidth=2)
self.ax.legend()
# Dessiner le tuyau
self.tuyau = plt.Rectangle((0, -5), self.L.get(), 10, linewidth=1, edgecolor='black', facecolor='none', linestyle='--')
self.ax.add_patch(self.tuyau)
# Lancer l'animation
self.ani = FuncAnimation(self.fig, self.update, frames=np.linspace(0, 20, 200), interval=50, blit=False, repeat=True)
def update(self, t):
# Mettre à jour les données
x = np.linspace(0, self.L.get(), 1000)
y1 = onde_plane_progresive(x, self.A1.get(), self.k1.get(), self.omega1.get(), self.c1.get(), self.phi1.get(), self.x01.get(), t)
y2 = onde_plane_regressive(x, self.A2.get(), self.k2.get(), self.omega2.get(), self.c2.get(), self.phi2.get(), self.x02.get(), t)
y_total = y1 + y2
# Mettre à jour les lignes
self.line1.set_data(x, y1)
self.line2.set_data(x, y2)
self.line_total.set_data(x, y_total)
# Mettre à jour les limites
self.ax.set_xlim(0, self.L.get())
self.ax.set_ylim(-max(self.A1.get(), self.A2.get())*2, max(self.A1.get(), self.A2.get())*2)
self.tuyau.set_width(self.L.get())
self.tuyau.set_height(max(self.A1.get(), self.A2.get())*4)
self.tuyau.set_xy((0, -max(self.A1.get(), self.A2.get())*2))
self.ax.set_title(f"Propagation d'ondes planes dans un tuyau (t = {t:.2f} s)")
return self.line1, self.line2, self.line_total, self.tuyau
def restart_animation(self):
self.ani.event_source.stop()
self.ani = FuncAnimation(self.fig, self.update, frames=np.linspace(0, 20, 200), interval=50, blit=False, repeat=True)
self.canvas.draw()
# Lancer l'application
if __name__ == "__main__":
root = tk.Tk()
app = OndeApp(root)
root.mainloop()Explication Ondes Progressives Ondes Regressives Battements
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
import tkinter as tk
from tkinter import ttk
# Paramètres par défaut
L = 50.0 # Longueur du tuyau
N = 1000
x = np.linspace(0, L, N)
# Fonction pour générer une onde plane progressive
def onde_plane_progressive(x, A, k, omega, c, phi, x0, t):
return A * np.cos(k * (x - x0) - omega * t + phi)
def onde_plane_regressive(x, A, k, omega, c, phi, x0, t):
return A * np.cos(-k * (x - L + x0) - omega * t - phi)
# Classe pour l'application Tkinter
class OndeApp:
def __init__(self, root):
self.root = root
self.root.title("Simulation d'ondes planes dans un tuyau")
# Variables de contrôle
self.A1 = tk.DoubleVar(value=1.0)
self.A2 = tk.DoubleVar(value=1.0)
self.k1 = tk.DoubleVar(value=2*np.pi*100/L)
self.k2 = tk.DoubleVar(value=2*np.pi*100/L)
self.omega1 = tk.DoubleVar(value=1.0)
self.omega2 = tk.DoubleVar(value=1.0)
self.c1 = tk.DoubleVar(value=1.0)
self.c2 = tk.DoubleVar(value=1.0)
self.phi1 = tk.DoubleVar(value=0.0)
self.phi2 = tk.DoubleVar(value=0.0)
self.x01 = tk.DoubleVar(value=0.0)
self.x02 = tk.DoubleVar(value=0.0)
self.L = tk.DoubleVar(value=L)
# Création de l'interface
self.create_widgets()
self.create_plot()
def create_widgets(self):
# Frame pour les contrôles
control_frame = ttk.LabelFrame(self.root, text="Paramètres des ondes")
control_frame.pack(side=tk.LEFT, fill=tk.BOTH, expand=False, padx=5, pady=5)
# Paramètres Onde 1
ttk.Label(control_frame, text="Onde 1").grid(row=0, column=0, columnspan=2, sticky=tk.W, pady=2)
ttk.Label(control_frame, text="Amplitude:").grid(row=1, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=5.0, variable=self.A1, orient=tk.HORIZONTAL).grid(row=1, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Vecteur d'onde (k):").grid(row=2, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=20.0, variable=self.k1, orient=tk.HORIZONTAL).grid(row=2, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Pulsation (ω):").grid(row=3, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=10.0, variable=self.omega1, orient=tk.HORIZONTAL).grid(row=3, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Célérité:").grid(row=4, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=5.0, variable=self.c1, orient=tk.HORIZONTAL).grid(row=4, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Déphasage:").grid(row=5, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0, to=2*np.pi, variable=self.phi1, orient=tk.HORIZONTAL).grid(row=5, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Position départ:").grid(row=6, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0, to=self.L.get(), variable=self.x01, orient=tk.HORIZONTAL).grid(row=6, column=1, sticky=tk.EW)
# Paramètres Onde 2
ttk.Label(control_frame, text="Onde 2").grid(row=7, column=0, columnspan=2, sticky=tk.W, pady=2)
ttk.Label(control_frame, text="Amplitude:").grid(row=8, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=5.0, variable=self.A2, orient=tk.HORIZONTAL).grid(row=8, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Vecteur d'onde (k):").grid(row=9, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=20.0, variable=self.k2, orient=tk.HORIZONTAL).grid(row=9, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Pulsation (ω):").grid(row=10, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=10.0, variable=self.omega2, orient=tk.HORIZONTAL).grid(row=10, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Célérité:").grid(row=11, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0.1, to=5.0, variable=self.c2, orient=tk.HORIZONTAL).grid(row=11, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Déphasage:").grid(row=12, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0, to=2*np.pi, variable=self.phi2, orient=tk.HORIZONTAL).grid(row=12, column=1, sticky=tk.EW)
ttk.Label(control_frame, text="Position départ:").grid(row=13, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=0, to=self.L.get(), variable=self.x02, orient=tk.HORIZONTAL).grid(row=13, column=1, sticky=tk.EW)
# Paramètres du tuyau
ttk.Label(control_frame, text="Longueur du tuyau:").grid(row=14, column=0, sticky=tk.W)
ttk.Scale(control_frame, from_=10, to=100, variable=self.L, orient=tk.HORIZONTAL).grid(row=14, column=1, sticky=tk.EW)
# Bouton pour redémarrer l'animation
ttk.Button(control_frame, text="Redémarrer", command=self.restart_animation).grid(row=15, column=0, columnspan=2, pady=5)
def create_plot(self):
# Frame pour le graphique
plot_frame = ttk.Frame(self.root)
plot_frame.pack(side=tk.RIGHT, fill=tk.BOTH, expand=True)
# Figure Matplotlib
self.fig, self.ax = plt.subplots(figsize=(15, 8))
self.ax.set_xlim(0, self.L.get())
self.ax.set_ylim(-5, 5)
self.ax.set_xlabel("Position dans le tuyau (m)")
self.ax.set_ylabel("Amplitude")
self.ax.set_title("Propagation d'ondes planes dans un tuyau")
self.ax.grid(True)
# Canvas Tkinter
self.canvas = FigureCanvasTkAgg(self.fig, master=plot_frame)
self.canvas.draw()
self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=True)
# Animation
self.line1, = self.ax.plot([], [], label='Onde 1', color='blue', alpha=0.7)
self.line2, = self.ax.plot([], [], label='Onde 2', color='red', alpha=0.7)
self.line_total, = self.ax.plot([], [], label='Somme des ondes', color='green', linewidth=2)
self.ax.legend()
# Dessiner le tuyau
self.tuyau = plt.Rectangle((0, -5), self.L.get(), 10, linewidth=1, edgecolor='black', facecolor='none', linestyle='--')
self.ax.add_patch(self.tuyau)
# Lancer l'animation
self.ani = FuncAnimation(self.fig, self.update, frames=np.linspace(0, 20, 200), interval=50, blit=False, repeat=True)
def update(self, t):
# Mettre à jour les données
x = np.linspace(0, self.L.get(), 1000)
y1 = onde_plane_progressive(x, self.A1.get(), self.k1.get(), self.omega1.get(), self.c1.get(), self.phi1.get(), self.x01.get(), t)
y2 = onde_plane_regressive(x, self.A2.get(), self.k2.get(), self.omega2.get(), self.c2.get(), self.phi2.get(), self.x02.get(), t)
y_total = y1 + y2
# Mettre à jour les lignes
self.line1.set_data(x, y1)
self.line2.set_data(x, y2)
self.line_total.set_data(x, y_total)
# Mettre à jour les limites
self.ax.set_xlim(0, self.L.get())
self.ax.set_ylim(-max(self.A1.get(), self.A2.get())*2, max(self.A1.get(), self.A2.get())*2)
self.tuyau.set_width(self.L.get())
self.tuyau.set_height(max(self.A1.get(), self.A2.get())*4)
self.tuyau.set_xy((0, -max(self.A1.get(), self.A2.get())*2))
self.ax.set_title(f"Propagation d'ondes planes dans un tuyau (t = {t:.2f} s)")
return self.line1, self.line2, self.line_total, self.tuyau
def restart_animation(self):
self.ani.event_source.stop()
self.ani = FuncAnimation(self.fig, self.update, frames=np.linspace(0, 20, 200), interval=50, blit=False, repeat=True)
self.canvas.draw()
# Lancer l'application
if __name__ == "__main__":
root = tk.Tk()
app = OndeApp(root)
root.mainloop()Explication Ondes Stationnaires
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.widgets import Slider
# =====================================================
# Paramètres initiaux
# =====================================================
A0 = 1.0
f0 = 2.0
c0 = 2.0
L0 = 10.0
R0 = 1.0
# =====================================================
# Figure
# =====================================================
fig, ax = plt.subplots(figsize=(12,6))
plt.subplots_adjust(left=0.1, bottom=0.35)
x = np.linspace(0, 20, 2000)
incident_line, = ax.plot([], [], lw=2,
color='blue',
label='Onde incidente')
reflected_line, = ax.plot([], [], lw=2,
color='red',
label='Onde réfléchie')
total_line, = ax.plot([], [], lw=3,
color='black',
label='Somme')
ax.legend()
ax.grid(True)
ax.set_xlim(0,20)
ax.set_ylim(-3,3)
ax.set_xlabel("Position x")
ax.set_ylabel("Amplitude")
# =====================================================
# Sliders
# =====================================================
ax_f = plt.axes([0.15,0.25,0.7,0.03])
ax_c = plt.axes([0.15,0.20,0.7,0.03])
ax_L = plt.axes([0.15,0.15,0.7,0.03])
ax_A = plt.axes([0.15,0.10,0.7,0.03])
ax_R = plt.axes([0.15,0.05,0.7,0.03])
slider_f = Slider(ax_f,"Fréquence",0.1,10,valinit=f0)
slider_c = Slider(ax_c,"Célérité",0.2,20,valinit=c0)
slider_L = Slider(ax_L,"Longueur L",1,20,valinit=L0)
slider_A = Slider(ax_A,"Amplitude",0.1,3,valinit=A0)
slider_R = Slider(ax_R,"Réflexion R",-1,1,valinit=R0)
# =====================================================
# Animation
# =====================================================
def update(frame):
t = frame*0.01
f = slider_f.val
c = slider_c.val
L = slider_L.val
A = slider_A.val
R = slider_R.val
omega = 2*np.pi*f
k = omega/c
# ---------------------------
# Onde incidente
# ---------------------------
front_incident = c*t
ui = np.zeros_like(x)
mask_incident = x <= front_incident
ui[mask_incident] = (
A*np.cos(k*x[mask_incident] - omega*t)
)
# ---------------------------
# Onde réfléchie
# ---------------------------
ur = np.zeros_like(x)
if front_incident > L:
tr = t - L/c
reflected_front = L - c*tr
mask_reflected = x >= reflected_front
ur[mask_reflected] = (
A*R*np.cos(
k*(2*L - x[mask_reflected])
- omega*t
)
)
# ---------------------------
# Somme
# ---------------------------
utot = ui + ur
incident_line.set_data(x,ui)
reflected_line.set_data(x,ur)
total_line.set_data(x,utot)
ax.set_xlim(0,max(20,L+2))
ax.set_title(
f"L={L:.1f} "
f"f={f:.2f} Hz "
f"c={c:.2f} m/s "
f"R={R:.2f}"
)
return incident_line, reflected_line, total_line
ani = FuncAnimation(
fig,
update,
interval=30,
blit=True
)
plt.show()Illustration Cohérence Laser
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# ============================================================
# PARAMÈTRES MODIFIABLES
# ============================================================
fs = 20000 # fréquence d'échantillonnage (Hz)
f_laser = 4000 # fréquence du laser (Hz)
mod_freq = 50 # modulation lente (Hz)
noise_level = 0.1 # bruit
N = 2048 # taille FFT
# Axe temporel centré entre -t_max et +t_max
t = np.arange(N) / fs
t_centered = np.fft.fftshift(t - t.mean())
# Axe fréquentiel centré
freqs = np.fft.fftshift(np.fft.fftfreq(N, 1/fs))
# ============================================================
# SIGNAL LUMINEUX (LASER MODULÉ)
# ============================================================
def generate_laser_signal(frame):
signal = (1 + 0.5*np.sin(2*np.pi*mod_freq*t + 0.01*frame)) \
* np.sin(2*np.pi*f_laser*t)
signal += noise_level * np.random.randn(len(t))
return signal
# ============================================================
# FFT + IFFT CENTRÉES
# ============================================================
def compute_fft(signal):
spectrum = np.fft.fft(signal)
return np.fft.fftshift(spectrum)
def compute_ifft(spectrum_shifted):
# remettre le spectre dans l'ordre normal
spectrum = np.fft.ifftshift(spectrum_shifted)
# IFFT
signal = np.fft.ifft(spectrum)
# recentrer le signal temporel
return np.real(np.fft.fftshift(signal))
# ============================================================
# FIGURE AVEC DEUX SOUS-FENÊTRES
# ============================================================
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7))
# Sous-fenêtre 1 : spectre centré
line_fft, = ax1.plot(freqs, np.zeros(N), color="blue")
ax1.set_title("Analyse spectrale centrée du signal lumineux (FFT)")
ax1.set_xlabel("Fréquence (Hz)")
ax1.set_ylabel("Amplitude")
ax1.set_xlim(-fs/2, fs/2)
ax1.set_ylim(0, 200)
# Sous-fenêtre 2 : signal temporel centré
line_ifft, = ax2.plot(t_centered, np.zeros_like(t_centered), color="red")
ax2.set_title("Signal temporel reconstruit (centré autour de 0)")
ax2.set_xlabel("Temps (s)")
ax2.set_ylabel("Amplitude")
ax2.set_ylim(-2, 2)
# ============================================================
# ANIMATION
# ============================================================
def update(frame):
signal = generate_laser_signal(frame)
spectrum = compute_fft(signal)
reconstructed = compute_ifft(spectrum)
line_fft.set_ydata(np.abs(spectrum))
line_ifft.set_ydata(reconstructed)
return line_fft, line_ifft
ani = FuncAnimation(fig, update, frames=1000, interval=30, blit=True)
plt.tight_layout()
plt.show()
Retroaction Positive Instabilité
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# ============================================================
# PARAMETRES
# ============================================================
# ALI
G_ALI = 1e5
f_ALI = 10 # Hz
w_ALI = 2*np.pi*f_ALI
# Filtres
K = 1
f0 = 1000 # Hz
w0 = 2*np.pi*f0
Q = 1/np.sqrt(2)
# ============================================================
# AXE FREQUENTIEL
# ============================================================
f = np.logspace(0, 7, 5000)
w = 2*np.pi*f
# ============================================================
# FONCTION UTILITAIRE
# ============================================================
def bode_from_num_den(num, den, w):
sys = signal.TransferFunction(num, den)
w, mag, phase = signal.bode(sys, w)
return mag, phase
def nyquist_from_num_den(num, den, w):
s = 1j * w
num_eval = np.polyval(num, s)
den_eval = np.polyval(den, s)
H = num_eval / den_eval
return np.real(H), np.imag(H)
# ============================================================
# 1) G(p)
# ============================================================
num_G = [G_ALI * w_ALI]
den_G = [1, w_ALI]
mag_G, phase_G = bode_from_num_den(num_G, den_G, w)
# ============================================================
# 2) FILTRES CANONIQUES
# ============================================================
filters = {}
# ------------------------------------------------------------
# Passe-bas 1er ordre
# ------------------------------------------------------------
filters["PB1"] = (
[K*w0],
[1, w0]
)
# ------------------------------------------------------------
# Passe-haut 1er ordre
# ------------------------------------------------------------
filters["PH1"] = (
[K, 0],
[1, w0]
)
# ------------------------------------------------------------
# Passe-bas 2e ordre
# ------------------------------------------------------------
filters["PB2"] = (
[K*w0**2],
[1, w0/Q, w0**2]
)
# ------------------------------------------------------------
# Passe-bande 2e ordre
# ------------------------------------------------------------
filters["PBand2"] = (
[K*w0/Q, 0],
[1, w0/Q, w0**2]
)
# ------------------------------------------------------------
# Passe-haut 2e ordre
# ------------------------------------------------------------
filters["PH2"] = (
[K, 0, 0],
[1, w0/Q, w0**2]
)
# ============================================================
# 3) T(p)=G(p)/(1+G(p)H(p))
# ============================================================
closed_loops = {}
for name, (num_H, den_H) in filters.items():
numG = np.array(num_G)
denG = np.array(den_G)
numH = np.array(num_H)
denH = np.array(den_H)
# GH
num_open = np.polymul(numG, numH)
den_open = np.polymul(denG, denH)
# T(p)=G/(1+GH)
num_T = np.polymul(numG, denH)
den_T = np.polyadd(den_open, num_open)
closed_loops[name] = (num_T, den_T)
# ============================================================
# FIGURE 1 : BODE GAIN
# ============================================================
fig_gain, axes_gain = plt.subplots(3, 1, figsize=(12, 16))
# ------------------------------------------------------------
# G(p)
# ------------------------------------------------------------
axes_gain[0].semilogx(f, mag_G, linewidth=2)
axes_gain[0].set_title("Bode Gain - G(p)")
axes_gain[0].set_ylabel("Gain (dB)")
axes_gain[0].grid(True, which="both")
# ------------------------------------------------------------
# Filtres
# ------------------------------------------------------------
for name, (num, den) in filters.items():
mag, phase = bode_from_num_den(num, den, w)
axes_gain[1].semilogx(f, mag, label=name)
axes_gain[1].set_title("Bode Gain - Filtres canoniques")
axes_gain[1].set_ylabel("Gain (dB)")
axes_gain[1].grid(True, which="both")
axes_gain[1].legend()
# ------------------------------------------------------------
# Boucles fermées
# ------------------------------------------------------------
for name, (num, den) in closed_loops.items():
mag, phase = bode_from_num_den(num, den, w)
axes_gain[2].semilogx(f, mag, label=name)
axes_gain[2].set_title("Bode Gain - Boucle fermée T(p)")
axes_gain[2].set_xlabel("Fréquence (Hz)")
axes_gain[2].set_ylabel("Gain (dB)")
axes_gain[2].grid(True, which="both")
axes_gain[2].legend()
plt.tight_layout()
# ============================================================
# FIGURE 2 : BODE PHASE
# ============================================================
fig_phase, axes_phase = plt.subplots(3, 1, figsize=(12, 16))
# ------------------------------------------------------------
# G(p)
# ------------------------------------------------------------
axes_phase[0].semilogx(f, phase_G, linewidth=2)
axes_phase[0].set_title("Bode Phase - G(p)")
axes_phase[0].set_ylabel("Phase (deg)")
axes_phase[0].grid(True, which="both")
# ------------------------------------------------------------
# Filtres
# ------------------------------------------------------------
for name, (num, den) in filters.items():
mag, phase = bode_from_num_den(num, den, w)
axes_phase[1].semilogx(f, phase, label=name)
axes_phase[1].set_title("Bode Phase - Filtres canoniques")
axes_phase[1].set_ylabel("Phase (deg)")
axes_phase[1].grid(True, which="both")
axes_phase[1].legend()
# ------------------------------------------------------------
# Boucles fermées
# ------------------------------------------------------------
for name, (num, den) in closed_loops.items():
mag, phase = bode_from_num_den(num, den, w)
axes_phase[2].semilogx(f, phase, label=name)
axes_phase[2].set_title("Bode Phase - Boucle fermée T(p)")
axes_phase[2].set_xlabel("Fréquence (Hz)")
axes_phase[2].set_ylabel("Phase (deg)")
axes_phase[2].grid(True, which="both")
axes_phase[2].legend()
plt.tight_layout()
# ============================================================
# FIGURE 3 : NYQUIST
# ============================================================
fig_nyquist, axes_nyquist = plt.subplots(3, 1, figsize=(10, 18))
# ------------------------------------------------------------
# G(p)
# ------------------------------------------------------------
re, im = nyquist_from_num_den(num_G, den_G, w)
axes_nyquist[0].plot(re, im, linewidth=2)
axes_nyquist[0].plot(re, -im, linestyle='--')
axes_nyquist[0].plot(-1, 0, 'rx', markersize=10)
axes_nyquist[0].set_title("Nyquist - G(p)")
axes_nyquist[0].set_xlabel("Partie réelle")
axes_nyquist[0].set_ylabel("Partie imaginaire")
axes_nyquist[0].grid(True)
axes_nyquist[0].axis('equal')
# ------------------------------------------------------------
# Filtres canoniques
# ------------------------------------------------------------
for name, (num, den) in filters.items():
re, im = nyquist_from_num_den(num, den, w)
axes_nyquist[1].plot(re, im, label=name)
axes_nyquist[1].plot(re, -im, linestyle='--')
axes_nyquist[1].plot(-1, 0, 'rx', markersize=10)
axes_nyquist[1].set_title("Nyquist - Filtres canoniques")
axes_nyquist[1].set_xlabel("Partie réelle")
axes_nyquist[1].set_ylabel("Partie imaginaire")
axes_nyquist[1].grid(True)
axes_nyquist[1].legend()
axes_nyquist[1].axis('equal')
# ------------------------------------------------------------
# Boucles fermées T(p)
# ------------------------------------------------------------
for name, (num, den) in closed_loops.items():
re, im = nyquist_from_num_den(num, den, w)
axes_nyquist[2].plot(re, im, label=name)
axes_nyquist[2].plot(re, -im, linestyle='--')
axes_nyquist[2].plot(-1, 0, 'rx', markersize=10)
axes_nyquist[2].set_title("Nyquist - Boucle fermée T(p)")
axes_nyquist[2].set_xlabel("Partie réelle")
axes_nyquist[2].set_ylabel("Partie imaginaire")
axes_nyquist[2].grid(True)
axes_nyquist[2].legend()
axes_nyquist[2].axis('equal')
plt.tight_layout()
# ============================================================
# AFFICHAGE
# ============================================================
plt.show()Portraits de phase v1
# Créé par ennio2, le 29/04/2026 en Python 3.7
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
## Reduction de l ordre pour utiliser la fonction solve_ivp de SciPy qui sait résoudre uniquement des équations du premier ordre.
## y′′=f(t,y,y′), on pose v=y′ Donc :
## {y′=v, v′=y′′ => {y′=v, v′=Acos(ωt)−av−ω0^2)y}
## Pour des ordres supérieurs comme pour une equation differentielle d’ordre n, on introduire n variables pour obtenir un système de n équations du premier ordre. Exemple :
## y′′′⇒3 équations
## Méthode générale de création d'une fonction deréduction d'ordre
#def reduction_d_ordre(t, Y):
# # y, v = Y (commande de déballage de liste)
# Y = [y, v]
# dydt = v
# dvdt = A*np.cos(w*t) - a*v - w0**2*y
# return [dydt, dvdt]
# =========================
# DIFFERENTS TYPES D'EQUATIONS
# =========================
# 1. Oscillateur linéaire non amorti non forcé
def free(t, Y, w0):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = - w0**2 * (y)
return [yp, ypp]
### 11. Oscillateur RLC à Resistance compensée par une Résistance négative (souvent plutot modélisé par une équation de Van Der Pol comme la résistane négative ne peut pas compenser exactement la résistance du circuit RLC)
##def negative_resistance(t, Y, w0, a):
## y, yp = Y
## ypp = -w0**2 * y + a * yp # a = 0 → compensation par une résistance négative
## return [yp, ypp]
# 2. Oscillateur non linéaire (Borda) amorti
def damped_oscillator_Borda(t, Y, a, w0):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = -a * yp - w0**2 * np.sin(y)
return [yp, ypp]
# 3. Oscillateur non linéaire (Duffing simple)
def duffing(t, Y, a, w0, beta):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = -a * yp - w0**2 * y - beta * y**3
return [yp, ypp]
# 4. Oscillateur libre amorti
def damped_oscillator(t, Y, a, w0):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = -a * yp - w0**2 * y
return [yp, ypp]
# 5. Oscillateur forcé linéaire
def forced_oscillator(t, Y, a, w0, A, w):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = A * np.cos(w * t) - a * yp - w0**2 * y
return [yp, ypp]
# 6. Oscillateur auto-entretenu (contre-réaction) (modélisé par l'équation de Van Der Pol)
def Van_Der_Pol_oscillator(t, Y, w0, E0, y00):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = E0 * (1 - (y/y00)**2) * yp - w0**2 * y
return [yp, ypp]
### 61. Oscillateur de Wien
##def wien_oscillator(t, Y, w0, mu):
## y, yp = Y
## ypp = -w0**2 * y + mu * (1 - y**2) * yp
## return [yp, ypp]
# 7. Forçage non sinusoïdal
def square_forcing(t, Y, a, w0, A):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
force = A * np.sign(np.sin(t))
ypp = force - a * yp - w0**2 * y
return [yp, ypp]
# =========================
# CONFIGURATION DES CAS ET DEFINITION DES VARIABLES
# =========================
cases = [
{
"name": "Oscillateur non amorti linéaire",
"func": free,
"params": {"w0": 1.5},
"y0": [10, 0]
},
## {
## "name": "Oscillateur à résistance négative",
## "func": negative_resistance,
## "params": {"w0": 1.5, "a" : 0},
## "y0": [10, 0]
## },
{
"name": "Oscillateur amorti non Linéaire (Borda)",
"func": damped_oscillator_Borda,
"params": {"a": 0, "w0": 2.0},
"y0": [10, 0]
},
{
"name": "Duffing (non linéaire)",
"func": duffing,
"params": {"a": 0, "w0": 1.5, "beta": 0.5},
"y0": [10, 0]
},
{
"name": "Oscillateur amorti lineaire",
"func": damped_oscillator,
"params": {"a": 0.5, "w0": 2.0},
"y0": [1, 0]
},
{
"name": "Oscillateur amorti forcé",
"func": forced_oscillator,
"params": {"a": 0.3, "w0": 2.0, "A": 1.0, "w": 1.5},
"y0": [1, 0]
},
{
"name": "Oscillateur auto-entretenu (modélisation par l'équation de Van Der Pol)",
"func": Van_Der_Pol_oscillator,
"y0": [1, 0],
"params": {"w0": 2.0, "E0": 1.0, "y00" : 1.0} # y00 doit etre egal à la position en 0 definie par les conditions initiales
},
## {
## "name": "Oscillateur de Wien",
## "func": wien_oscillator,
## "y0": [1, 0],
## "params": {"w0": 2.0, "mu": 1.0}
## },
{
"name": "Forçage carré",
"func": square_forcing,
"params": {"a": 0.3, "w0": 2.0, "A": 1.0},
"y0": [0, 0]
}
]
# Temps
t_span = (0, 20)
t_eval = np.linspace(*t_span, 1000)
# =========================
# OBTENTION DES SOLUTIONS ET PLOT
# =========================
#création d'un solveur générique
# L'utilisation d'une fonction intermédiare avec lamba ou en définissant proprement une fonction est nécessaire car solve_ivp necessite en paramètres une fonction avec deux paramètres uniquement. Il faut donc adapter la fonction definie à la fonction d'entrée demadée par solve_ivp
### Solveur générique (version avec lambda)
##def solve_ode(system_func, y0, t_span, t_eval, params):
## sol = solve_ivp(
## lambda t, Y: system_func(t, Y, **params),
## t_span,
## y0,
## t_eval=t_eval
## )
## return sol
## Solveur numérique (version plus lisible)
def solve_ode(system_func, y0, t_span, t_eval, params):
def wrapped_system(t, Y):
return system_func(t, Y, **params)
return solve_ivp(wrapped_system, t_span, y0, t_eval=t_eval)
m = 1 # la masse doit etre définie pour passer aux energies
figs = []
for i, case in enumerate(cases):
w0 = (case["params"]["w0"])
sol = solve_ode(case["func"], case["y0"], t_span, t_eval, case["params"])
# une figure par case
fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(5, 1, figsize=(8, 6))
# subplot 1 : y(t)
ax1.plot(sol.t, sol.y[0], label=case["name"])
ax1.legend()
ax1.set_xlabel("Temps")
ax1.set_ylabel("y(t)")
ax1.set_title("Postion = f(t)")
ax1.grid()
# subplot 2 : y vs dy/dt
dy_dt = np.gradient(sol.y[0], sol.t)
ax2.plot(sol.y[0], dy_dt, label=case["name"], color='r', linewidth=2.0)
ax2.legend()
ax2.set_xlabel("Position")
ax2.set_ylabel("y(t)")
ax2.set_title("Portraits de phase")
ax2.grid()
# subplot 3 : Ec = f(t)
ax3.plot(sol.t, 1/2*m*(dy_dt)**2, color='m', label=case["name"])
ax3.legend()
ax3.set_xlabel("Temps")
ax3.set_ylabel("Ec")
ax3.set_title("Ec = f(t)")
ax3.grid()
# subplot 4 : Ep = f(t)
ax4.plot(sol.t, w0**2/m*(sol.y[0])**2, label=case["name"], color='g', linewidth=2.0)
ax4.legend()
ax4.set_xlabel("Temps")
ax4.set_ylabel("Ep")
ax4.set_title("Ep = f(t)")
ax4.grid()
# subplot 5 : Ep
ax5.plot(w0**2/m*(sol.y[0])**2, 1/2*m*(dy_dt)**2, label=case["name"], color='g', linewidth=2.0)
ax5.legend()
ax5.set_xlabel("Ep")
ax5.set_ylabel("Ec")
ax5.set_title("Ec = f(E(p)")
ax5.grid()
fig.suptitle(case["name"])
plt.tight_layout()
# fig.suptitle(case["name"])
figs.append(fig)
# L'affichage de toutes les figures d'un coup doit etre hors de la boucle for
plt.show()
Portraits de phase v2
# Créé par ennio2, le 29/04/2026 en Python 3.7
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
## Reduction de l ordre pour utiliser la fonction solve_ivp de SciPy qui sait résoudre uniquement des équations du premier ordre.
## y′′=f(t,y,y′), on pose v=y′ Donc :
## {y′=v, v′=y′′ => {y′=v, v′=Acos(ωt)−av−ω0^2)y}
## Pour des ordres supérieurs comme pour une equation differentielle d’ordre n, on introduire n variables pour obtenir un système de n équations du premier ordre. Exemple :
## y′′′⇒3 équations
## Méthode générale de création d'une fonction deréduction d'ordre
#def reduction_d_ordre(t, Y):
# # y, v = Y (commande de déballage de liste)
# Y = [y, v]
# dydt = v
# dvdt = A*np.cos(w*t) - a*v - w0**2*y
# return [dydt, dvdt]
# =========================
# DIFFERENTS TYPES D'EQUATIONS
# =========================
# 1. Oscillateur linéaire non amorti non forcé
def free(t, Y, w0):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = - w0**2 * (y)
return [yp, ypp]
### 11. Oscillateur RLC à Resistance compensée par une Résistance négative (souvent plutot modélisé par une équation de Van Der Pol comme la résistane négative ne peut pas compenser exactement la résistance du circuit RLC)
##def negative_resistance(t, Y, w0, a):
## y, yp = Y
## ypp = -w0**2 * y + a * yp # a = 0 → compensation par une résistance négative
## return [yp, ypp]
# 2. Oscillateur non linéaire (Borda) amorti
def damped_oscillator_Borda(t, Y, a, w0):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = -a * yp - w0**2 * np.sin(y)
return [yp, ypp]
# 3. Oscillateur non linéaire (Duffing simple)
def duffing(t, Y, a, w0, beta):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = -a * yp - w0**2 * y - beta * y**3
return [yp, ypp]
# 4. Oscillateur libre amorti
def damped_oscillator(t, Y, a, w0):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = -a * yp - w0**2 * y
return [yp, ypp]
# 5. Oscillateur forcé linéaire
def forced_oscillator(t, Y, a, w0, A, w):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = A * np.cos(w * t) - a * yp - w0**2 * y
return [yp, ypp]
# 6. Oscillateur auto-entretenu (contre-réaction) (modélisé par l'équation de Van Der Pol)
def Van_Der_Pol_oscillator(t, Y, w0, E0, y00):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
ypp = E0 * (1 - (y/y00)**2) * yp - w0**2 * y
return [yp, ypp]
### 61. Oscillateur de Wien
##def wien_oscillator(t, Y, w0, mu):
## y, yp = Y
## ypp = -w0**2 * y + mu * (1 - y**2) * yp
## return [yp, ypp]
# 7. Forçage non sinusoïdal
def square_forcing(t, Y, a, w0, A):
# y, yp = Y (commande de déballage de liste)
[y, yp] = Y
force = A * np.sign(np.sin(t))
ypp = force - a * yp - w0**2 * y
return [yp, ypp]
# =========================
# CONFIGURATION DES CAS ET DEFINITION DES VARIABLES
# =========================
cases = [
{
"name": "Oscillateur non amorti linéaire",
"func": free,
"params": {"w0": 1.5},
"y0": [10, 0]
},
## {
## "name": "Oscillateur à résistance négative",
## "func": negative_resistance,
## "params": {"w0": 1.5, "a" : 0},
## "y0": [10, 0]
## },
{
"name": "Oscillateur amorti non Linéaire (Borda)",
"func": damped_oscillator_Borda,
"params": {"a": 0, "w0": 2.0},
"y0": [10, 0]
},
{
"name": "Duffing (non linéaire)",
"func": duffing,
"params": {"a": 0, "w0": 1.5, "beta": 0.5},
"y0": [10, 0]
},
{
"name": "Oscillateur amorti lineaire",
"func": damped_oscillator,
"params": {"a": 0.5, "w0": 2.0},
"y0": [1, 0]
},
{
"name": "Oscillateur amorti forcé",
"func": forced_oscillator,
"params": {"a": 0.3, "w0": 2.0, "A": 1.0, "w": 1.5},
"y0": [1, 0]
},
{
"name": "Oscillateur auto-entretenu (modélisation par l'équation de Van Der Pol)",
"func": Van_Der_Pol_oscillator,
"y0": [1, 0],
"params": {"w0": 2.0, "E0": 1.0, "y00" : 1.0} # y00 doit etre egal à la position en 0 definie par les conditions initiales
},
## {
## "name": "Oscillateur de Wien",
## "func": wien_oscillator,
## "y0": [1, 0],
## "params": {"w0": 2.0, "mu": 1.0}
## },
{
"name": "Forçage carré",
"func": square_forcing,
"params": {"a": 0.3, "w0": 2.0, "A": 1.0},
"y0": [0, 0]
}
]
# Temps
t_span = (0, 20)
t_eval = np.linspace(*t_span, 1000)
# =========================
# OBTENTION DES SOLUTIONS ET PLOT
# =========================
#création d'un solveur générique
# L'utilisation d'une fonction intermédiare avec lamba ou en définissant proprement une fonction est nécessaire car solve_ivp necessite en paramètres une fonction avec deux paramètres uniquement. Il faut donc adapter la fonction definie à la fonction d'entrée demadée par solve_ivp
### Solveur générique (version avec lambda)
##def solve_ode(system_func, y0, t_span, t_eval, params):
## sol = solve_ivp(
## lambda t, Y: system_func(t, Y, **params),
## t_span,
## y0,
## t_eval=t_eval
## )
## return sol
## Solveur numérique (version plus lisible)
def solve_ode(system_func, y0, t_span, t_eval, params):
def wrapped_system(t, Y):
return system_func(t, Y, **params)
return solve_ivp(wrapped_system, t_span, y0, t_eval=t_eval)
m = 1 # la masse doit etre définie pour passer aux energies
figs = []
for i, case in enumerate(cases):
w0 = (case["params"]["w0"])
sol = solve_ode(case["func"], case["y0"], t_span, t_eval, case["params"])
# une figure par case
fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(5, 1, figsize=(8, 6))
# subplot 1 : y(t)
ax1.plot(sol.t, sol.y[0], label=case["name"])
ax1.legend()
ax1.set_xlabel("Temps")
ax1.set_ylabel("y(t)")
ax1.set_title("Postion = f(t)")
ax1.grid()
# subplot 2 : y vs dy/dt
dy_dt = np.gradient(sol.y[0], sol.t)
ax2.plot(sol.y[0], dy_dt, label=case["name"], color='r', linewidth=2.0)
ax2.legend()
ax2.set_xlabel("Position")
ax2.set_ylabel("y(t)")
ax2.set_title("Portraits de phase")
ax2.grid()
# subplot 3 : Ec = f(t)
ax3.plot(sol.t, 1/2*m*(dy_dt)**2, color='m', label=case["name"])
ax3.legend()
ax3.set_xlabel("Temps")
ax3.set_ylabel("Ec")
ax3.set_title("Ec = f(t)")
ax3.grid()
# subplot 4 : Ep = f(t)
ax4.plot(sol.t, w0**2/m*(sol.y[0])**2, label=case["name"], color='g', linewidth=2.0)
ax4.legend()
ax4.set_xlabel("Temps")
ax4.set_ylabel("Ep")
ax4.set_title("Ep = f(t)")
ax4.grid()
# subplot 5 : Ep
ax5.plot(w0**2/m*(sol.y[0])**2, 1/2*m*(dy_dt)**2, label=case["name"], color='g', linewidth=2.0)
ax5.legend()
ax5.set_xlabel("Ep")
ax5.set_ylabel("Ec")
ax5.set_title("Ec = f(E(p)")
ax5.grid()
fig.suptitle(case["name"])
plt.tight_layout()
# fig.suptitle(case["name"])
figs.append(fig)
# L'affichage de toutes les figures d'un coup doit etre hors de la boucle for
plt.show()