5 Diskretisieren: Maxwell im Daumenkino

Ein Daumenkino erzeugt Bewegung aus Standbildern: Jede Seite unterscheidet sich ein kleines bisschen von der vorigen, und beim Durchblättern läuft der Film. Eine Simulation ist ein Daumenkino, das sich selbst zeichnet: Aus Bild \(n\) berechnet sie Bild \(n+1\) — immer wieder, tausendfach.

In Kapitel 3 haben wir geprüft, ob ein fertiger Puls die Maxwell-Gleichungen erfüllt — den Film hatten wir von Hand gezeichnet, die Gleichungen durften nur nicken. Jetzt drehen wir den Spieß um: Die Gleichungen bekommen nur das erste Bild (die Felder jetzt) und müssen alle weiteren selbst zeichnen. Das klingt nach einer Kleinigkeit — schließlich sagen die Gleichungen ja gerade, wie sich die Felder ändern. Aber der erste, völlig naheliegende Bauplan scheitert spektakulär, und die Reparatur ist eine der elegantesten Ideen der numerischen Physik. Mit ihr beginnt der Simulationsmotor dieses Buchs: FDTD — Finite Differenzen im Zeitbereich.

Lernziele

Nach diesem Kapitel kannst du …

… die zwei Wirbelgleichungen in Update-Regeln umstellen („morgen = heute + Δt mal Änderungsrate”) und eine solche Regel an Zahlen von Hand durchrechnen,
… vorführen, warum der naive Gleichschritt (Euler) für Wellen explodiert — und warum kleinere Zeitschritte ihn nicht retten,
… das Leapfrog-Prinzip erklären: E zu ganzen, B zu halben Zeitschritten — und warum das eine zentrale Differenz ist,
… das Yee-Versatz-Argument wiedergeben: B sitzt räumlich zwischen den E-Punkten,
… den Zehn-Zeilen-Kern einer 1D-Feldsimulation lesen — und ihr Ergebnis gegen die Lichtgeschwindigkeit messen.

5.1 Von der Gleichung zur Update-Regel

Ausgangspunkt ist das 1D-Gleichungspaar aus Kapitel 3, diesmal nach den Zeitableitungen aufgelöst:

\[ \frac{\partial B}{\partial t} = -\frac{\partial E}{\partial x}, \qquad \frac{\partial E}{\partial t} = -c^2\,\frac{\partial B}{\partial x}. \]

So herum gelesen sind das keine Bedingungen mehr, sondern Kochrezepte: Das räumliche Gefälle des einen Feldes sagt, wie schnell sich das andere ändert. Und „wie schnell es sich ändert” heißt im Daumenkino: was im nächsten Bild steht. Mit der Nachbar-Differenz aus Kapitel 2 — diesmal in der Zeit — wird aus jeder Gleichung eine Update-Regel:

\[ E_{\text{neu}} \;=\; E_{\text{alt}} \;+\; \Delta t \cdot \left( -c^2\, \frac{\partial B}{\partial x} \right), \]

und genauso für \(B\). Eine Handrechnung an einem konkreten Punkt, mit Gitterabstand \(\Delta x = 1\) cm und Zeitschritt \(\Delta t = \Delta x / c \approx 33{,}4\) ps (warum genau dieser Zeitschritt eine ausgezeichnete Wahl ist, klärt Kapitel 8):

import numpy as np
from scipy.constants import c

DX = 0.01                 # m
DT = DX / c               # s
print(f"Zeitschritt Δt = {DT*1e12:.2f} ps\n")

E_alt = 100.0             # V/m an unserem Punkt
B_rechts, B_links = 0.2e-6, 0.5e-6        # T, die beiden B-Nachbarn

dB_dx = (B_rechts - B_links) / DX
print(f"∂B/∂x ≈ (0,2 µT − 0,5 µT) / 1 cm = {dB_dx*1e6:+.0f} µT/m")

E_neu = E_alt + DT * (-c**2 * dB_dx)
print(f"E_neu = {E_alt} + Δt·(−c²·∂B/∂x) = {E_neu:.0f} V/m")

Zeitschritt Δt = 33.36 ps

∂B/∂x ≈ (0,2 µT − 0,5 µT) / 1 cm = -30 µT/m
E_neu = 100.0 + Δt·(−c²·∂B/∂x) = 190 V/m

Zwei Nachbarwerte angeschaut, eine Multiplikation — und der Punkt weiß, was er im nächsten Bild anzeigt. Das ganze Verfahren besteht aus nichts anderem, millionenfach wiederholt. Es bleibt nur eine Designfrage offen, und sie hat es in sich: Wo und wann genau sollen \(E\) und \(B\) auf dem Gitter wohnen?

5.2 Der naheliegende Plan — und sein Scheitern

Der erste Entwurf liegt auf der Hand: beide Felder an denselben Gitterpunkten, beide zu denselben Zeiten. Pro Zeitschritt berechnen wir beide räumlichen Ableitungen (zentral, über \(2\Delta x\)) aus den alten Werten und schreiben beide Felder fort. Sauber, symmetrisch, direkt aus dem Lehrbuch der Schrittverfahren — es ist das klassische Euler-Verfahren, hier „Gleichschritt” genannt, weil beide Felder im selben Takt marschieren.

Ob dieser Plan taugt, soll ein Experiment entscheiden. Damit der Code gleich nur noch vorführt, was hier beschlossen wird, zuerst der vollständige Versuchsaufbau:

Die Bühne. Eine 10 m lange Vakuum-Strecke, zerlegt in \(N = 1001\) Gitterpunkte im Abstand \(\Delta x = 1\) cm. Der Zeitschritt ist das \(\Delta t = \Delta x/c \approx 33\) ps aus der Handrechnung oben; gerechnet werden 220 Schritte, zusammen gut 7 ns.

Die Anfangsbedingung — das erste Daumenkino-Bild: ein ruhender Gauß-Hügel im E-Feld bei \(x = 5\) m (charakteristische Breite 0,5 m) und \(B = 0\) überall. Das ist der einfachste denkbare Start: ein einziges glattes Feldgebilde, keine Quelle, kein Antrieb. Was dieser Anfangszustand physikalisch anstellt, heben wir uns für den Schluss des Kapitels auf — für den Test zählt nur, dass beide Verfahren dasselbe glatte Feld weiterrechnen müssen.

Die Randbedingung. Irgendetwas muss an den Enden der Strecke passieren, denn die äußersten Punkte haben nur einen Nachbarn. Wir wählen die einfachste Möglichkeit: Die Update-Zeilen fassen nur die inneren Punkte an (im Code am Slice [1:-1] zu erkennen), die äußersten E-Werte werden nie verändert und bleiben für immer null. Ein festgenagelter E-Wert wirkt wie eine ideale Metallwand, also wie ein Spiegel (Kapitel 6 führt das als „PEC” richtig ein und spielt damit). Hier sind die Wände reine Kulisse: In 220 Schritten kommt eine Welle höchstens 2,2 m weit, von \(x = 5\) m aus erreicht also nichts die Ränder bei 0 und 10 m.

Die Messgröße und das Erfolgskriterium. Wir verfolgen nicht die Feldform, sondern eine einzige Zahl pro Zeitschritt: die Gesamtenergie des Felds. Felder tragen Energie (Kapitel 4); ihre Dichte ist \(\tfrac{\varepsilon_0}{2}E^2\) für das elektrische und \(\tfrac{1}{2\mu_0}B^2\) für das magnetische Feld (beide in J/m³ — vertieft in Kapitel 18). Die Funktion energie() summiert beide Anteile über die ganze Strecke. Der Grund für diese Wahl: Energie ist ein unbestechlicher Schiedsrichter. Im verlustfreien Vakuum zwischen verlustfreien Spiegeln darf die Gesamtenergie sich nicht ändern — egal, was die Felder im Detail treiben, und ohne dass wir die richtige Lösung kennen müssten. Aufgetragen wird sie relativ zum Anfangswert; die Solllinie ist also eine Waagerechte bei 1, und ein Verfahren, dessen Kurve steigt, erschafft Energie aus dem Nichts — disqualifiziert.

Die Kandidaten. Der Gleichschritt tritt nicht allein an: Als Vergleichsmaßstab läuft eine zweite Variante mit, die \(E\) und \(B\) räumlich wie zeitlich versetzt anordnet (im Code leapfrog_schritt; ihr B-Array ist einen Eintrag kürzer, weil die B-Punkte zwischen den E-Punkten sitzen). Wie dieser Versatz funktioniert und warum, ist das Thema des nächsten Abschnitts — hier darf die Variante nur schon einmal zeigen, was sie kann.

Vorhersage-Punkt

Bevor du die Zelle liest: Beide Verfahren starten mit demselben Gauß-Hügel, und die Energie-Kurven beider werden über 220 Schritte aufgezeichnet (logarithmische Achse!). Die Solllinie kennst du: konstant bei 1. Was erwartest du für die beiden Kurven — und falls eine abweicht: nach oben oder unten, allmählich oder plötzlich?

# von oben: DX, DT; c (Import)   (Handrechnungs-Zelle)
from scipy.constants import epsilon_0, mu_0

N = 1001
X = np.arange(N) * DX     # np.arange(n): die Zahlen 0, 1, …, n−1


def anfangspuls():
    e = np.exp(-((X - 5.0) / 0.5)**2)      # Gauß-Puls im E-Feld
    return e, np.zeros(N - 1)              # B: zwischen den E-Punkten


def energie(e, b):
    # np.sum: summiert alle Array-Elemente zu einer Zahl auf
    return 0.5 * (epsilon_0 * np.sum(e**2) + np.sum(b**2) / mu_0) * DX


def euler_schritt(e, b):
    """Gleichschritt: beide Felder am selben Ort, beide aus ALTEN Werten."""
    e_neu, b_neu = e.copy(), b.copy()      # .copy(): echte Kopie statt Verweis
    e_neu[1:-1] = e[1:-1] - c**2 * DT * (b[2:] - b[:-2]) / (2*DX)
    b_neu[1:-1] = b[1:-1] - DT * (e[2:] - e[:-2]) / (2*DX)
    return e_neu, b_neu


def leapfrog_schritt(e, b):
    """Versetzt: B zwischen den E-Punkten, B-Update zwischen den E-Updates."""
    b -= DT / DX * np.diff(e)              # np.diff: e[i+1] − e[i], elementweise
    e[1:-1] -= c**2 * DT / DX * np.diff(b)
    return e, b


e_eu, _ = anfangspuls(); b_eu = np.zeros(N)       # Euler: B am selben Ort
e_lf, b_lf = anfangspuls()
e0 = energie(*anfangspuls())

verlauf_eu, verlauf_lf = [], []
for n in range(220):
    e_eu, b_eu = euler_schritt(e_eu, b_eu)
    e_lf, b_lf = leapfrog_schritt(e_lf, b_lf)
    verlauf_eu.append(energie(e_eu, b_eu[:-1]) / e0)
    verlauf_lf.append(energie(e_lf, b_lf) / e0)

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(6.4, 3.4))
ax.semilogy(verlauf_eu, label="Euler-Gleichschritt")
ax.semilogy(verlauf_lf, label="Leapfrog (versetzt)")
ax.set_xlabel("Zeitschritt"); ax.set_ylabel("Energie / Anfangsenergie")
ax.legend(); plt.tight_layout(); plt.show()

print(f"Energie nach 220 Schritten — Euler: {verlauf_eu[-1]:.1e},"
      f"  Leapfrog: {verlauf_lf[-1]:.6f}")

Abbildung 5.1: Gesamtenergie über 220 Zeitschritte, logarithmisch. Der Euler-Gleichschritt explodiert um mehr als 30 Größenordnungen; der Leapfrog auf versetzten Stützstellen hält die Energie konstant. Der Knick bei Schritt ~100 ist kein Zufall — Übung 5.2 berechnet ihn.

Energie nach 220 Schritten — Euler: 1.7e+35,  Leapfrog: 1.000100

Der Gleichschritt sprengt jede Skala: 35 Größenordnungen Energiegewinn aus dem Nichts. Das ist kein Programmierfehler und keine Frage der Genauigkeit — das Verfahren selbst ist für Wellen unbrauchbar.

Naheliegende Vermutung

Vermutung: „Dann eben kleinere Zeitschritte — kleine Schritte, kleine Fehler, irgendwann stimmt es.”

Warum sie naheliegt: Bei allem bisher in diesem Buch wurde der Fehler kleiner, wenn die Schrittweite schrumpfte — das war geradezu unser Wahrheitskriterium (Kapitel 2).

Was stattdessen stimmt: Der Euler-Gleichschritt verstärkt jede Schwingung bei jedem Schritt um einen Faktor größer als eins — egal wie klein \(\Delta t\) ist. Kleinere Schritte verkleinern zwar den Faktor, aber dafür braucht dieselbe Simulationszeit mehr Schritte: Das Produkt bleibt exponentielles Wachstum, nur langsamer. Der Unterschied zwischen „ungenau” und „instabil” ist der zentrale Begriff dieses Kapitels: Ungenauigkeit schrumpft mit feinerem Gitter — Instabilität nicht. Sie ist ein Konstruktionsfehler, und man behebt sie nicht mit Geduld, sondern mit einem besseren Bauplan.

5.3 Der bessere Bauplan: doppelt versetzt

Die orangefarbene Kurve im Energie-Plot gehört zu diesem besseren Bauplan, und sein Kern ist ein einziger Gedanke: Alle Differenzen sollen zentral sein — auch die in der Zeit. Der Euler-Schritt schätzt die Zeitableitung einseitig (aus „jetzt” in Richtung „gleich”), und einseitige Differenzen waren schon in Kapitel 2 die schlechtere Wahl. Zentral wäre: Werte vor und nach dem Moment der Auswertung. Aber wie soll das gehen, ohne die Zukunft zu kennen?

Mit Versatz. Man lässt die beiden Felder abwechselnd springen:

\(E\) lebt zu den ganzen Zeitschritten \(0,\ 1,\ 2,\ \dots\)
\(B\) lebt zu den halben Zeitschritten \(\tfrac{1}{2},\ \tfrac{3}{2},\ \tfrac{5}{2},\ \dots\)

Beim Update von \(B^{1/2}\) auf \(B^{3/2}\) liegt das benutzte \(E^1\) zeitlich genau in der Mitte — die Differenz ist zentral, ganz ohne Hellseherei. Und beim Update von \(E^1\) auf \(E^2\) liegt \(B^{3/2}\) wieder genau in der Mitte. Die beiden Felder überspringen einander wie beim Bockspringen — daher der Name Leapfrog.

Derselbe Trick noch einmal im Raum: Die E-Update-Regel braucht \(\partial B/\partial x\) am Ort des E-Punktes. Sitzen die B-Punkte zwischen den E-Punkten (um eine halbe Zelle versetzt), dann ist \((B_{i+1/2} - B_{i-1/2})/\Delta x\) eine zentrale Differenz direkt am E-Punkt \(i\) — mit Spannweite \(\Delta x\) statt \(2\Delta x\), also viermal genauer als der Gleichschritt-Stencil. Dieses doppelt versetzte Gitter heißt nach seinem Erfinder Yee-Gitter (Kane Yee, 1966), und es ist bis heute das Herz jedes FDTD-Programms:

Abbildung 5.2: Das Yee-Schema. Oben der Zeitversatz: E (blau) lebt zu ganzen, B (rot) zu halben Zeitschritten; jeder Pfeil-Dreierpack ist ein Update, der benutzte Querwert liegt immer zeitlich mittig. Unten der Raumversatz: B sitzt zwischen den E-Punkten, die Differenz der beiden B-Nachbarn ist eine zentrale Differenz direkt am E-Punkt.

Im Code ist der ganze Bauplan kürzer als seine Erklärung — das ist der Kern, den du in der Energie-Zelle oben schon als leapfrog_schritt benutzt hast:

b -= DT / DX * np.diff(e)                  # B springt (halber Takt)
e[1:-1] -= c**2 * DT / DX * np.diff(b)     # dann springt E (ganzer Takt)

Diese zwei Zeilen sind die dichtesten des ganzen Buchs — in ihnen stecken die NumPy-Mechanik, das versetzte Gitter und die Randbehandlung gleichzeitig. Deshalb falten wir sie einmal komplett auseinander; dieselbe Rechnung, nur mit benannten Zwischenschritten:

delta_e = np.diff(e)        # Nachbar-Differenzen von e: eine pro Lücke
b = b - DT / DX * delta_e   # B-Update: elementweise, gleiche Länge

delta_b = np.diff(b)        # Nachbar-Differenzen von b
e_innen = e[1:-1]           # die inneren E-Punkte (ohne Rand)
e[1:-1] = e_innen - c**2 * DT / DX * delta_b

Was np.diff mechanisch tut: nichts als Nachbar-Differenzen. np.diff([a, b, c, d]) ergibt [b−a, c−b, d−c] — aus vier Werten werden drei Differenzen, denn es gibt eine pro Lücke, nicht eine pro Eintrag; das Ergebnis ist immer um eins kürzer als die Eingabe. Zwei Dinge tut np.diff ausdrücklich nicht: Es teilt nicht durch \(\Delta x\) (dieser Faktor steckt bei uns im Vorfaktor DT / DX), und es ist kein zentraler Differenzenquotient — der hieße np.gradient (Kapitel 2), behielte die Länge bei und verrechnete übernächste Nachbarn.

Warum die simple Nachbar-Differenz hier trotzdem zentral ist, zeigt die Buchhaltung der Plätze. Trag einmal auf, wer wo wohnt und wozu welche Differenz gehört:

Ort:          0       Δx      2Δx      3Δx      4Δx
e:           e[0]    e[1]    e[2]     e[3]     e[4]      (N Werte)
b:               b[0]    b[1]     b[2]     b[3]          (N−1 Werte)
np.diff(b):          d[0]     d[1]     d[2]              (N−2 Werte)
gehört zu:           e[1]     e[2]     e[3]

Die Differenz d[0] = b[1] − b[0] benutzt zwei B-Werte, die symmetrisch um den Punkt e[1] liegen — je eine halbe Zelle links und rechts. Vom E-Punkt aus gesehen ist das die zentrale Differenz aus Abbildung 5.2; die Zentrierung kommt nicht von der Funktion diff, sondern von der Lage der B-Punkte. Und die Längen gehen wie von selbst auf: e hat \(N\) Einträge, das dazwischen wohnende b hat \(N-1\), dessen Differenzen-Array np.diff(b) hat \(N-2\) — genau einen Wert für jeden inneren E-Punkt.

Damit erklärt sich auch das [1:-1] (der Slice „alles außer dem ersten und letzten Eintrag”, Kapitel 2): Links steht ein Ausschnitt mit \(N-2\) Plätzen, rechts ein Ergebnis mit \(N-2\) Werten — die Zuweisung schreibt sie elementweise in die inneren Positionen e[1] bis e[N-2]. Die beiden Randpunkte rührt die Zeile bewusst nicht an, und das ist kein NumPy-Zufall, sondern Physik: e[0] hat keinen B-Nachbarn links, e[-1] keinen rechts — die zentrale Differenz ist dort gar nicht bildbar. Randpunkte brauchen eine eigene Regel (eine Randbedingung): In diesem Kapitel bleiben sie einfach null (festgenagelt = Spiegelwand), in Kapitel 6 lernen sie Wellen hinauszulassen. Merke: Maxwell regiert innen, der Rand wird extra regiert.

Das Gitter und die Mathematik greifen ineinander wie Zahnräder. Mehr ist FDTD im Kern nicht; alles Weitere (Quellen, Materialien, Ränder) sind Anbauten an diese zwei Zeilen — sie kommen in Kapitel 6.

5.4 Die erste lebende Welle

Zeit für den Moment, auf den Teil I hingearbeitet hat: Wir geben dem Daumenkino ein erstes Bild und lassen die Maxwell-Gleichungen zeichnen. Die Bühne ist dieselbe wie beim Energie-Test (10 m, \(\Delta x = 1\) cm, Spiegelränder, Gauß-Hügel bei 5 m, \(B = 0\)) — nur das Messprogramm ist neu: Wir fotografieren das E-Feld bei Schritt 0, 150 und 300, und ab Schritt 150 notieren wir zusätzlich in jedem Schritt, wo das Maximum der rechten Bildhälfte sitzt. Durch diese Ort-über-Zeit-Punkte legen wir am Ende eine Gerade: Ihre Steigung ist die gemessene Ausbreitungsgeschwindigkeit, und die vergleichen wir mit \(c\).

Vorhersage-Punkt

Bevor du die Zelle liest: Das erste Bild ist ein Gauß-Berg im E-Feld — und B ist überall null. Es gibt also keine Information darüber, in welche Richtung dieser Puls laufen „soll” (erinnere dich an Kapitel 4: Eine laufende Welle hat \(B = E/c\), hier fehlt das B komplett). Was wird passieren? Läuft der Puls nach rechts, nach links, bleibt er stehen, zerfließt er? Lege dich fest.

# von oben: anfangspuls(), leapfrog_schritt(), X, N, DT
e, b = anfangspuls()
schnappschuesse = {0: e.copy()}
orte, zeiten = [], []

for n in range(1, 351):
    e, b = leapfrog_schritt(e, b)
    if n in (150, 300):
        schnappschuesse[n] = e.copy()
    if n >= 150:                       # rechter Puls, sauber getrennt
        # np.argmax: der INDEX des größten Werts (nicht der Wert selbst)
        orte.append(X[N//2 + np.argmax(e[N//2:])])
        zeiten.append(n * DT)

fig, ax = plt.subplots(figsize=(7.0, 3.2))
for n, farbe in ((0, "gray"), (150, "tab:blue"), (300, "tab:red")):
    ax.plot(X, schnappschuesse[n], color=farbe, label=f"t = {n*DT*1e9:.1f} ns")
ax.set_xlabel("x (m)"); ax.set_ylabel("E (V/m)")
ax.legend(); plt.tight_layout(); plt.show()

# np.polyfit(x, y, 1): Gerade durch die Messpunkte; [0] ist die Steigung
tempo = np.polyfit(zeiten, orte, 1)[0]
print(f"gemessenes Tempo: {tempo:,.0f} m/s   (c = {c:,.0f} m/s,"
      f" Abweichung {abs(tempo-c)/c:.0e})")

Abbildung 5.3: Ein E-Puls ohne Magnetfeld kann sich nicht entscheiden — und nimmt beide Richtungen: Er zerfällt in zwei Pulse halber Höhe, die mit Lichtgeschwindigkeit auseinanderlaufen.

gemessenes Tempo: 299,792,458 m/s   (c = 299,792,458 m/s, Abweichung 4e-16)

Der Puls teilt sich: zwei Kopien halber Höhe, eine nach links, eine nach rechts. Niemand hat dem Programm das beigebracht — es kennt nur die zwei Update-Zeilen. Die Physik dahinter steckt im Kapitel-4-Befund: Ein nach rechts laufender Puls braucht \(B = +E/c\), ein nach links laufender \(B = -E/c\). Unser Start (\(E\) da, \(B = 0\)) ist exakt die Summe aus je einer halben Portion von beidem — und dank Superposition läuft jede Hälfte ungestört los. Wer einen Puls will, der nur nach rechts läuft, muss das passende B mitliefern: Übung 5.4.

Und das Tempo? Die Maximum-Verfolgung liefert die Lichtgeschwindigkeit auf Maschinengenauigkeit — nicht ungefähr, sondern exakt. Das ist kein Standard, sondern ein Sonderfall: Beim gewählten Zeitschritt \(\Delta t = \Delta x/c\) rückt die Welle pro Schritt exakt eine Zelle weiter, und das Schema transportiert sie fehlerfrei („magischer Zeitschritt” — was bei anderen Zeitschritten passiert und warum größere sofort explodieren, ist das Thema von Kapitel 8).

In dieser HTML-Fassung kannst du dem Daumenkino direkt zusehen — die Steuerleiste unter dem Bild spielt den Film ab, und mit den Einzelschritt-Knöpfen lässt sich der Moment der Teilung in Ruhe durchblättern:

Code der Animation (nur in der HTML-Fassung)

# von oben: anfangspuls(), leapfrog_schritt(), X, DT
# matplotlib.animation: baut aus einer Zeichenfunktion einen Film
from matplotlib import animation
# IPython.display.HTML: bettet den fertigen Player in die Seite ein
from IPython.display import HTML

e_a, b_a = anfangspuls()
filmbilder = [e_a.copy()]
for n in range(1, 301):
    e_a, b_a = leapfrog_schritt(e_a, b_a)
    if n % 5 == 0:                      # jedes 5. Bild reicht fürs Auge
        filmbilder.append(e_a.copy())

fig_a, ax_a = plt.subplots(figsize=(7.0, 3.0))
linie, = ax_a.plot(X, filmbilder[0], color="tab:blue")
ax_a.set_xlim(0, X[-1])
ax_a.set_ylim(-0.15, 1.05)
ax_a.set_xlabel("x (m)")
ax_a.set_ylabel("E (V/m)")

def zeichne(i):
    linie.set_ydata(filmbilder[i])      # nur die Kurve austauschen
    ax_a.set_title(f"Zeitschritt {5 * i}")
    return [linie]

anim = animation.FuncAnimation(fig_a, zeichne,
                               frames=len(filmbilder), interval=50)
plt.close(fig_a)              # sonst erschiene zusätzlich ein Standbild
HTML(anim.to_jshtml(default_mode="loop"))

Und jetzt du. Die folgende Zelle ist keine Abbildung, sondern ein Labor: Sie läuft in deinem Browser, und du kannst den Code direkt ändern und mit dem Run-Knopf (oder Strg+Enter) neu ausführen — nichts muss installiert werden, und kaputtmachen kannst du auch nichts (neu laden stellt alles wieder her). Drei Vorschläge zum Ausprobieren: Mach den Puls schmaler (breite = 0.1 — bleibt das Tempo gleich?), lass den Film länger laufen (schritte = 700 — was machen die festgenagelten Ränder mit den Pulshälften?), und verschiebe den Startort.

Das Kapitel-Programm programme/kap05/kap05_leapfrog_1d.py fasst alles zusammen — Energievergleich, Pulsteilung, Tempomessung — mit assert-Schranken (Euler-Energie > 10³, Leapfrog-Energie auf 1 % konstant, Tempo auf \(10^{-6}\) genau).

Merkkasten

Eine Simulation ist ein selbstzeichnendes Daumenkino: Update-Regel = Gleichung, nach der Zeitableitung aufgelöst, plus Nachbar-Differenzen.
Instabil ≠ ungenau: Der Euler-Gleichschritt verstärkt Wellen bei jedem \(\Delta t\) exponentiell — kleinere Schritte retten ihn nicht.
Leapfrog: E zu ganzen, B zu halben Zeitschritten → die Zeitdifferenz wird zentral, das Verfahren stabil.
Yee-Gitter: B wohnt zwischen den E-Punkten → die Raumdifferenz ist zentral mit Spannweite \(\Delta x\).
Der Kern sind zwei Zeilen NumPy; ein E-Puls ohne B teilt sich in zwei Hälften (\(B = \pm E/c\)), die exakt mit \(c\) laufen.

Roter Faden

Hier sind drei Fäden zusammengelaufen: die zentralen Differenzen und das „Fehler schrumpft mit Δx²”-Kriterium aus Kapitel 2, die zwei Wirbelgleichungen aus Kapitel 3 und der \(B = E/c\)-Zusammenhang aus Kapitel 4, der die Pulsteilung erklärt. Nach vorn: Kapitel 6 baut um den Zehn-Zeilen-Kern ein vollständiges Labor — Quellen, Verluste, Spiegel, Materialwechsel. Kapitel 7 zeigt, dass unsere Update-Schleife heimlich eine Matrix-Multiplikation ist. Und Kapitel 8 beantwortet die hier offen gebliebene Frage, welcher Zeitschritt erlaubt ist — und was der „magische” so magisch macht.

Übungen

Ü 5.1 (Verstehen). Die Leapfrog-Zeitdifferenz benutzt nur zwei Werte, \(B^{n+1/2}\) und \(B^{n-1/2}\) — trotzdem nennt das Kapitel sie „zentral”. Um welchen Zeitpunkt ist sie zentriert, und welcher Stelle im Update entspricht das?

Musterlösung zu Ü 5.1

Sie ist um den ganzen Zeitschritt \(n\) zentriert: \(B^{n+1/2}\) liegt eine halbe Schrittweite danach, \(B^{n-1/2}\) eine halbe davor — genau wie bei der zentralen Differenz aus Kapitel 2, nur dass die „Nachbarn” hier in der Zeit liegen. Und \(n\) ist exakt der Moment, an dem das \(E\)-Feld lebt, dessen Update diese Differenz speist. Zentral heißt also nicht „drei Punkte”, sondern „symmetrisch um den Auswertepunkt”.

Ü 5.2 (Verstehen). Im Energie-Plot bleibt die Euler-Kurve rund 100 Schritte lang scheinbar flach und schießt dann hoch. Die schnellstwachsende Störung wird pro Schritt um den Faktor \(\sqrt{2}\) verstärkt und startet als Rundungsrauschen bei etwa \(10^{-16}\). Berechne, nach wie vielen Schritten sie die Größenordnung 1 erreicht — passt das zum Knick?

Musterlösung zu Ü 5.2

Gesucht ist \(n\) mit \(10^{-16}\cdot(\sqrt{2})^n = 1\):

n_knick = 16 / np.log10(np.sqrt(2))
print(f"n = 16 / log10(√2) ≈ {n_knick:.0f} Schritte")

n = 16 / log10(√2) ≈ 106 Schritte

Etwa 106 Schritte — genau dort knickt die Kurve hoch. Die Explosion war also von Anfang an da, nur unsichtbar klein: Der glatte Puls selbst wächst kaum, aber das unvermeidliche Rundungsrauschen der Gleitkommazahlen enthält die am schnellsten wachsenden (kurzwelligsten) Störungen, und \(\sqrt{2}^{\,n}\) holt 16 Größenordnungen in ~100 Schritten auf. Instabilität kündigt sich nicht an — sie war nur noch nicht fällig.

Ü 5.3 (Verändern). Halbiere im Kapitel-Programm den Zeitschritt (\(\Delta t = \Delta x/(2c)\), also \(S = 0{,}5\)) — bei gleicher Simulationsdauer (doppelt so viele Schritte). Sage vorher: Rettet das den Euler-Gleichschritt? Und bleibt der Leapfrog-Puls noch „magisch” perfekt?

Musterlösung zu Ü 5.3

# von oben: DX, N, anfangspuls(), energie(), e0; c (Import)
DT2 = DX / (2 * c)                         # S = 0,5

e_eu, _ = anfangspuls(); b_eu = np.zeros(N)
e_lf, b_lf = anfangspuls()
for n in range(440):                       # gleiche Dauer: 2 × 220 Schritte
    e_neu, b_neu = e_eu.copy(), b_eu.copy()
    e_neu[1:-1] = e_eu[1:-1] - c**2*DT2*(b_eu[2:] - b_eu[:-2])/(2*DX)
    b_neu[1:-1] = b_eu[1:-1] - DT2*(e_eu[2:] - e_eu[:-2])/(2*DX)
    e_eu, b_eu = e_neu, b_neu
    b_lf -= DT2/DX * np.diff(e_lf)
    e_lf[1:-1] -= c**2*DT2/DX * np.diff(b_lf)

print(f"Euler   : Energie ×{energie(e_eu, b_eu[:-1])/e0:.1e}")
print(f"Leapfrog: Energie ×{energie(e_lf, b_lf)/e0:.6f},"
      f"  Amplitude rechts: {np.max(e_lf[N//2:]):.4f}")

Euler   : Energie ×4.7e+11
Leapfrog: Energie ×1.000025,  Amplitude rechts: 0.5000

Euler explodiert weiterhin (langsamer pro Schritt, aber es sind doppelt so viele — Instabilität bleibt Instabilität). Der Leapfrog bleibt stabil und energietreu, aber die Pulsform ist nicht mehr bitgenau perfekt: Die Amplitude weicht minimal von 0,5 ab, weil bei \(S < 1\) verschiedene Wellenlängen geringfügig verschieden schnell laufen (numerische Dispersion — das Thema von Kapitel 8). Stabil und brauchbar ist \(S = 0{,}5\) trotzdem.

Ü 5.4 (Übertragen). Baue einen Anfangszustand, dessen Puls nur nach rechts läuft. (Hinweis aus Kapitel 4: Eine rechtslaufende Welle hat \(B = +E/c\) — und denke daran, dass die B-Punkte eine halbe Zelle versetzt sitzen.)

Musterlösung zu Ü 5.4

Das Magnetfeld wird mit derselben Glockenform belegt, ausgewertet an den B-Positionen (halbe Zelle versetzt) und durch \(c\) geteilt:

# von oben: X, DX, DT; c (Import)
x_b = X[:-1] + DX/2                        # die Orte der B-Punkte
e = np.exp(-((X - 5.0)/0.5)**2)
b = np.exp(-((x_b + c*DT/2 - 5.0)/0.5)**2) / c   # B = E/c, rechtslaufend

fig, ax = plt.subplots(figsize=(7.0, 3.0))
ax.plot(X, e, "gray", label="t = 0")
for n in range(1, 301):
    b -= DT/DX * np.diff(e)
    e[1:-1] -= c**2 * DT/DX * np.diff(b)
ax.plot(X, e, "tab:red", label=f"t = {300*DT*1e9:.1f} ns")
ax.set_xlabel("x (m)"); ax.set_ylabel("E (V/m)")
ax.legend(); plt.tight_layout(); plt.show()
print(f"Amplitude: {np.max(e):.4f}   Ort des Maximums: {X[np.argmax(e)]:.2f} m")

Amplitude: 1.0000   Ort des Maximums: 8.00 m

Ein einziger Puls voller Höhe, der bei \(x = 8\) m angekommen ist — nichts läuft nach links. (Der kleine Zusatzterm \(c\Delta t/2\) im B-Argument berücksichtigt, dass das B-Feld im Leapfrog eine halbe Zeit-Stufe vor dem Start lebt; ohne ihn bliebe ein winziger Linksläufer übrig.)

Das Kleingedruckte

Die Ränder unseres Gitters halten \(E = 0\) fest — das ist physikalisch eine ideale Metallwand. Sobald ein Puls sie erreicht, wird er reflektiert; in diesem Kapitel haben wir vorher gestoppt. Ränder, die Wellen schlucken statt spiegeln, sind ein eigenes Thema (Kapitel 6 und 9).
„B = 0 am Start” heißt im Leapfrog genau genommen \(B(-\Delta t/2) = 0\) — eine halbe Zeitstufe vor dem ersten Bild. Dieser Versatz verschiebt die geteilten Pulse um eine halbe Zelle; sichtbar wird das erst bei Präzisionsmessungen (deshalb beginnt die Tempomessung im Programm erst nach der sauberen Trennung).
Der „magische” Zeitschritt \(\Delta t = \Delta x/c\) ist eine 1D-Spezialität — im Yee-Gitter in 2D und 3D gibt es ihn nicht, dort ist numerische Dispersion unvermeidlich und will gemessen werden (Kapitel 8/9).
Diskrete Energie: Die hier summierte Energie mischt \(E\) und \(B\) von leicht verschiedenen Zeitpunkten; deshalb schwankt sie im Promillebereich statt exakt konstant zu sein. Eine exakt erhaltene Energie kann man für das Schema definieren — für unsere Diagnose reicht die einfache Summe.