Dwuwymiarowa symulacja światła w ośrodku Lorentza — NumPy¶
Ten notebook buduje od podstaw edukacyjny solver różnic skończonych. Porównujemy dwa zsynchronizowane eksperymenty:
- falę biegnącą przez próżnię,
- tę samą falę przechodzącą przez warstwę materiału opisaną modelem Lorentza.
Celem jest pokazanie związku między lokalną odpowiedzią materiału $P(x,y,t)$ a opóźnieniem pola $E(x,y,t)$. Model jest skalarny, bezwymiarowy i jakościowy — nie jest pełnym solverem równań Maxwella.
1. Zależności i katalog wynikowy¶
Potrzebujemy tylko NumPy, Matplotlib i Pillow. Agg pozwala zapisywać grafiki również na serwerze bez środowiska graficznego.
from dataclasses import asdict, dataclass
from pathlib import Path
import csv
import json
import matplotlib
import matplotlib.animation as animation
import matplotlib.pyplot as plt
import numpy as np
import imageio_ffmpeg
from IPython.display import HTML, display
matplotlib.rcParams["animation.ffmpeg_path"] = imageio_ffmpeg.get_ffmpeg_exe()
matplotlib.rcParams["animation.bitrate"] = 500
PROJECT_ROOT = Path.cwd().parent if Path.cwd().name == "notebooks" else Path.cwd()
OUTPUT_DIR = PROJECT_ROOT / "outputs/notebook-numpy"
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
print(f"NumPy: {np.__version__}")
print(f"Wyniki: {OUTPUT_DIR.resolve()}")
NumPy: 2.5.1 Wyniki: P:\coding\case_studies\Lorentz_Light_Simulation\outputs\notebook-numpy
2. Model fizyczny¶
Pole elektryczne opisujemy uproszczonym skalarnym równaniem falowym:
$$ \frac{\partial^2 E}{\partial t^2}=c^2\nabla^2E-\beta\frac{\partial^2P}{\partial t^2}+S. $$
Zbiorową odpowiedź związanych ładunków w materiale przybliżamy tłumionym oscylatorem Lorentza:
$$ \frac{\partial^2P}{\partial t^2}+\gamma\frac{\partial P}{\partial t}+\omega_0^2P=\alpha E. $$
alpha określa siłę pobudzania materiału, gamma jego tłumienie, omega_0 częstotliwość rezonansową, a beta wpływ polaryzacji na pole. Wszystkie wielkości pozostają w jednostkach symulacji.
3. Konfiguracja eksperymentu¶
Preset quick() jest odpowiedni do pracy interaktywnej. Źródło znajduje się kilka komórek od lewego brzegu, materiał leży pośrodku, a sonda za materiałem.
@dataclass(frozen=True)
class SimulationConfig:
nx: int = 240
ny: int = 90
steps: int = 900
dx: float = 1.0
dt: float = 0.30
c: float = 1.0
source_x: int = 5
source_omega: float = 0.08
source_t0: float = 170.0
source_width: float = 115.0
source_amplitude: float = 0.08
material_start: int = 95
material_end: int = 150
omega_0: float = 0.42
gamma: float = 0.02
alpha: float = 0.08
beta: float = 0.65
frame_stride: int = 5
probe_x: int = 185
def validate(self):
if self.c * self.dt / self.dx > 1 / np.sqrt(2):
raise ValueError("Niestabilny krok czasu: naruszony warunek CFL")
ordered = 0 < self.source_x < self.material_start < self.material_end < self.probe_x < self.nx
if not ordered:
raise ValueError("Wymagana kolejność: źródło < materiał < sonda < brzeg")
@dataclass
class SimulationResult:
config: SimulationConfig
frames_e: list
frames_p: list
frame_steps: np.ndarray
probe: np.ndarray
max_abs_e: np.ndarray
has_material: bool
cfg = SimulationConfig()
cfg.validate()
courant = cfg.c * cfg.dt / cfg.dx
print(f"Liczba Couranta: {courant:.3f}; limit 2D: {1/np.sqrt(2):.3f}")
Liczba Couranta: 0.300; limit 2D: 0.707
4. Dyskretyzacja przestrzeni i stabilność¶
Laplasjan przybliżamy pięciopunktowym stencil'em:
$$ \nabla^2E_{i,j}\approx\frac{E_{i+1,j}+E_{i-1,j}+E_{i,j+1}+E_{i,j-1}-4E_{i,j}}{\Delta x^2}. $$
Dla jawnego schematu falowego 2D obowiązuje warunek CFL:
$$\frac{c\Delta t}{\Delta x}\leq\frac{1}{\sqrt{2}}.$$
W osi Y używamy warunku periodycznego. Ponieważ źródło ma identyczną wartość na całej wysokości, rozwiązanie pozostaje falą płaską i nie pojawia się dyfrakcja od krawędzi apertury.
def laplacian(field):
lap = np.zeros_like(field)
lap[1:-1, :] = (
field[2:, :] + field[:-2, :]
+ np.roll(field[1:-1, :], 1, axis=1)
+ np.roll(field[1:-1, :], -1, axis=1)
- 4.0 * field[1:-1, :]
)
return lap
5. Aktualizacja w czasie i otwarte brzegi X¶
Centralne różnice skończone prowadzą do aktualizacji polaryzacji:
$$P^{n+1}=2P^n-P^{n-1}-\gamma\Delta t(P^n-P^{n-1})-\omega_0^2\Delta t^2P^n+\alpha\Delta t^2E^n.$$
Następnie obliczamy $P_{tt}$ i aktualizujemy $E$. Na lewym i prawym brzegu stosujemy warunek absorbujący Mur pierwszego rzędu. Przybliża on domenę otwartą i dobrze działa dla naszej fali padającej prostopadle. Nie jest jednak tak ogólny jak PML.
def run_simulation(config, material=True):
config.validate()
shape = (config.nx, config.ny)
e_prev = np.zeros(shape, dtype=np.float32)
e = np.zeros(shape, dtype=np.float32)
p_prev = np.zeros(shape, dtype=np.float32)
p = np.zeros(shape, dtype=np.float32)
material_mask = np.zeros(shape, dtype=np.float32)
if material:
material_mask[config.material_start:config.material_end, :] = 1.0
source_profile = np.zeros(shape, dtype=np.float32)
source_profile[config.source_x, :] = 1.0
frames_e, frames_p, frame_steps = [], [], []
probe = np.empty(config.steps, dtype=np.float32)
max_abs_e = np.empty(config.steps, dtype=np.float32)
courant2 = (config.c * config.dt / config.dx) ** 2
mur = (config.c * config.dt - config.dx) / (config.c * config.dt + config.dx)
for step in range(config.steps):
envelope = np.exp(-((step - config.source_t0) / config.source_width) ** 2)
pulse = config.source_amplitude * np.sin(config.source_omega * step) * envelope
p_next = (
2.0 * p - p_prev
- config.gamma * config.dt * (p - p_prev)
- config.omega_0**2 * config.dt**2 * p
+ material_mask * config.alpha * config.dt**2 * e
)
p_tt_dt2 = p_next - 2.0 * p + p_prev
e_next = (
2.0 * e - e_prev
+ courant2 * laplacian(e)
- config.beta * p_tt_dt2
+ source_profile * np.float32(pulse)
)
# Otwarte brzegi Mur w osi X.
e_next[0, :] = e[1, :] + mur * (e_next[1, :] - e[0, :])
e_next[-1, :] = e[-2, :] + mur * (e_next[-2, :] - e[-1, :])
p_next[[0, -1], :] = 0.0
if not np.isfinite(e_next).all():
raise FloatingPointError(f"Niestabilność numeryczna w kroku {step}")
e_prev, e = e, e_next.astype(np.float32, copy=False)
p_prev, p = p, p_next.astype(np.float32, copy=False)
probe[step] = np.mean(e[config.probe_x, :])
max_abs_e[step] = np.max(np.abs(e))
if step % config.frame_stride == 0 or step == config.steps - 1:
frames_e.append(e.copy())
frames_p.append(p.copy())
frame_steps.append(step)
return SimulationResult(
config, frames_e, frames_p, np.asarray(frame_steps),
probe, max_abs_e, material
)
6. Dwa zsynchronizowane eksperymenty¶
Przebieg próżniowy jest próbą kontrolną. Jedyną różnicą między symulacjami jest maska materiału, dlatego przesunięcia sygnału nie można przypisać innemu źródłu ani innej siatce.
print("Przebieg kontrolny: próżnia...")
vacuum = run_simulation(cfg, material=False)
print("Przebieg z ośrodkiem Lorentza...")
medium = run_simulation(cfg, material=True)
print(f"Zapisano {len(medium.frames_e)} klatek z każdego przebiegu.")
Przebieg kontrolny: próżnia... Przebieg z ośrodkiem Lorentza... Zapisano 181 klatek z każdego przebiegu.
7. Kontrola pojedynczej klatki¶
Kolory pokazują znak chwilowej wartości pola, nie intensywność. Czerwony i niebieski odpowiadają przeciwnym fazom. Intensywność byłaby związana z $E^2$. Panel $P$ pokazuje lokalne drgania związanych ładunków, a nie drugą wiązkę światła.
index = len(medium.frames_e) * 2 // 3
limit_e = max(np.percentile(np.abs(medium.frames_e), 99.5), np.percentile(np.abs(vacuum.frames_e), 99.5))
p_crop = medium.frames_p[index][cfg.material_start:cfg.material_end, :]
limit_p = max(np.percentile(np.abs(medium.frames_p), 99.5), 1e-6)
fig = plt.figure(figsize=(10, 6.6), constrained_layout=True)
grid = fig.add_gridspec(2, 2, width_ratios=[cfg.nx, cfg.material_end-cfg.material_start])
ax_e = fig.add_subplot(grid[0, 0])
ax_p = fig.add_subplot(grid[0, 1])
ax_v = fig.add_subplot(grid[1, 0])
fig.add_subplot(grid[1, 1]).axis("off")
ax_e.imshow(medium.frames_e[index].T, origin="lower", cmap="RdBu_r", vmin=-limit_e, vmax=limit_e, aspect="auto")
ax_p.imshow(p_crop.T, origin="lower", cmap="RdBu_r", vmin=-limit_p, vmax=limit_p, extent=(cfg.material_start, cfg.material_end, 0, cfg.ny), aspect="auto")
ax_v.imshow(vacuum.frames_e[index].T, origin="lower", cmap="RdBu_r", vmin=-limit_e, vmax=limit_e, aspect="auto")
ax_e.set_title("E — ośrodek Lorentza")
ax_p.set_title("Polaryzacja P")
ax_v.set_title("E — próżnia")
for ax in (ax_e, ax_v):
for boundary in (cfg.material_start, cfg.material_end):
ax.axvline(boundary, color="black", ls="--", lw=0.7)
for ax in (ax_e, ax_p, ax_v):
ax.set(xlabel="x", ylabel="y")
fig.suptitle(f"Krok {medium.frame_steps[index]}")
plt.show()
8. Pomiar opóźnienia¶
Sonda zapisuje średnie $E$ na całej wysokości za materiałem. Maksimum wyszukujemy tylko w oknie pierwszego przejścia, ponieważ późniejsze dzwonienie oscylatora lub słabe odbicia mogłyby być większe.
Edukacyjną estymatę efektywnego współczynnika wyznaczamy z:
$$n_{eff}\approx 1+\frac{c\,\Delta t_{delay}}{L}.$$
Nie jest to skalibrowana stała materiałowa: impuls jest również deformowany i częściowo odbijany.
expected = int(cfg.source_t0 + (cfg.probe_x-cfg.source_x)*cfg.dx/(cfg.c*cfg.dt))
start = max(0, int(expected - 2*cfg.source_width))
stop = min(cfg.steps, int(expected + 4*cfg.source_width))
vacuum_peak = start + int(np.argmax(np.abs(vacuum.probe[start:stop])))
medium_peak = start + int(np.argmax(np.abs(medium.probe[start:stop])))
delay_steps = medium_peak - vacuum_peak
delay_time = delay_steps * cfg.dt
slab_length = (cfg.material_end-cfg.material_start) * cfg.dx
n_eff = 1 + cfg.c * delay_time / slab_length
print(f"Maksimum w próżni: krok {vacuum_peak}")
print(f"Maksimum z materiałem: krok {medium_peak}")
print(f"Opóźnienie: {delay_steps} kroków = {delay_time:.3f} jednostek czasu")
print(f"Edukacyjne n_eff ~ {n_eff:.3f}")
times = np.arange(cfg.steps) * cfg.dt
fig, ax = plt.subplots(figsize=(10, 4.5), constrained_layout=True)
ax.plot(times, vacuum.probe, label="próżnia")
ax.plot(times, medium.probe, label="ośrodek Lorentza")
ax.axvline(vacuum_peak*cfg.dt, color="C0", ls="--", alpha=0.6)
ax.axvline(medium_peak*cfg.dt, color="C1", ls="--", alpha=0.6)
ax.set(xlabel="czas symulacji", ylabel="średnie E", title="Sygnał sondy za warstwą")
ax.grid(alpha=0.25)
ax.legend()
fig.savefig(OUTPUT_DIR/"probe_comparison.png", dpi=160)
plt.show()
Maksimum w próżni: krok 797 Maksimum z materiałem: krok 840 Opóźnienie: 43 kroków = 12.900 jednostek czasu Edukacyjne n_eff ~ 1.235
9. Eksport danych CSV i metadanych JSON¶
CSV umożliwia niezależną analizę sygnałów, natomiast JSON zapisuje parametry potrzebne do odtworzenia eksperymentu.
with (OUTPUT_DIR/"probe_signals.csv").open("w", newline="", encoding="utf-8") as handle:
writer = csv.writer(handle)
writer.writerow(["step", "time", "vacuum_E", "material_E", "material_max_abs_E"])
for step, (vac, mat, maximum) in enumerate(zip(vacuum.probe, medium.probe, medium.max_abs_e, strict=True)):
writer.writerow([step, step*cfg.dt, float(vac), float(mat), float(maximum)])
summary = {
"backend": "numpy",
"config": asdict(cfg),
"metrics": {
"vacuum_peak_step": vacuum_peak,
"material_peak_step": medium_peak,
"delay_steps": delay_steps,
"delay_time": delay_time,
"effective_index": n_eff,
},
}
(OUTPUT_DIR/"summary.json").write_text(json.dumps(summary, indent=2), encoding="utf-8")
print("Zapisano CSV i JSON.")
Zapisano CSV i JSON.
10. Animacja osadzona w notebooku¶
Panel polaryzacji jest przycięty do materiału. Dolny panel próżniowy ma dokładnie taki sam rozmiar i skalę pola $E$ jak panel górny.
Animacja jest kodowana jako skompresowane wideo H.264 i osadzana bezpośrednio w wyniku komórki. Nie zapisujemy osobnego pliku GIF ani MP4. Dodatkowo ograniczamy liczbę klatek, żeby notebook wyeksportowany do HTML pozostał lekki.
DISPLAY_ANIMATION = True
if DISPLAY_ANIMATION:
polarization = [f[cfg.material_start:cfg.material_end, :] for f in medium.frames_p]
target_video_frames = min(64, len(medium.frames_e))
video_indices = np.linspace(
0, len(medium.frames_e) - 1, num=target_video_frames, dtype=int
)
fig = plt.figure(figsize=(8, 5.2), dpi=80, constrained_layout=True)
grid = fig.add_gridspec(2, 2, width_ratios=[cfg.nx, cfg.material_end-cfg.material_start])
ax_e, ax_p = fig.add_subplot(grid[0, 0]), fig.add_subplot(grid[0, 1])
ax_v, ax_empty = fig.add_subplot(grid[1, 0]), fig.add_subplot(grid[1, 1])
ax_empty.axis("off")
im_e = ax_e.imshow(medium.frames_e[0].T, origin="lower", cmap="RdBu_r", vmin=-limit_e, vmax=limit_e, aspect="auto")
im_p = ax_p.imshow(polarization[0].T, origin="lower", cmap="RdBu_r", vmin=-limit_p, vmax=limit_p, extent=(cfg.material_start,cfg.material_end,0,cfg.ny), aspect="auto")
im_v = ax_v.imshow(vacuum.frames_e[0].T, origin="lower", cmap="RdBu_r", vmin=-limit_e, vmax=limit_e, aspect="auto")
ax_e.set_title("E — ośrodek Lorentza")
ax_p.set_title("Polaryzacja P")
ax_v.set_title("E — próżnia")
for ax in (ax_e, ax_v):
for boundary in (cfg.material_start, cfg.material_end):
ax.axvline(boundary, color="black", ls="--", lw=0.7)
for ax in (ax_e, ax_p, ax_v):
ax.set(xlabel="x", ylabel="y")
title = fig.suptitle("")
def update(frame_index):
i = int(video_indices[frame_index])
im_e.set_data(medium.frames_e[i].T)
im_p.set_data(polarization[i].T)
im_v.set_data(vacuum.frames_e[i].T)
title.set_text(f"Symulacja Lorentza 2D — krok {medium.frame_steps[i]}")
return im_e, im_p, im_v, title
movie = animation.FuncAnimation(fig, update, frames=len(video_indices), interval=50, blit=False)
embedded_animation = HTML(movie.to_html5_video())
plt.close(fig)
display(embedded_animation)
print(f"Osadzono {len(video_indices)} z {len(medium.frames_e)} klatek jako H.264.")
Osadzono 64 z 181 klatek jako H.264.
11. Jak interpretować wynik i co zbadać dalej¶
Oczekiwane obserwacje:
- przed dotarciem fali do materiału $P=0$,
- polaryzacja oscyluje wyłącznie wewnątrz warstwy,
- przebieg z materiałem dociera do sondy później niż próba próżniowa,
- pionowe pasy nie rozszerzają się w Y, ponieważ symulujemy falę płaską,
- słabe odbicie od granicy materiału jest zjawiskiem fizycznym modelu; odbicie od końca domeny jest artefaktem numerycznym.
Najbardziej wartościowe rozszerzenia to badanie zbieżności dla kilku wartości $\Delta x$, pomiar prędkości numerycznej w próżni, rozdzielenie prędkości fazowej i grupowej oraz zastąpienie warunku Mur warstwą PML.