import numpy as np
import matplotlib.pyplot as plt
DX = 0.01 # Gitterabstand (m)
def leapfrog(s, nx, schritte, e0=None, u0=None, quelle=None):
"""Der Kern aus Kapitel 5-7 (u = c·B), Ränder PEC, weiche Quelle bei i=5."""
e = np.zeros(nx) if e0 is None else e0.copy()
u = np.zeros(nx - 1) if u0 is None else u0.copy()
for n in range(schritte):
u -= s * np.diff(e)
e[1:-1] -= s * np.diff(u)
if quelle is not None:
e[5] += quelle(n)
return e8 Stabilität und numerische Dispersion
Zwei Rechnungen aus den letzten Kapiteln sind offen geblieben. In Kapitel 5 lief der Puls nicht ungefähr mit Lichtgeschwindigkeit, sondern exakt — auf sechzehn Stellen; wir haben das dem „magischen Zeitschritt” \(\Delta t = \Delta x/c\) zugeschrieben und die Erklärung vertagt. Und in Übung 7.4 hat NumPy gemessen, dass der größte Eigenwert der Schritt-Matrix bis \(S = 1\) stur auf \(1{,}000000\) bleibt und bei \(S = 1{,}05\) auf \(1{,}87\) springt — warum kippt es ausgerechnet bei der Courant-Zahl \(S = c\,\Delta t/\Delta x = 1\), und warum so schlagartig?
Dieses Kapitel beantwortet beides mit einer einzigen Formel. Sie ist kurz genug für einen Kaffeebecher-Aufdruck, und sie hat zwei Gesichter: Eines verbietet zu große Zeitschritte (sonst Explosion), das andere bestraft zu grobe Gitter (die Wellen laufen dann zu langsam — und zwar gemessen zu langsam). Am Ende weißt du nicht nur, warum Simulationen explodieren, sondern auch, wie fein dein Gitter sein muss, damit du den Ergebnissen trauen darfst.
Lernziele
Nach diesem Kapitel kannst du …
- … anschaulich begründen, warum \(S \le 1\) sein muss (die Welle darf pro Zeitschritt höchstens eine Zelle weit kommen),
- … die gefährlichste Gitterwelle — den Zickzack mit \(\lambda = 2\Delta x\) — von Hand durch den Zeitschritt schicken und aus ihrer \(2\times 2\)-Matrix die CFL-Grenze exakt bei \(S = 1\) herleiten,
- … numerische Dispersion vorführen und deuten: Kurze Gitterwellen laufen zu langsam, ein scharfer Puls zieht eine Wellenschleppe,
- … die Dispersionsrelation \(\sin(\omega\Delta t/2) = S\,\sin(k\Delta x/2)\) herleiten und gegen Messungen bei 4, 8 und 16 Punkten pro Wellenlänge prüfen,
- … erklären, was den magischen Zeitschritt magisch macht — und mit Faustregeln ein Gitter wählen, dem man trauen kann.
8.1 Das Tempolimit, anschaulich
Ein Blick auf die zwei Update-Zeilen aus Kapitel 5 genügt für eine erste, sehr handfeste Beobachtung:
u -= S * np.diff(e) # u sieht nur seine zwei e-Nachbarn
e[1:-1] -= S * np.diff(u) # e sieht nur seine zwei u-NachbarnJeder Gitterpunkt spricht pro Zeitschritt nur mit seinen direkten Nachbarn. Eine Information — sagen wir: „hier ist gerade ein Puls losgelaufen” — kommt im Programm also pro Schritt höchstens eine Zelle weit, egal was wir rechnen. Die echte Welle dagegen legt in der Zeit \(\Delta t\) die Strecke \(c\,\Delta t\) zurück. Wenn \(c\,\Delta t > \Delta x\) ist, also \(S > 1\), dann ist die Wahrheit schneller als das Schema: Das Programm müsste Feldwerte an Orten ändern, von denen es noch gar nichts wissen kann. Das kann nicht gut gehen — die Bedingung
\[S = \frac{c\,\Delta t}{\Delta x} \le 1\]
heißt Courant-Bedingung (nach Richard Courant, auch CFL-Bedingung nach Courant, Friedrichs und Lewy, 1928 — älter als jeder Computer).
Das Argument ist ehrlich, aber unvollständig. Es sagt, dass es jenseits der Grenze nicht gehen kann — aber nicht, wie das Versagen aussieht, und auch nicht, ob nicht vielleicht schon bei \(S = 0{,}9\) Schluss ist. Dafür müssen wir hinsehen — mit diesem Versuchsaufbau:
Die Bühne ist die vertraute 1D-Strecke: 400 Zellen à 1 cm, also 4 m, mit festen Rändern (PEC, Kapitel 6) — die Funktion leapfrog unten ist wörtlich der Kern aus den Kapiteln 5–7 in der \(u = c\,B\)-Schreibweise, in der \(S\) der einzige Koeffizient ist. Anfangsbedingung: ein glatter Gauß-Puls im E-Feld in der Mitte bei \(x = 2\) m (Breite 12 cm), \(u = 0\) — er wird sich wie in Kapitel 5 in zwei Hälften teilen, aber darum geht es diesmal nicht. Der einzige Streichholz-Moment: Wir setzen \(S = 1{,}05\) — der Zeitschritt ist fünf Prozent größer, als die Courant-Bedingung erlaubt. Gemessen wird das Feldbild nach 60 und nach 120 Schritten, einmal im Überblick und einmal unter der Lupe (einzelne Gitterpunkte!), dazu zwei Kennzahlen: der Maximalbetrag des Felds und die Nachbar-Korrelation — das mittlere Produkt benachbarter Feldwerte, normiert auf den mittleren quadrierten Wert. Sie ist unser Form-Detektor: \(+1\) hieße „Nachbarn sind praktisch gleich” (glattes Feld), \(-1\) hieße „jeder Punkt ist das Negativ seines Nachbarn”.
# von oben: leapfrog(), DX
nx = 400
x = np.arange(nx) * DX
e_start = np.exp(-((x - 2.0) / 0.12)**2) # glatter Puls bei x = 2 m
e_60 = leapfrog(1.05, nx, 60, e0=e_start)
e_120 = leapfrog(1.05, nx, 120, e0=e_start)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7.2, 3.0))
ax1.plot(x, e_start, "gray", lw=1.2, label="Start")
ax1.plot(x, e_60, "tab:red", lw=0.8, label="nach 60 Schritten")
ax1.set_xlim(1.4, 2.6)
ax1.set_xlabel("x (m)"); ax1.set_ylabel("E (V/m)")
ax1.legend(fontsize=8)
ax2.plot(x, e_60, "o-", color="tab:red", ms=3.5, lw=0.7)
ax2.set_xlim(1.94, 2.06)
ax2.set_xlabel("x (m)")
ax2.set_title("Lupe: jeder Punkt ein Gitterwert", fontsize=9)
plt.tight_layout(); plt.show()
zacken = e_120[150:250] # Region um den alten Puls
korrelation = np.mean(zacken[1:] * zacken[:-1]) / np.mean(zacken**2)
print(f"nach 120 Schritten: max|E| = {np.max(np.abs(e_120)):.1e} V/m")
print(f"Nachbar-Korrelation: {korrelation:+.3f} (perfekter Zickzack: -1)")
nach 120 Schritten: max|E| = 3.5e+17 V/m
Nachbar-Korrelation: -1.009 (perfekter Zickzack: -1)Weder wachsender Puls noch Zufallsrauschen: Das Feld organisiert sich in die kürzeste Welle, die das Gitter darstellen kann — Plus, Minus, Plus, Minus, ein Vorzeichenwechsel pro Zelle, also \(\lambda = 2\Delta x\). Der Form-Detektor schlägt voll aus: Die Nachbar-Korrelation liegt bei \(-1\), jeder Punkt ist das Negativ seines Nachbarn. Nach 120 Schritten hat dieser Zickzack die Größenordnung \(10^{17}\) erreicht, während der eigentliche Puls die ganze Zeit brav geblieben ist.
Die Explosion ist also kein allgemeines Chaos, sondern das Werk einer einzigen Welle. Das ist eine gute Nachricht: Eine einzelne Welle können wir von Hand durch den Zeitschritt schicken.
Bevor wir den Verdächtigen verhören, der Tathergang als Film. Links das Feld, zur Lesbarkeit in jedem Bild auf sein Maximum normiert (sonst wäre nach der Explosion nichts mehr zu erkennen); rechts der Pegel \(\max|E|\) auf der Logarithmus-Skala, mit dem Fahrplan \(10^{-15}\cdot 1{,}877^n\) als gestrichelter Geraden — die \(1{,}877\) wird im nächsten Abschnitt hergeleitet, hier darf sie schon vorhersagen. Das Verstörende: Bis etwa Schritt 55 zeigt das linke Bild den völlig braven Puls, während rechts das Unheil längst hochgezählt wird — das Rundungsrauschen der Größe \(10^{-15}\) ist nur noch zu leise, um sichtbar zu sein. Dann überholt der Zickzack den Puls, und binnen weniger Bilder gehört ihm die ganze Szene. Exponentielles Wachstum kommt nicht schleichend — es kommt pünktlich.
Code der Animation (nur in der HTML-Fassung)
# von oben: np, plt, leapfrog(), e_start, x, nx
from matplotlib import animation
from IPython.display import HTML
S_F = 1.05
a_f = abs(1 - 2 * S_F**2)
EIGENWERT = a_f + np.sqrt(a_f**2 - 1) # 1,877 (nächster Abschnitt)
pegel_alle = np.array([np.max(np.abs(leapfrog(S_F, nx, n, e0=e_start)))
for n in range(121)])
TAKTE = range(0, 121, 4)
filme = [leapfrog(S_F, nx, n, e0=e_start) for n in TAKTE]
fig_a, (lnks, rchts) = plt.subplots(1, 2, figsize=(7.6, 3.2),
gridspec_kw={"width_ratios": [1.3, 1]})
form, = lnks.plot(x, filme[0] / np.max(np.abs(filme[0])), lw=0.8,
color="tab:red")
lnks.set_xlim(1.2, 2.8); lnks.set_ylim(-1.15, 1.15)
lnks.set_xlabel("x (m)"); lnks.set_ylabel("E (auf 1 normiert)")
n_alle = np.arange(121)
rchts.semilogy(n_alle, 1e-15 * EIGENWERT**n_alle, "k--", lw=1.0,
label="Fahrplan $10^{-15}\\cdot 1{,}877^n$")
spur, = rchts.semilogy([], [], color="tab:red", lw=1.6,
label="$\\max|E|$ gemessen")
rchts.set_xlim(0, 120); rchts.set_ylim(1e-16, 1e19)
rchts.set_xlabel("Schritt n")
rchts.legend(fontsize=7, loc="upper left")
def zeichne(j):
n = TAKTE[j]
form.set_ydata(filme[j] / np.max(np.abs(filme[j])))
spur.set_data(n_alle[:n + 1], pegel_alle[:n + 1])
lnks.set_title(f"Schritt {n} — max|E| = {pegel_alle[n]:.1e}",
fontsize=10)
return [form, spur]
anim = animation.FuncAnimation(fig_a, zeichne, frames=len(filme),
interval=110)
plt.close(fig_a)
HTML(anim.to_jshtml(default_mode="loop"))8.2 Der Hauptverdächtige: die Zickzackwelle
Die Zickzackwelle hat an den \(e\)-Punkten die Werte \(E, -E, E, -E,
\dots\) und an den dazwischen sitzenden \(u\)-Punkten die Werte \(U, -U,
U, -U, \dots\) — zwei Zahlen genügen, um sie komplett zu beschreiben. Was macht np.diff mit so einem Muster? Die Handrechnung mit \(E = 1\):
zick = np.array([1.0, -1.0, 1.0, -1.0, 1.0])
print("Zickzack: ", zick)
print("np.diff: ", np.diff(zick))Zickzack: [ 1. -1. 1. -1. 1.]
np.diff: [-2. 2. -2. 2.]Jede Differenz ist \(-1 - 1 = -2\) oder \(1 - (-1) = +2\), im Wechsel: Aus dem Zickzack wird wieder ein Zickzack, nur mit dem Faktor \(-2\). Die Form bleibt unter dem Update also erhalten; was sich ändert, sind nur die zwei Amplituden \(E\) und \(U\). Setzen wir das in die Update-Zeilen ein (u -= S*np.diff(e) macht aus dem \(e\)-Zickzack der Stärke \(E\) einen Beitrag \(+2SE\) zum \(u\)-Zickzack, und entsprechend für die zweite Zeile, in der bereits das neue \(u\) steht):
\[U' = U + 2S\,E, \qquad E' = E - 2S\,U' = (1 - 4S^2)\,E - 2S\,U.\]
Ein Zeitschritt ist für diese Welle also eine \(2\times 2\)-Matrix, die auf das Zahlenpaar \((E, U)\) wirkt:
\[\begin{pmatrix}E'\\U'\end{pmatrix} = \underbrace{\begin{pmatrix}1-4S^2 & -2S\\ 2S & 1\end{pmatrix}}_{M(S)} \begin{pmatrix}E\\U\end{pmatrix}.\]
Ob die Welle wächst, entscheiden die Eigenwerte von \(M\) — die Streckfaktoren, mit denen wir schon in Übung 7.4 gearbeitet haben: Gibt es ein \(\lambda\) mit \(|\lambda| > 1\), wird der zugehörige Anteil bei jedem Schritt um diesen Faktor verstärkt. Für eine \(2\times 2\)-Matrix \(\bigl(\begin{smallmatrix}a&b\\c&d \end{smallmatrix}\bigr)\) führt die Eigenwert-Bedingung auf eine gewöhnliche quadratische Gleichung:
\[\lambda^2 - T\,\lambda + D = 0, \qquad T = a + d \;\text{(die \textbf{Spur})}, \quad D = ad - bc \;\text{(die Determinante)}.\]
Beides können wir für \(M(S)\) direkt hinschreiben — die Handrechnung, hier für \(S = 1{,}05\), jeder Zwischenwert eine nachprüfbare Zahl:
S = 1.05
a, b_, c_, d = 1 - 4*S**2, -2*S, 2*S, 1.0
T = a + d # Spur: -3,41 + 1 = -2,41
D = a*d - b_*c_ # Det.: -3,41 + 4,41 = 1
diskriminante = (T/2)**2 - D # 1,4520 - 1 = 0,4520
lam1 = T/2 + np.sqrt(diskriminante) # -1,205 + 0,672
lam2 = T/2 - np.sqrt(diskriminante) # -1,205 - 0,672
print(f"Spur T = {T:.4f}, Determinante D = {D:.4f}")
print(f"Eigenwerte: {lam1:.4f} und {lam2:.4f}, Produkt = {lam1*lam2:.4f}")Spur T = -2.4100, Determinante D = 1.0000
Eigenwerte: -0.5327 und -1.8773, Produkt = 1.0000Zwei Befunde stecken in diesen Zahlen. Erstens ist die Determinante exakt 1 — und das gilt für jedes \(S\), denn \(D = (1-4S^2)\cdot 1 - (-2S)(2S) = 1 - 4S^2 + 4S^2 = 1\). Weil das Produkt der Eigenwerte gleich der Determinante ist, gilt also immer \(\lambda_1\lambda_2 = 1\): Was die Matrix in einer Richtung streckt, muss sie in der anderen stauchen. Es gibt nur zwei Möglichkeiten, diese Bilanz zu erfüllen — entweder ein echtes Streck-Stauch-Paar wie hier (\(-1{,}877\) und \(-0{,}533\)), oder beide Eigenwerte haben den Betrag exakt 1.
Zweitens: Welcher der beiden Fälle eintritt, entscheidet allein die Spur. Die Wurzel \(\sqrt{(T/2)^2 - 1}\) liefert nur dann eine reelle Zahl — also einen echten Streckfaktor —, wenn \(|T| \ge 2\) ist. Für \(|T| < 2\) gibt es keinen reellen Eigenwert: keinen Streckfaktor, keine Richtung, die einfach nur verstärkt würde. Das Zahlenpaar \((E, U)\) wird dann bei jedem Schritt gedreht statt gestreckt — \(E\) und \(U\) schaukeln ineinander wie Auslenkung und Schwung eines Pendels, der Betrag bleibt. (NumPy bestätigt das: np.linalg.eigvals liefert für \(S \le 1\) ein komplexes Zahlenpaar — Drehungen statt Streckungen —, dessen Betrag np.abs exakt 1 ist; so kam Übung 7.4 zu ihrer Messung \(1{,}000000\).)
Mit \(T = 2 - 4S^2\) wird aus \(|T| \le 2\) die Bedingung \(-2 \le 2 - 4S^2\), also \(S^2 \le 1\). Da ist sie — die CFL-Grenze, exakt bei \(S = 1\), nicht ungefähr:
\[\max_S |\lambda| = 1 \;\text{ für } S \le 1, \qquad |\lambda| = \frac{|T| + \sqrt{T^2-4}}{2} > 1 \;\text{ für } S > 1.\]
# von oben: leapfrog() ungenutzt; T-Formel aus dem Text
def verstaerkung(s):
spur = 2 - 4 * s**2
if abs(spur) <= 2:
return 1.0
return (abs(spur) + np.sqrt(spur**2 - 4)) / 2
s_achse = np.linspace(0, 1.2, 400)
fig, ax = plt.subplots(figsize=(6.4, 3.2))
ax.plot(s_achse, [verstaerkung(s) for s in s_achse], "tab:blue",
label="Zickzack-Welle (Spurformel)")
ax.plot([0.5, 1.0, 1.05], [1.0, 1.0, 1.874893], "ko", ms=6,
label="gemessen in Ü 7.4 (volle Matrix)")
ax.axvline(1.0, color="gray", ls=":", lw=0.8)
ax.set_xlabel("Courant-Zahl S = cΔt/Δx")
ax.set_ylabel("max|λ| pro Zeitschritt")
ax.legend()
plt.tight_layout(); plt.show()
print(f"Spurformel bei S = 1,05: {verstaerkung(1.05):.6f}"
f" (Ü 7.4 maß: 1.874893)")
Spurformel bei S = 1,05: 1.877328 (Ü 7.4 maß: 1.874893)Die Spurformel sagt \(1{,}8773\), Übung 7.4 hat an der vollen Matrix \(1{,}8749\) gemessen — eine Abweichung von gut einem Promille. Sie ist kein Fehler, sondern Fensterglas zum Verständnis: Auf dem endlichen Gitter mit seinen zwei Metallwänden passt der perfekte Zickzack nicht ganz hinein, die randnächsten Punkte tanzen leicht aus der Reihe. Je größer das Gitter, desto näher rückt der Messwert an die Formel (bei \(n = 320\) Zellen sind es schon \(1{,}8772\)).
Und warum explodiert dann jede Simulation, auch die mit einem glatten, harmlosen Puls? Das hat Übung 5.2 schon beantwortet: Das Rundungsrauschen der Gleitkommazahlen enthält jede Welle mikroskopisch — auch den Zickzack, mit einer Amplitude um \(10^{-15}\). Bei \(S = 1{,}05\) wird er pro Schritt ver-\(1{,}877\)-facht: \(10^{-15} \cdot 1{,}877^{120} \approx 7\cdot 10^{17}\). Gemessen haben wir oben \(3{,}5\cdot 10^{17}\) — die Explosion lief nach Fahrplan.
WarnungNaheliegende Vermutung
„Dann nehme ich zur Sicherheit \(S = 0{,}5\). Ein ordentlicher Sicherheitsabstand zur Grenze — und kleinere Zeitschritte sind ja sowieso genauer.”
Warum sie naheliegt: Überall sonst stimmt das. In Kapitel 2 wurde jede Ableitung genauer, je kleiner die Schrittweite; Sicherheitsabstand zu einer Abbruchkante ist generell eine vernünftige Idee; und „kleineres \(\Delta t\) = feinere Zeitauflösung = besser” klingt wie ein Naturgesetz.
Was stattdessen stimmt: Stabil ist \(S = 0{,}5\) — aber nicht genauer, sondern ungenauer. Unterhalb der Grenze laufen die Gitterwellen zu langsam, und zwar umso falscher, je kleiner \(S\) ist. Der Sicherheitsabstand wird mit verfälschten Geschwindigkeiten bezahlt; in 1D ist die Grenze \(S = 1\) selbst der genaueste Punkt. Der Rest des Kapitels führt das vor — und misst es aus.
8.3 Unter dem Limit: Wellen im falschen Tempo
# von oben: leapfrog(), DX
nx = 2000
x = np.arange(nx) * DX
sigma, x0 = 0.03, 3.0 # schmaler Puls (3 Zellen breit)
ergebnis = {}
for s in (1.0, 0.5):
e0 = np.exp(-((x - x0) / sigma)**2)
# Rechtsläufer wie in Ü 5.4: u = c·B = E an den u-Positionen,
# eine halbe Zeitstufe vor dem Start ausgewertet
x_u = x[:-1] + DX / 2
u0 = np.exp(-((x_u + s * DX / 2 - x0) / sigma)**2)
ergebnis[s] = leapfrog(s, nx, int(12.0 / DX / s), e0=e0, u0=u0)
fig, ax = plt.subplots(figsize=(7.0, 3.2))
ax.plot(x, np.exp(-((x - x0) / sigma)**2), "gray", lw=1.2, label="Start (3 m)")
ax.plot(x, ergebnis[1.0], "tab:blue", lw=1.2, label="S = 1")
ax.plot(x, ergebnis[0.5], "tab:orange", lw=1.2, label="S = 0,5")
ax.set_xlim(2.0, 16.5)
ax.set_xlabel("x (m)"); ax.set_ylabel("E (V/m)")
ax.legend()
plt.tight_layout(); plt.show()
print(f"Spitze nach 12 m: S = 1: {np.max(ergebnis[1.0]):.6f}")
print(f" S = 0,5: {np.max(ergebnis[0.5]):.4f}")
print(f"Schleppe (max |E| weit hinter dem Puls, S = 0,5): "
f"{np.max(np.abs(ergebnis[0.5][700:1400])):.4f}")
Spitze nach 12 m: S = 1: 1.000000
S = 0,5: 0.5343
Schleppe (max |E| weit hinter dem Puls, S = 0,5): 0.0380Mit \(S = 1\) kommt der Puls perfekt an: Spitze \(1{,}000000\), und hinter ihm ist das Feld auf Maschinengenauigkeit null — das ist der magische Zeitschritt aus Kapitel 5, und seine Erklärung sind wir gleich schuldig. Mit \(S = 0{,}5\) dagegen ist die Spitze auf \(0{,}53\) gesackt, und hinter dem Puls kräuselt eine Wellenschleppe her.
Die Schleppe verrät den Mechanismus. Ein schmaler Puls ist — das werden wir in Kapitel 10 mit der Fourier-Analyse präzise machen, die Idee genügt hier — ein Chor aus vielen Sinuswellen verschiedener Wellenlängen, die exakt im Takt starten. Im Vakuum bleibt der Chor für immer im Takt, weil alle Wellen mit demselben \(c\) laufen. Auf dem Gitter mit \(S < 1\) offenbar nicht: Die langen Wellen laufen fast richtig und bilden den (gerundeten) Hauptpuls, die kurzen hinken hinterher und tröpfeln als Schleppe nach. Wenn unterschiedliche Wellenlängen unterschiedlich schnell laufen, heißt das Dispersion — hier numerische Dispersion, denn sie ist kein Naturphänomen, sondern ein Artefakt des Gitters. (Echte, physikalische Dispersion gibt es auch — Glas macht genau das mit Licht, daher zerlegt ein Prisma Weiß in Farben; davon handelt Kapitel 16. Das Gitter äfft sie nur nach.)
In dieser HTML-Fassung kannst du der Schleppe beim Entstehen zusehen — der \(S = 0{,}5\)-Puls verliert von Schritt zu Schritt Höhe, und die kurzen Wellen tröpfeln sichtbar hinter ihm her:
Code der Animation (nur in der HTML-Fassung)
# von oben: np, plt, DX; FuncAnimation/HTML wie in Kapitel 4/5
from matplotlib import animation
from IPython.display import HTML
s_a, nx_a = 0.5, 2000
x_a = np.arange(nx_a) * DX
e_a = np.exp(-((x_a - 3.0) / 0.03) ** 2)
xu_a = x_a[:-1] + DX / 2
u_a = np.exp(-((xu_a + s_a * DX / 2 - 3.0) / 0.03) ** 2)
filmbilder = [e_a.copy()]
for n in range(2400):
u_a -= s_a * np.diff(e_a)
e_a[1:-1] -= s_a * np.diff(u_a)
if n % 40 == 39:
filmbilder.append(e_a.copy())
fig_a, ax_a = plt.subplots(figsize=(7.0, 2.8))
linie, = ax_a.plot(x_a, filmbilder[0], color="tab:orange", lw=1.2)
ax_a.set_xlim(2.0, 16.5)
ax_a.set_ylim(-0.35, 1.05)
ax_a.set_xlabel("x (m)")
ax_a.set_ylabel("E (V/m)")
def zeichne(i):
linie.set_ydata(filmbilder[i])
ax_a.set_title(f"S = 0,5 — Schritt {i * 40}")
return [linie]
anim = animation.FuncAnimation(fig_a, zeichne,
frames=len(filmbilder), interval=80)
plt.close(fig_a)
HTML(anim.to_jshtml(default_mode="loop"))Und das Tempolimit selbst steht hier als Live-Zelle bereit (editierbar, Run-Knopf oder Strg+Enter). Die eine Stellschraube \(S\) regiert beide Welten: S = 1.0 ist der magische Schritt, S = 0.5 die Schleppe — und S = 1.005 schon die verbotene Seite. Sage vorher an, wie schnell es dort wächst (die Zickzack-Theorie aus diesem Kapitel sagt es exakt voraus — der Vergleich steht im Titel):
„Kurze Wellen sind zu langsam” ist bisher eine Beobachtung. Jetzt machen wir sie zu einer Formel mit Zahlen, die man messen kann.
8.4 Die Dispersionsrelation des Gitters
Die Zickzack-Rechnung war so bequem, weil die Welle unter dem Update ihre Form behielt und nur Amplituden übrig blieben. Derselbe Trick funktioniert für jede Gitterwelle. Wir setzen als Lösungs-Ansatz eine laufende Welle an, wie wir sie seit Kapitel 4 schreiben — an den \(e\)-Punkten und den (räumlich wie zeitlich versetzten) \(u\)-Punkten jeweils mit eigener Amplitude:
\[e(x, t) = E\,\sin(kx - \omega t), \qquad u(x, t) = U\,\sin(kx - \omega t).\]
Gesucht ist, welche Kombinationen aus Wellenzahl \(k\) und Kreisfrequenz \(\omega\) das Update zulässt. Nehmen wir die erste Update-Zeile an einem \(u\)-Punkt am Ort \(x\): Sie verknüpft die \(u\)-Werte zu den Zeiten \(t \pm \Delta t/2\) mit den \(e\)-Werten an den Orten \(x \pm \Delta x/2\) zur Zeit \(t\):
\[U\sin\!\bigl(\Phi - \tfrac{\omega\Delta t}{2}\bigr) - U\sin\!\bigl(\Phi + \tfrac{\omega\Delta t}{2}\bigr) = -S\,\Bigl[E\sin\!\bigl(\Phi + \tfrac{k\Delta x}{2}\bigr) - E\sin\!\bigl(\Phi - \tfrac{k\Delta x}{2}\bigr)\Bigr],\]
wobei \(\Phi = kx - \omega t\) die Phase in der Mitte des Schritts ist. Auf beiden Seiten steht dieselbe Struktur: eine Differenz zweier Sinuswerte, deren Argumente symmetrisch um \(\Phi\) liegen. Dafür gibt es das passende Werkzeug aus der Schulmathematik, das Additionstheorem \(\sin(\Phi \pm \delta) = \sin\Phi\cos\delta \pm \cos\Phi\sin\delta\) — in der Differenz fällt der \(\sin\Phi\)-Anteil weg:
\[\sin(\Phi + \delta) - \sin(\Phi - \delta) = 2\cos\Phi\,\sin\delta.\]
Damit wird die linke Seite zu \(-2U\cos\Phi\,\sin(\omega\Delta t/2)\) und die rechte zu \(-2SE\cos\Phi\,\sin(k\Delta x/2)\) — und nun passiert das Schönste an der ganzen Rechnung: \(\cos\Phi\) kürzt sich weg. Ort und Zeit verschwinden aus der Gleichung; was an einem Punkt zu einem Zeitpunkt gilt, gilt automatisch überall und immer. Übrig bleibt
\[U\,\sin\frac{\omega\Delta t}{2} = S\,E\,\sin\frac{k\Delta x}{2}.\]
Die zweite Update-Zeile liefert mit derselben Rechnung das Spiegelbild \(E\,\sin(\omega\Delta t/2) = S\,U\,\sin(k\Delta x/2)\). Beide zusammen erzwingen \(U^2 = E^2\), also \(U = \pm E\) (Plus: Rechtsläufer, Minus: Linksläufer — der \(B = \pm E/c\)-Befund aus Kapitel 4 und 5, vom Gitter bestätigt), und vor allem die Dispersionsrelation des Gitters:
\[\boxed{\;\sin\frac{\omega\Delta t}{2} = S\,\sin\frac{k\Delta x}{2}\;}\]
Das ist die versprochene Formel mit den zwei Gesichtern. Sie verknüpft drei reine Zahlen (alles Skalare): den Zeit-Drehwinkel \(\omega\Delta t\), den Orts-Drehwinkel \(k\Delta x\) und die Courant-Zahl \(S\). Für jede Wellenzahl \(k\), die man hineinsteckt, legt sie fest, mit welcher Frequenz \(\omega\) — und damit welchem Tempo \(v = \omega/k\) — diese Welle auf dem Gitter läuft. Lesen wir sie dreimal:
Lesart 1 — die CFL-Grenze, noch einmal. Die linke Seite ist ein Sinus und kann höchstens 1 sein. Die rechte Seite wird am größten für die kürzeste Welle, den Zickzack (\(k\Delta x = \pi\), der Sinus wird 1): Dort verlangt die Gleichung \(\sin(\omega\Delta t/2) = S\). Für \(S > 1\) gibt es kein reelles \(\omega\), das das leistet — die Zickzackwelle kann nicht schwingen, ihr bleibt nur das, was wir oben in der \(2\times 2\)-Rechnung gefunden haben: wachsen. Spurkriterium und Dispersionsrelation sind dieselbe Physik in zwei Sprachen (die Eigenwert-Drehung pro Schritt ist der Winkel \(\omega\Delta t\)).
Lesart 2 — der magische Zeitschritt. Setze \(S = 1\): Dann steht da \(\sin(\omega\Delta t/2) = \sin(k\Delta x/2)\), also \(\omega\Delta t = k\Delta x\) — und nach Division durch \(\Delta t\): \(\omega = k\,\Delta x/\Delta t = ck\). Jede Welle, egal wie kurz, läuft exakt mit \(v = \omega/k = c\). Keine Dispersion, kein Formverlust, der Chor bleibt im Takt — darum kam der Puls oben unversehrt an, darum maß Kapitel 5 das Tempo auf Maschinengenauigkeit. Anschaulich: Bei \(S = 1\) rückt die Welle pro Schritt exakt eine Zelle weiter; das Schema muss nichts interpolieren, es darf durchreichen.
Lesart 3 — Dispersion unterhalb der Grenze. Für \(S < 1\) löst man nach \(\omega\) auf (\(\arcsin\) macht den Sinus rückgängig) und teilt durch \(k\):
\[\frac{v}{c} = \frac{\omega}{ck} = \frac{2\arcsin\!\bigl(S\,\sin\tfrac{k\Delta x}{2}\bigr)}{S\,k\Delta x}.\]
Weil der Sinus unterwegs „durchhängt” (er wächst langsamer als sein Argument), kommt für jedes \(k > 0\) etwas kleiner als 1 heraus — und zwar umso deutlicher, je größer \(k\Delta x\) ist, je weniger Gitterpunkte sich eine Wellenlänge also gönnt. Genau das ist die beobachtete Regel „kurze Wellen hinken”.
Wie groß ist der Effekt in Zahlen? Die Handrechnung für eine ordentlich aufgelöste Welle — 8 Punkte pro Wellenlänge — bei \(S = 0{,}5\). Die Quelle gibt die Frequenz vor, also rechnen wir von \(\omega\) nach \(k\) (gleiche Formel, rückwärts gelesen — jeder Zwischenwert zum Nachrechnen):
S = 0.5
w_dt = 2 * np.pi * S / 8 # Zeit-Drehwinkel der Quelle: 0,39270
schritt1 = np.sin(w_dt / 2) # sin(0,19635) = 0,19509
schritt2 = schritt1 / S # durch S = 0,39018
# np.arcsin: Umkehrung des Sinus (Bogenmass), Gegenstueck zu np.arccos
k_dx = 2 * np.arcsin(schritt2) # Orts-Drehwinkel = 0,80166
print(f"omega*dt = {w_dt:.5f} -> k*dx = {k_dx:.5f}")
print(f"Gitter-Wellenlänge: 2*pi/(k*dx) = {2*np.pi/k_dx:.4f} Zellen (statt 8)")
print(f"Tempo: v/c = (omega*dt)/(S*k*dx) = {w_dt/(S*k_dx):.5f}")omega*dt = 0.39270 -> k*dx = 0.80166
Gitter-Wellenlänge: 2*pi/(k*dx) = 7.8378 Zellen (statt 8)
Tempo: v/c = (omega*dt)/(S*k*dx) = 0.97972Die Welle, die physikalisch 8 Zellen lang sein müsste, wird auf dem Gitter auf \(7{,}84\) Zellen gestaucht und läuft mit \(0{,}980\,c\) — zwei Prozent zu langsam. Klingt harmlos; ob es das ist, klärt der Praxis-Abschnitt am Ende.
8.5 Messung gegen Theorie
Eine hergeleitete Formel bekommt bei uns keinen Vertrauensvorschuss (Kapitel 2 und 3 haben es vorgemacht): Wir messen nach. Der Versuchsaufbau kommt aus dem Wellenlabor von Kapitel 6 — eine weiche Quelle speist einen Sinus mit sanft anschwellender Amplitude ein, und wenn der Wellenzug eingeschwungen ist, vermessen wir ihn. Die Quelle diktiert die Frequenz; was sich das Gitter aussucht, ist die Wellenlänge. Die messen wir über die Abstände der Nulldurchgänge (linear interpoliert zwischen den zwei Gitterpunkten mit Vorzeichenwechsel, gemittelt über hunderte Durchgänge), denn aus \(v = f\lambda\) folgt direkt \(v/c = \lambda_\text{Gitter}/\lambda_0\).
# von oben: leapfrog(), DX
def tempo_theorie(s, n_lambda):
"""v/c laut Dispersionsrelation für eine Quelle mit λ₀ = n_lambda·Δx."""
w_dt = 2 * np.pi * s / n_lambda
k_dx = 2 * np.arcsin(np.sin(w_dt / 2) / s)
return w_dt / (s * k_dx)
def tempo_messung(s, n_lambda, nx=2600):
"""Sinusquelle einschwingen lassen, λ über Nulldurchgänge messen."""
w_dt = 2 * np.pi * s / n_lambda
anlauf = int(4 * 2 * np.pi / w_dt) # 4 Perioden sanft anfahren
quelle = lambda n: min(1.0, n / anlauf)**2 * np.sin(w_dt * n)
schritte = int((nx - 600) / s) # Front bleibt vor dem Rand
e = leapfrog(s, nx, schritte, quelle=quelle)
seg = e[400:int(s * schritte) - 600] # eingeschwungenes Fenster
# np.sign: Vorzeichen (-1/0/+1); np.flatnonzero: Indizes der Treffer
i = np.flatnonzero(np.sign(seg[:-1]) * np.sign(seg[1:]) < 0)
null = i - seg[i] / (seg[i + 1] - seg[i]) # linear interpoliert
lambda_gitter = 2 * np.mean(np.diff(null)) # in Zellen
return lambda_gitter / n_lambda # = v/c
print("Punkte/λ gemessen Theorie Abweichung")
for n_lambda in (4, 8, 16):
mess = tempo_messung(0.5, n_lambda)
theo = tempo_theorie(0.5, n_lambda)
print(f" {n_lambda:2} {mess:.4f} {theo:.4f} {abs(mess-theo):.1e}")Punkte/λ gemessen Theorie Abweichung
4 0.9013 0.9011 2.4e-04
8 0.9797 0.9797 1.3e-05
16 0.9951 0.9951 1.2e-06# von oben: tempo_theorie(), tempo_messung()
n_achse = np.linspace(3.0, 24, 300)
fig, ax = plt.subplots(figsize=(6.6, 3.4))
ax.plot(n_achse, [tempo_theorie(0.5, n) for n in n_achse], "tab:blue",
label="Theorie: Dispersionsrelation (S = 0,5)")
ax.plot([4, 8, 16], [tempo_messung(0.5, n) for n in (4, 8, 16)], "ko",
ms=7, label="gemessen (Nulldurchgänge)")
ax.axhline(1.0, color="gray", ls=":", lw=0.9, label="Lichtgeschwindigkeit c")
ax.set_xlabel("Punkte pro Wellenlänge")
ax.set_ylabel("v / c")
ax.legend(loc="lower right")
plt.tight_layout(); plt.show()
print(f"Gegenprobe magischer Zeitschritt — S = 1, 8 Punkte/λ: "
f"v/c = {tempo_messung(1.0, 8):.6f}")
Gegenprobe magischer Zeitschritt — S = 1, 8 Punkte/λ: v/c = 1.000000Messung und Theorie decken sich auf besser als ein Promille — die Herleitung hält. Und die Gegenprobe bestätigt Lesart 2: Bei \(S = 1\) liefert dieselbe Messmaschine \(v/c = 1{,}000000\).
Zwei Feinheiten am Rand der Abbildung verdienen einen Satz. Erstens endet die Theoriekurve links bei 3 Punkten pro Wellenlänge: Für noch kürzere Wellen verlangt die Relation \(\sin(k\Delta x/2) = \sin(\omega\Delta t/2)/S > 1\) — unmöglich, das Gitter verweigert den Transport (die Quelle erzeugt dann nur ein örtlich gefangenes Zappeln). Zweitens fällt der Fehler mit jeder Verdopplung der Auflösung auf etwa ein Viertel — \(9{,}9\,\%\), \(2{,}0\,\%\), \(0{,}49\,\%\). Das ist die Handschrift der zweiten Ordnung aus Kapitel 2: FDTD besteht aus zentralen Differenzen, und deren Fehler schrumpft mit dem Quadrat der Schrittweite.
8.6 Gitterwahl in der Praxis: Faustregeln
Wie fein muss das Gitter nun sein? Die ehrliche Antwort: Das hängt nicht nur von der Wellenlänge ab, sondern auch davon, wie weit die Welle laufen soll — der Tempofehler wirkt wie eine Uhr, die pro Tag drei Minuten nachgeht: Für eine Verabredung morgen ist das egal, nach einem Monat fehlt eine Stunde. Eine Welle, die mit \(v/c = 0{,}98\) über \(L\) Wellenlängen läuft, kommt um \(L\,(c/v - 1) \approx L \cdot 0{,}02\) Schwingungsperioden zu spät an; nach 50 Wellenlängen ist das eine ganze Periode — die Welle ist am Ziel komplett aus dem Takt, jede Interferenz-Vorhersage (Kapitel 14) wäre Unsinn. Daraus ergeben sich die Faustregeln:
| Regel | Begründung |
|---|---|
| Mindestens 10–20 Punkte pro kürzester Wellenlänge | 10 Punkte ≈ 1,3 % Tempofehler (bei \(S = 0{,}5\)) — gut für kurze Läufe |
| Lange Laufstrecken → feineres Gitter (Fehler wächst linear mit der Strecke) | \(L\) Wellenlängen × Tempofehler = Phasenfehler am Ziel (Ü 8.4) |
| \(S\) so nah an die Stabilitätsgrenze wie erlaubt | je größer \(S\), desto schwächer die Dispersion; in 1D ist \(S = 1\) exakt |
| Misstrauen bei scharfen Kanten | ein scharfer Puls enthält kurze Wellen — sie werden falsch transportiert (die Schleppe!) |
Ein Beispiel mit echten Einheiten: WLAN bei \(2{,}4\) GHz hat \(\lambda = 12{,}5\) cm. „16 Punkte pro Wellenlänge” heißt dann \(\Delta x \approx 8\) mm, und \(S = 1\) erzwingt \(\Delta t = \Delta x/c \approx 26\) ps — die Gitterwahl beginnt immer bei der kürzesten Wellenlänge im Problem, nie beim Speicherplatz.
Zum Schluss die Einschränkung, die schon mehrfach anklang: Das Privileg \(S = 1\) ist ein 1D-Sonderfall. In zwei und drei Dimensionen kann eine Welle auch diagonal übers Gitter laufen, und diagonal sind die Gitterpunkte weiter voneinander entfernt als achsenparallel — keine einzelne Zeitschrittwahl passt für beide Richtungen gleichzeitig. Die Stabilitätsgrenze rückt auf \(S \le 1/\sqrt{2}\) (2D) bzw. \(1/\sqrt{3}\) (3D), einen magischen Zeitschritt gibt es nicht mehr, und die Dispersion hängt zusätzlich von der Laufrichtung ab. Kapitel 9 führt das vor — die Werkzeuge von hier reichen dafür aus.
8.7 Das Kapitel-Programm
programme/kap08/kap08_dispersion.py bündelt alle vier Befunde eigenständig und prüft sie mit assert-Schranken: die Zickzack-Eigenwerte gegen die Spurformel (und gegen den Ü-7.4-Messwert), die Explosion als Zickzack (Nachbar-Korrelation \(< -0{,}9\), Amplitude \(> 10^{12}\)), den Pulsvergleich (bei \(S = 1\) Spitze auf \(10^{-9}\) exakt, bei \(S = 0{,}5\) Spitze \(< 0{,}6\) und Schleppe \(> 0{,}02\)) und die Tempomessung gegen die Dispersionsrelation (Abweichung \(< 5\cdot 10^{-4}\)).
Roter Faden
Drei alte Fäden sind hier zusammengelaufen: Die Fehler-schrumpft-quadratisch-Messung aus Kapitel 2 erklärt das \(1/N^2\) der Tempotabelle; der magische Zeitschritt aus Kapitel 5 ist eingelöst (\(S = 1 \Rightarrow \omega = ck\)); und die Eigenwert-Messung aus Übung 7.4 hat ihre Formel bekommen — der \(1{,}87\)-Sprung war die Zickzackwelle. Nach vorn: Kapitel 9 geht in zwei Dimensionen (\(S \le 1/\sqrt{2}\), richtungsabhängige Dispersion), und Kapitel 10 macht mit der Fourier-Analyse den „Chor aus Sinuswellen” quantitativ, der hier die Schleppe erklärt hat. Kapitel 14 (Interferenz) lebt davon, dass Phasen am Ziel stimmen — die Laufstrecken-Faustregel von hier ist dort Geschäftsgrundlage. Und in Kapitel 16 begegnet uns Dispersion als Physik: Glas tut mit Licht, was das Gitter mit unseren Wellen tut — nur dass es dort kein Artefakt ist, sondern Materialeigenschaft.
Übungen
Ü 8.1 (Verstehen). Wie schnell läuft die Zickzackwelle (\(\lambda = 2\Delta x\), also \(k\Delta x = \pi\)) bei \(S = 0{,}5\)? Rechne von Hand mit der Dispersionsrelation — als Kontrollwerte: \(\sin(\pi/2) = 1\) und \(\arcsin(1/2) = \pi/6\).
HinweisMusterlösung zu Ü 8.1
Die Relation verlangt \(\sin(\omega\Delta t/2) = S\cdot\sin(\pi/2) = 0{,}5\), also \(\omega\Delta t/2 = \pi/6\) und damit \(\omega\Delta t = \pi/3\). Das Tempo:
\[\frac{v}{c} = \frac{\omega\Delta t}{S\,k\Delta x} = \frac{\pi/3}{0{,}5\cdot\pi} = \frac{2}{3}.\]
Die kürzeste Gitterwelle läuft mit nur zwei Dritteln der Lichtgeschwindigkeit — der krasseste Fall der Tempolüge, noch unterhalb des linken Kurvenendes der Messabbildung (dort begann die Kurve erst bei 3 Punkten pro Wellenlänge; der Zickzack hat 2 und existiert nur als Anfangszustand, nicht als anregbare Welle). Die Kontrolle in einer Zeile:
print(f"v/c = {2*np.arcsin(0.5)/(0.5*np.pi):.6f} (2/3 = {2/3:.6f})")v/c = 0.666667 (2/3 = 0.666667)Ü 8.2 (Verstehen). Das Vakuum hat keine Dispersion — alle Wellenlängen laufen exakt mit \(c\). Das Gitter hat welche. Was besitzt das Gitter, das dem Vakuum fehlt und das diesen Unterschied erzwingt?
HinweisMusterlösung zu Ü 8.2
Einen eingebauten Maßstab. Das Gitter hat eine ausgezeichnete Länge (\(\Delta x\)) und eine ausgezeichnete Zeit (\(\Delta t\)); jede Welle wird daran gemessen — genau dafür steht die Größe „Punkte pro Wellenlänge”, und die Dispersionsrelation enthält \(k\) und \(\omega\) nur in den Kombinationen \(k\Delta x\) und \(\omega\Delta t\). Eine Welle mit 4 Punkten pro Wellenlänge ist auf dem Gitter etwas objektiv anderes als eine mit 16. Das Vakuum dagegen besitzt keinerlei eigene Länge oder Zeit: Es gibt nichts, womit es eine Wellenlänge vergleichen könnte, also kann es Wellen verschiedener Länge nicht verschieden behandeln — gleiches Tempo für alle ist die einzige maßstabsfreie Möglichkeit. (Glas hat wieder einen Maßstab — seine Atome und deren Resonanzen. Darum hat Glas Dispersion; Kapitel 16.)
Ü 8.3 (Verändern). Erweitere die Tempomessung um den Punkt 32 Punkte pro Wellenlänge und vergleiche den Fehler \(1 - v/c\) mit dem bei 16 Punkten. Welchen Quotienten erwartest du nach der „zweiten Ordnung” aus Kapitel 2 — und welchen misst du?
HinweisMusterlösung zu Ü 8.3
Verdopplung der Auflösung sollte den Fehler vierteln (Fehler \(\propto 1/N^2\)):
# von oben: tempo_messung() (Abschnitt „Messung gegen Theorie")
fehler_16 = 1 - tempo_messung(0.5, 16)
fehler_32 = 1 - tempo_messung(0.5, 32)
print(f"Fehler bei 16 Punkten/λ: {fehler_16:.5f}")
print(f"Fehler bei 32 Punkten/λ: {fehler_32:.5f}")
print(f"Quotient: {fehler_16/fehler_32:.2f} (zweite Ordnung: 4)")Fehler bei 16 Punkten/λ: 0.00488
Fehler bei 32 Punkten/λ: 0.00121
Quotient: 4.04 (zweite Ordnung: 4)Gemessen wird ein Quotient von rund 4 — die zweite Ordnung der zentralen Differenzen, jetzt zum dritten Mal im Buch bestätigt (Kapitel 2: einzelne Ableitung; hier: das Wellentempo). Wer den Fehler einer FDTD-Rechnung abschätzen will, kann also schlicht zweimal rechnen: Halbiert sich der Unterschied zur feineren Rechnung auf ein Viertel, ist man im sauberen Konvergenz-Regime.
Ü 8.4 (Übertragen). Du planst eine Simulation, in der ein Funksignal 50 Wellenlängen weit läuft, und erlaubst am Ziel höchstens ein Zehntel Periode Phasenfehler (\(S = 0{,}5\)). Schreibe einen kleinen Gitterplaner: Wie viele Punkte pro Wellenlänge brauchst du? Reicht die 10-Punkte-Faustregel?
HinweisMusterlösung zu Ü 8.4
Eine Welle mit Tempo \(v\) braucht für die Strecke \(L = 50\lambda\) die Zeit \(L/v\) statt \(L/c\) — die Verspätung, gemessen in Schwingungsperioden \(T = \lambda/c\), ist \(50\,(c/v - 1)\). Der Planer probiert Auflösungen durch, bis sie unter \(0{,}1\) fällt:
# von oben: tempo_theorie() (Abschnitt „Messung gegen Theorie")
for n in range(4, 60):
v = tempo_theorie(0.5, n)
verspaetung = 50 * (1 / v - 1) # in Perioden
if verspaetung < 0.1:
print(f"{n} Punkte pro Wellenlänge: Verspätung "
f"{verspaetung:.3f} Perioden (v/c = {v:.5f})")
break25 Punkte pro Wellenlänge: Verspätung 0.099 Perioden (v/c = 0.99802)Es braucht 25 Punkte pro Wellenlänge — die 10-Punkte-Faustregel reicht hier nicht (sie ergäbe \(50 \cdot 0{,}0129 \approx 0{,}64\) Perioden: das Signal käme mehr als eine halbe Schwingung verschoben an, bei Interferenz-Fragen der Unterschied zwischen Verstärkung und Auslöschung). Die Faustregel taugt für kurze Wege; bei langen zählt das Produkt aus Strecke und Tempofehler. Zwei Auswege, beide aus diesem Kapitel: feiner auflösen — oder \(S\) an die Grenze schieben (in 1D löst \(S = 1\) das Problem exakt; in 2D/3D, wo es keinen magischen Zeitschritt gibt, bleibt nur das feinere Gitter).