import numpy as np
import matplotlib.pyplot as plt
DX = 0.01 # Gitterabstand (m), Δy = Δx
S_2D = 1 / np.sqrt(2) # warum genau dieser Wert: gleich!
def felder(nx, ny):
"""Leere TMz-Felder: e auf den Kreuzungen, ux/uy auf den Kanten."""
e = np.zeros((nx, ny)) # E_z bei (i, j)
ux = np.zeros((nx, ny - 1)) # c·B_x bei (i, j+1/2) — eine Spalte weniger
uy = np.zeros((nx - 1, ny)) # c·B_y bei (i+1/2, j) — eine Zeile weniger
return e, ux, uy
def schritt(e, ux, uy, s):
"""Ein Yee-Zeitschritt in 2D; Ränder bleiben fest (PEC)."""
# axis-Argument: np.diff bildet Nachbar-Differenzen ENTLANG einer
# Richtung — axis=0 läuft über i (x), axis=1 über j (y)
ux -= s * np.diff(e, axis=1) # Faraday, x-Komponente
uy += s * np.diff(e, axis=0) # Faraday, y-Komponente
e[1:-1, 1:-1] += s * (np.diff(uy, axis=0)[:, 1:-1]
- np.diff(ux, axis=1)[1:-1, :]) # Ampère
return e, ux, uy9 FDTD in der Fläche: das volle Yee-Gitter
Vier Kapitel lang lebte unsere Welt auf einer Schnur. Das war kein Geiz, sondern Methode — jede Idee (Leapfrog, Quellen, Ränder, Stabilität, Dispersion) ließ sich in 1D sauber sezieren. Aber eine Schnur kennt nur links und rechts. Alles, was Funk und Optik spannend macht, braucht Richtung und Form: der Schatten hinter einer Wand, die Keule einer Antenne, die Beugung am Spalt, die Linse. Nichts davon existiert in 1D.
Dieses Kapitel macht den Sprung in die Fläche — und er ist kleiner, als er klingt: Der komplette 2D-Kern sind drei np.diff-Zeilen statt zwei. Dafür bekommen wir zum ersten Mal Feldbilder im Wortsinn (die Farbdarstellung, die Kapitel 4 versprochen hat), eine neue Stabilitätsgrenze an einer überraschenden Stelle, eine Dispersion, die von der Laufrichtung abhängt — und endlich die in Kapitel 6 zugesagte „volle Konstruktion” des reflexionsfreien Randes: die PML.
Lernziele
Nach diesem Kapitel kannst du …
- … das 2D-Yee-Gitter aufbauen (\(E_z\) auf den Kreuzungen, die beiden \(B\)-Komponenten auf den Kanten) und den Zeitschritt in drei
np.diff-Zeilen schreiben, - … eine Ringwelle simulieren, als Farbbild lesen und ihr Tempo gegen \(c\) messen,
- … erklären und vorführen, warum die Stabilitätsgrenze in 2D bei \(S \le 1/\sqrt{2}\) liegt — und welche Welle dort kippt (das Schachbrett),
- … die richtungsabhängige Dispersion messen: Bei \(S = 1/\sqrt{2}\) laufen Wellen diagonal exakt, entlang der Achsen zu langsam,
- … offene Ränder bauen: von der \(\sigma\)-Rampe zur Split-Feld-PML, mit gemessener Restreflexion.
9.1 Die halbe Welt genügt: das TMz-Gitter
In 3D haben die Felder sechs Komponenten (\(E_x, E_y, E_z\) und \(B_x, B_y, B_z\)). Wir betrachten eine Welt, die sich in \(z\)-Richtung nicht ändert — alles hängt nur von \(x\) und \(y\) ab, wie bei einem sehr langen geraden Draht oder einem Wellenbad zwischen zwei Glasplatten. Schreibt man die beiden Wirbelgleichungen aus Kapitel 3 komponentenweise auf und streicht alle \(z\)-Ableitungen, zerfallen die sechs Komponenten in zwei Dreiergruppen, die nichts voneinander wissen: \((E_z, B_x, B_y)\) und \((B_z, E_x, E_y)\). Es genügt, eine davon zu rechnen — wir nehmen die erste (sie heißt traditionell TM, „transversal magnetisch”: das E-Feld steht senkrecht auf der Rechenebene, das B-Feld liegt in ihr). Von den Wirbelgleichungen bleibt dann:
\[\frac{\partial E_z}{\partial t} = c^2\left(\frac{\partial B_y}{\partial x} - \frac{\partial B_x}{\partial y}\right), \qquad \frac{\partial B_x}{\partial t} = -\frac{\partial E_z}{\partial y}, \qquad \frac{\partial B_y}{\partial t} = +\frac{\partial E_z}{\partial x}.\]
Die erste Zeile ist die \(z\)-Komponente von „Ampère mit Verschiebungsstrom” — der vertraute Wirbelausdruck \(\partial_x B_y - \partial_y B_x\) aus Kapitel 2. Die anderen beiden sind Faraday, Komponente für Komponente.
Wo wohnen die drei Größen auf dem Gitter? In 1D saß \(B\) zwischen den \(E\)-Punkten, damit jede Ableitung eine zentrale Differenz mit Spannweite \(\Delta x\) wird (das Yee-Argument aus Kapitel 5). Dasselbe Prinzip, jetzt auf einem Stadtplan: \(E_z\) wohnt auf den Kreuzungen \((i, j)\). \(B_x\) braucht die \(y\)-Ableitung von \(E_z\), also wohnt es auf halbem Weg zwischen zwei Kreuzungen in \(y\)-Richtung, bei \((i,\, j{+}\tfrac12)\). \(B_y\) braucht die \(x\)-Ableitung, also sitzt es bei \((i{+}\tfrac12,\, j)\). Jede Differenz im Update verbindet damit exakt die zwei nächsten Nachbarn der Zielgröße — alles bleibt zentral, nichts wird übersprungen.
Mit der Abkürzung aus Kapitel 7 — wir rechnen mit \(u_x = c\,B_x\) und \(u_y = c\,B_y\), die dieselbe Einheit wie \(E_z\) haben — bleibt wieder die Courant-Zahl \(S = c\,\Delta t/\Delta x\) als einziger Koeffizient (wir setzen \(\Delta y = \Delta x\)):
Das ist der ganze Sprung von 1D nach 2D: Die Update-Zeile für \(e\) zieht jetzt zwei Differenzen heran (den Wirbel \(\partial_x u_y - \partial_y u_x\) statt einer einzelnen Ableitung), und das Magnetfeld hat zwei Komponenten statt einer. Drei Zeilen, keine Matrix, kein Löser — Kapitel 7 lässt grüßen: explizit.
Das erste Experiment soll zeigen, was dieses Gitter mit einer punktförmigen Störung macht. Der Versuchsaufbau: Die Bühne ist ein Quadrat von 8 m × 8 m (\(801 \times 801\) Zellen à 1 cm), alle Felder starten auf null; die Ränder sind PEC-Wände, aber weit genug weg, dass im Beobachtungszeitraum nichts dort ankommt. Als Courant-Zahl nehmen wir \(S = 1/\sqrt{2}\) — warum das die 2D-Grenze ist, klärt der nächste Abschnitt; hier genügt: Es ist stabil. Die Anregung: Statt einer Anfangsbedingung benutzen wir eine weiche Quelle im Zentrum (Kapitel 6: aufaddieren statt überschreiben), die genau ein zeitliches Gauß-Päckchen abgibt — ein kurzes „Pling” an einem einzigen Gitterpunkt, danach herrscht dort wieder Ruhe. Gemessen wird dreierlei: Schnappschüsse des ganzen Felds bei Schritt 180, 340 und 500; als eingezeichnete Erwartung in jedem Bild ein Kreis um die Quelle mit Radius \(c \cdot t\) (mit \(t\) ab dem Quellen-Maximum); und das Tempo der Störung mit der Maximum-Verfolgung aus Kapitel 5 — wir notieren laufend den Ort des Wellenkamms auf einem Strahl nach rechts und legen eine Gerade durch die Ort-Zeit-Punkte.
# von oben: felder(), schritt(), DX, S_2D
nx = 801
mitte = nx // 2
def gauss(n, t0=30.0, tau=9.0):
return np.exp(-((n - t0) / tau)**2)
e, ux, uy = felder(nx, nx)
schnappschuesse, zeiten, kamm_orte = {}, [], []
for n in range(521):
e, ux, uy = schritt(e, ux, uy, S_2D)
e[mitte, mitte] += gauss(n)
if n in (180, 340, 500):
schnappschuesse[n] = e.copy()
if n >= 150 and n % 30 == 0:
zeiten.append(n)
kamm_orte.append(np.argmax(e[mitte, mitte:])) # Kamm auf x-Strahl
fig, achsen = plt.subplots(1, 3, figsize=(7.4, 2.7))
ausdehnung = (0, nx * DX, 0, nx * DX)
winkel = np.linspace(0, 2 * np.pi, 200)
for ax, n in zip(achsen, schnappschuesse):
bild = schnappschuesse[n]
w = np.abs(bild).max()
ax.imshow(bild.T, origin="lower", extent=ausdehnung, cmap="RdBu_r",
vmin=-w / 3, vmax=w / 3)
radius = S_2D * (n - 30) * DX # c·t seit dem Quellen-Pling
ax.plot(mitte * DX + radius * np.cos(winkel),
mitte * DX + radius * np.sin(winkel), "k:", lw=0.9)
ax.set_title(f"Schritt {n}", fontsize=9)
ax.set_xlabel("x (m)")
achsen[0].set_ylabel("y (m)")
plt.tight_layout(); plt.show()
tempo = np.polyfit(zeiten, kamm_orte, 1)[0] / S_2D # Zellen/Schritt → v/c
print(f"Kamm-Tempo (Maximum-Verfolgung wie in Kapitel 5): {tempo:.4f} c")
Kamm-Tempo (Maximum-Verfolgung wie in Kapitel 5): 0.9967 cEin Ring, der mit Lichtgeschwindigkeit wächst — die Maximum-Verfolgung liefert \(0{,}997\,c\), und der Kamm sitzt in jedem Schnappschuss auf dem \(c\cdot t\)-Kreis. Niemand hat dem Programm „Ring” beigebracht; es kennt nur Nachbar-Differenzen. Dass aus einem Punkt-Pling ein perfekter Kreis wird, ist die 2D-Ausgabe derselben Physik, die in 1D den Puls in zwei Hälften teilte: Die Störung läuft in jede Richtung gleich schnell davon.
Und die Farbe im Innern? Weiß — fast. Hinter dem Ring kehrt das Feld auf nahezu null zurück (hier bleibt ein Nachzügler von etwa einem halben Prozent des Kamms). Das winzige Nachleuchten ist kein Programmierfehler, sondern ehrliche 2D-Physik: Unsere Ebene ist ja ein Schnitt durch eine 3D-Welt, in der die „Punktquelle” eine unendlich lange Linie ist — und deren entfernte Abschnitte liefern Beiträge, die später eintreffen. In echtem 3D (Kugelwellen) gäbe es das nicht. Fürs Auge ist es unsichtbar; wir erwähnen es, damit du ihm nicht irgendwann misstraust.
In dieser HTML-Fassung wächst der Ring als Film — mit dem eingezeichneten \(c\cdot t\)-Kreis als ständigem Begleiter (die Einzelschritt-Knöpfe zeigen schön, wie der Kamm nie vom Kreis abweicht):
Code der Animation (nur in der HTML-Fassung)
# von oben: np, plt, felder(), schritt(), DX, S_2D
from matplotlib import animation
from IPython.display import HTML
nx_a = 401
mitte_a = nx_a // 2
e_a, ux_a, uy_a = felder(nx_a, nx_a)
filmbilder = []
for n in range(261):
e_a, ux_a, uy_a = schritt(e_a, ux_a, uy_a, S_2D)
e_a[mitte_a, mitte_a] += np.exp(-((n - 30.0) / 9.0) ** 2)
if n % 10 == 0 and n >= 40:
filmbilder.append((n, e_a.copy()))
fig_a, ax_a = plt.subplots(figsize=(5.2, 5.2))
w = np.abs(filmbilder[len(filmbilder) // 2][1]).max() / 3
ausdehnung_a = (0, nx_a * DX, 0, nx_a * DX)
bild = ax_a.imshow(filmbilder[0][1].T, origin="lower",
extent=ausdehnung_a, cmap="RdBu_r",
vmin=-w, vmax=w)
winkel_a = np.linspace(0, 2 * np.pi, 200)
kreis, = ax_a.plot([], [], "k:", lw=0.9)
ax_a.set_xlabel("x (m)"); ax_a.set_ylabel("y (m)")
def zeichne(i):
n, feld = filmbilder[i]
bild.set_data(feld.T)
radius = S_2D * (n - 30) * DX
kreis.set_data(mitte_a * DX + radius * np.cos(winkel_a),
mitte_a * DX + radius * np.sin(winkel_a))
ax_a.set_title(f"Schritt {n} — Kreis: c·t")
return [bild, kreis]
anim = animation.FuncAnimation(fig_a, zeichne,
frames=len(filmbilder), interval=90)
plt.close(fig_a)
HTML(anim.to_jshtml(default_mode="loop"))9.2 Das neue Tempolimit
In 1D hieß die Regel \(S \le 1\), und der „magische” Wert \(S = 1\) war sogar exakt. Gilt das hier weiter?
WarnungNaheliegende Vermutung
„\(S \le 1\) gilt auch in 2D. Das Informationsargument aus Kapitel 8 stimmt doch weiterhin: Pro Zeitschritt rückt die Information eine Zelle — in \(x\) wie in \(y\).”
Warum sie naheliegt: Pro Achse ist das Argument korrekt — die Update-Zeilen verbinden nach wie vor nur direkte Nachbarn, und entlang einer Gitterachse kommt eine Information höchstens eine Zelle pro Schritt voran. Nichts daran wird in 2D falsch.
Was stattdessen stimmt: Der Engpass liegt diagonal. Ein Zeitschritt erreicht von einem Punkt aus nur die vier direkten Nachbarn — zusammen eine Raute \(|\Delta i| + |\Delta j| \le 1\). Die echte Welle breitet sich aber als Kreis aus, und der größte Kreis, der in diese Raute passt, hat nur den Radius \(\Delta x/\sqrt{2}\) (von der Mitte senkrecht zur schrägen Rautenkante). Diagonal schafft das Schema pro Schritt also nur \(\Delta x/\sqrt{2}\) — und damit muss \(c\,\Delta t \le \Delta x/\sqrt{2}\) sein:
\[S \le \frac{1}{\sqrt{2}} \approx 0{,}707.\]
Wie in Kapitel 8 ist das erst ein Plausibilitätsargument. Den scharfen Beweis — und die Form des Versagens — liefert wieder die gefährlichste Gitterwelle.
In Kapitel 8 haben wir die Dispersionsrelation des 1D-Gitters hergeleitet; in der quadrierten Form lautete sie \(\sin^2(\omega\Delta t/2) = S^2\sin^2(k\Delta x/2)\). Führt man dieselbe Rechnung mit dem 2D-Ansatz \(\sin(k_x x + k_y y - \omega t)\) durch, passiert genau das, was die dritte Update-Zeile vermuten lässt: Sie enthält zwei Differenzen, und jede steuert ihren eigenen Sinus-Term bei — die Beiträge der beiden Richtungen addieren sich:
\[\sin^2\frac{\omega\Delta t}{2} = S^2\left(\sin^2\frac{k_x\Delta x}{2} + \sin^2\frac{k_y\Delta x}{2}\right).\]
(Probe: Für \(k_y = 0\) — eine Welle, die exakt in \(x\)-Richtung läuft — steht da wieder die 1D-Formel.) Die linke Seite kann höchstens \(1\) werden. Die rechte wird am größten, wenn beide Sinus gleich \(1\) sind: \(k_x\Delta x = k_y\Delta x = \pi\), also Zickzack in \(x\)- und \(y\)-Richtung zugleich — ein Schachbrett, bei dem jeder Gitterwert das Negativ aller vier Nachbarn ist. Für diese Welle verlangt die Relation \(\sin(\omega\Delta t/2) = S\sqrt{2}\), und das geht nur reell, solange \(S\sqrt{2} \le 1\). Die Raute hatte recht, und zwar exakt.
Auch die Spur-Rechnung aus Kapitel 8 lässt sich wörtlich übertragen: Für das Schachbrett liefern die beiden Richtungen je einen \(-4S^2\)-Beitrag zur Spur der Schrittmatrix, \(T = 2 - 8S^2\), und die Handrechnung für \(S = 0{,}72\) (knapp über der Grenze, denn \(8\cdot 0{,}72^2 = 4{,}1472\)):
\[T = -2{,}1472, \qquad |\lambda|_{\max} = \frac{2{,}1472 + \sqrt{2{,}1472^2 - 4}}{2} = \frac{2{,}1472 + 0{,}7813}{2} = 1{,}464.\]
Der Praxistest dazu, wieder mit vollständigem Aufbau: ein kleines Gitter von \(101 \times 101\) Zellen (mehr braucht es nicht — die Explosion ist lokal), als Anfangsbedingung ein glatter 2D-Gauß-Hügel in der Mitte, PEC-Ränder, \(S = 0{,}72\) — nur zwei Prozent über der Grenze \(1/\sqrt{2} \approx 0{,}707\). Über 120 Schritte protokollieren wir den Maximalbetrag des Felds (aus seinen letzten 20 Schritten lesen wir den mittleren Wachstumsfaktor pro Schritt ab) und am Ende den Form-Detektor aus Kapitel 8, die Nachbar-Korrelation — jetzt zweimal, einmal für \(x\)- und einmal für \(y\)-Nachbarn. Danach läuft zur Gegenprobe derselbe Aufbau mit \(S = 0{,}70\), knapp unter der Grenze.
# von oben: felder(), schritt()
n_klein = 101
xx, yy = np.meshgrid(np.arange(n_klein), np.arange(n_klein), indexing="ij")
huegel = np.exp(-(((xx - 50)**2 + (yy - 50)**2) / 12.0**2))
e, ux, uy = felder(n_klein, n_klein)
e += huegel
pegel = []
for n in range(120):
e, ux, uy = schritt(e, ux, uy, 0.72)
pegel.append(np.abs(e).max())
wachstum = (pegel[-1] / pegel[-21])**(1 / 20) # mittlerer Faktor, 20 Schritte
fig, ax = plt.subplots(figsize=(4.6, 4.0))
ax.imshow(e[38:63, 38:63].T, origin="lower", cmap="RdBu_r")
ax.set_xlabel("Zelle i"); ax.set_ylabel("Zelle j")
plt.tight_layout(); plt.show()
brett = e[40:60, 40:60]
print(f"max|E| nach 120 Schritten: {pegel[-1]:.1e}")
print(f"Wachstum pro Schritt: {wachstum:.3f} (Spurformel: 1.464)")
print(f"Nachbar-Korrelation: x-Richtung "
f"{np.mean(brett[1:, :] * brett[:-1, :]) / np.mean(brett**2):+.3f}, "
f"y-Richtung "
f"{np.mean(brett[:, 1:] * brett[:, :-1]) / np.mean(brett**2):+.3f}")
e2, ux2, uy2 = felder(n_klein, n_klein)
e2 += huegel
for n in range(120):
e2, ux2, uy2 = schritt(e2, ux2, uy2, 0.70)
print(f"Gegenprobe S = 0,70 (unter der Grenze): "
f"max|E| = {np.abs(e2).max():.3f} — stabil")
max|E| nach 120 Schritten: 5.4e+02
Wachstum pro Schritt: 1.453 (Spurformel: 1.464)
Nachbar-Korrelation: x-Richtung -1.010, y-Richtung -1.010
Gegenprobe S = 0,70 (unter der Grenze): max|E| = 0.243 — stabilMuster und Faktor, beides nach Fahrplan: Schachbrett (Nachbar-Korrelation \(\approx -1\) in beiden Richtungen) mit Wachstum \(1{,}45\) pro Schritt, und zwei Hundertstel darunter, bei \(S = 0{,}70\), passiert — nichts. Die 1D-Werkzeuge aus Kapitel 8 tragen unverändert; nur die Zahl an der Grenze ist neu.
Das schmale Fenster zwischen „nichts” und „Schachbrett” kannst du in dieser Live-Zelle selbst abtasten (editierbar, Run-Knopf oder Strg+Enter — und ja, hier rechnet dein Browser ein komplettes 2D-Yee-Gitter). Taste dich heran: S = 0.70, 0.707, 0.71, 0.72 — wo genau kippt es, und wie sieht der Täter aus?
9.3 Die Magie wandert in die Diagonale
In 1D gab es den magischen Zeitschritt: Bei \(S = 1\) liefen alle Wellen exakt mit \(c\). Diese Tür ist in 2D zu — \(S = 1\) ist jenseits der Stabilitätsgrenze. Was passiert am neuen Maximum \(S = 1/\sqrt{2}\)? Die Dispersionsrelation kennt jetzt nicht nur eine Wellenlänge, sondern auch eine Laufrichtung, und die beiden Extremfälle lohnen die Handrechnung.
Entlang einer Achse (\(k_y = 0\)) gilt die 1D-Formel mit \(S = 1/\sqrt{2} < 1\): Dispersion wie in Kapitel 8, kurze Wellen hinken. Entlang der Diagonalen ist \(k_x = k_y = k/\sqrt{2}\) (der Wellenvektor der Länge \(k\) verteilt sich gleichmäßig auf beide Achsen), die beiden gleichen Sinus-Terme addieren sich, und die Relation wird zu
\[\sin\frac{\omega\Delta t}{2} = S\sqrt{2}\,\sin\!\left(\frac{k\Delta x}{2\sqrt{2}}\right).\]
Jetzt setze \(S = 1/\sqrt{2}\) ein: Der Vorfaktor \(S\sqrt{2}\) wird exakt 1, also \(\omega\Delta t = k\Delta x/\sqrt{2}\) — und mit \(\Delta t = \Delta x/(c\sqrt{2})\) folgt \(\omega = c\,k\) für jede Wellenlänge. Das ist Wort für Wort die Rechnung des magischen Zeitschritts aus Kapitel 8, nur dass sie jetzt der Diagonale gehört: Bei \(S = 1/\sqrt{2}\) laufen diagonale Wellen exakt, achsenparallele zu langsam. Anschaulich passt das zur Raute: Diagonal ist der Schritt des Schemas gerade so groß wie der Weg der Welle — dort wird durchgereicht statt interpoliert.
# von oben: felder(), schritt(), DX, S_2D
def lambda_aus_strahl(strahl, abstand):
"""Wellenlänge über Nulldurchgänge (Methode aus Kapitel 8)."""
i = np.flatnonzero(np.sign(strahl[:-1]) * np.sign(strahl[1:]) < 0)
null = i - strahl[i] / (strahl[i + 1] - strahl[i])
return 2 * np.mean(np.diff(null)) * abstand
nx = 701
m = nx // 2
n_lambda = 5
w_dt = 2 * np.pi * S_2D / n_lambda
anlauf = int(4 * 2 * np.pi / w_dt)
e, ux, uy = felder(nx, nx)
for n in range(460):
e, ux, uy = schritt(e, ux, uy, S_2D)
e[m, m] += min(1.0, n / anlauf)**2 * np.sin(w_dt * n)
a, b = 60, 260 # Messfenster (Zellen)
strahl_achse = e[m, m + a:m + b]
strahl_diag = np.array([e[m + k, m + k] # Schritt entlang i = j:
for k in range(int(a / np.sqrt(2)),
int(b / np.sqrt(2)))])
v_achse = lambda_aus_strahl(strahl_achse, 1.0) / n_lambda
v_diag = lambda_aus_strahl(strahl_diag, np.sqrt(2)) / n_lambda
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7.4, 3.4),
gridspec_kw={"width_ratios": [1, 1.3]})
ii, jj = np.meshgrid(np.arange(nx), np.arange(nx), indexing="ij")
# Broadcasting/`np.hypot`: Abstand jedes Punkts vom Zentrum
rr = np.hypot(ii - m, jj - m)
bild = e * np.sqrt(rr + 1) # 1/√r-Abfall ausgeglichen
w = np.abs(bild[m, m + a:m + b]).max()
ax1.imshow(bild.T, origin="lower", cmap="RdBu_r", vmin=-w, vmax=w)
ax1.set_xlabel("Zelle i"); ax1.set_ylabel("Zelle j")
abstand = np.arange(a, b) * DX
ax2.plot(abstand, strahl_achse, "tab:blue", lw=1.0, label="Achsenstrahl")
ax2.plot(np.arange(int(a / np.sqrt(2)), int(b / np.sqrt(2)))
* np.sqrt(2) * DX, strahl_diag, "tab:orange", lw=1.0,
label="Diagonalstrahl")
ax2.set_xlim(1.2, 1.75) # Zoom: Drift der Kämme
ax2.set_xlabel("Abstand von der Quelle (m)"); ax2.set_ylabel("E (V/m)")
ax2.legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"5 Punkte/λ: Achse v/c = {v_achse:.4f} (1D-Formel: 0.9617)")
print(f" Diagonale v/c = {v_diag:.4f} (Diagonal-Formel: 1.0000)")
5 Punkte/λ: Achse v/c = 0.9613 (1D-Formel: 0.9617)
Diagonale v/c = 1.0001 (Diagonal-Formel: 1.0000)Die Messung bestätigt beide Formeln auf drei bis vier Nachkommastellen — und das Farbbild zeigt, was vier Prozent Richtungsunterschied anrichten: Die Wellenfronten bilden kein Kreis mehr, sondern ein abgerundetes Quadrat auf der Spitze. Damit ist der Themenfaden aus Kapitel 8 zu Ende erzählt: Numerische Dispersion macht kurze Wellen langsam, und in 2D tut sie es richtungsabhängig — ein simulierter Stern funkelt auf groben Gittern viereckig.
Für die Praxis heißt das: In 2D/3D gibt es keinen Zeitschritt, der alle Richtungen exakt macht. Die Faustregeln aus Kapitel 8 bleiben — \(S\) nah an die Grenze, 10–20 Punkte pro Wellenlänge, mehr bei langen Laufstrecken —, aber das bequeme „in 1D löst \(S = 1\) alles” ist Geschichte. Gegen Richtungsfehler hilft nur feiner auflösen (oder klügere Schemata, von denen Kapitel 11 erzählt).
9.4 Offene Ränder, jetzt wirklich: die PML
Bisher enden alle unsere 2D-Gitter in PEC-Wänden: Der Ring läuft hinaus, wird reflektiert und vermüllt die Szene. In 1D half der Mur-Trick (Randpunkt übernimmt den Nachbarwert) — aber er lebte davon, dass die Welle senkrecht und mit bekanntem Tempo ankommt. In 2D trifft sie unter jedem Winkel ein; ein Patentrezept pro Randpunkt gibt es nicht mehr. Kapitel 6 hat den besseren Weg schon vorgezeichnet: eine Rampe aus verlustigem Material, die die Welle schluckt, bevor sie die Wand erreicht — und versprochen, die „volle Konstruktion” in diesem Kapitel nachzuliefern. Wir lösen in zwei Stufen ein und messen jede: Ein Ring läuft in den Rand, danach notieren wir, was im Innenbereich an Feld übrig bleibt (relativ zur PEC-Wand als schlechtestem Fall).
Stufe 1, die \(\sigma\)-Rampe aus Kapitel 6: ein 20 Zellen breiter Rahmen, in dem die Leitfähigkeit quadratisch anwächst und nur das \(E\)-Update dämpft (die ca/cb-Koeffizienten aus Kapitel 6, jetzt als 2D-Felder). Das hilft — aber es bleibt ein messbares Echo: Ein nur elektrisch dämpfender Rahmen hat einen anderen Wellenwiderstand als das Vakuum davor (die Welle spürt den Materialwechsel), und schräg einfallende Wellen verbringen lange Wege in der flachen Rampenzone.
Stufe 2, die Split-Feld-PML (Bérenger 1994): zwei Ideen auf einmal. Erstens bekommt der Verlust einen magnetischen Zwilling — auch die \(u\)-Updates werden gedämpft, abgestimmt im Verhältnis \(\sigma_m/\mu_0 = \sigma/\varepsilon_0\), sodass der Wellenwiderstand der Schicht exakt dem des Vakuums gleicht (das fehlende Stück aus Kapitel 6). Zweitens wird \(E_z\) buchhalterisch in zwei Teilfelder zerlegt, \(E_z = E_{zx} + E_{zy}\): \(E_{zx}\) sammelt den \(\partial_x u_y\)-Anteil des Wirbels, \(E_{zy}\) den \(\partial_y u_x\)-Anteil — und jedes Teilfeld wird nur von „seinem” \(\sigma\)-Profil gedämpft (\(\sigma_x\) wächst nur zu den \(x\)-Rändern hin, \(\sigma_y\) nur zu den \(y\)-Rändern). Eine schräg einlaufende Welle wird so am \(x\)-Rand genau in dem Anteil gedämpft, der in \(x\)-Richtung läuft, während ihr Tangentialanteil ungestört weiterläuft, statt reflektiert zu werden. Physik ist diese Zerlegung nicht mehr — innerhalb der Schicht gilt nicht mehr Maxwell, sondern ein eigens konstruiertes Gleichungssystem, dessen einziger Daseinszweck das reflexionsfreie Schlucken ist.
# von oben: felder(), schritt(), gauss(), S_2D
def restfeld(rand, nx=241, dicke=20, schritte=400, s=S_2D, smax=0.25):
"""Ring läuft in den Rand; Restfeld im Innern nach dem Durchgang."""
m = nx // 2
r_e = np.arange(nx, dtype=float) # σ-Profil: quadratische
prof_e = (np.maximum(0, (dicke - r_e) / dicke)**2 # Rampe an
+ np.maximum(0, (r_e - (nx - 1 - dicke)) / dicke)**2) # beiden
r_u = np.arange(nx - 1) + 0.5 # ... und an den u-Plätzen
prof_u = (np.maximum(0, (dicke - r_u) / dicke)**2
+ np.maximum(0, (r_u - (nx - 1 - dicke)) / dicke)**2)
def koeff(sigma): # ca/cb wie in Kapitel 6
return (1 - sigma) / (1 + sigma), 1 / (1 + sigma)
if rand == "pml":
ezx, ezy = np.zeros((nx, nx)), np.zeros((nx, nx))
ux, uy = np.zeros((nx, nx - 1)), np.zeros((nx - 1, nx))
# [:, None] (Broadcasting): macht aus dem 1D-Profil eine Spalte,
# die NumPy über alle Spalten/Zeilen des 2D-Felds ausbreitet
cax, cbx = [k[:, None] for k in koeff(smax * prof_e)] # σx(x)
cay, cby = [k[None, :] for k in koeff(smax * prof_e)] # σy(y)
caux, cbux = [k[None, :] for k in koeff(smax * prof_u)]
cauy, cbuy = [k[:, None] for k in koeff(smax * prof_u)]
for n in range(schritte):
e = ezx + ezy # das physikalische Feld
ux = caux * ux - cbux * s * np.diff(e, axis=1)
uy = cauy * uy + cbuy * s * np.diff(e, axis=0)
ezx[1:-1, :] = (cax[1:-1] * ezx[1:-1, :]
+ cbx[1:-1] * s * np.diff(uy, axis=0))
ezy[:, 1:-1] = (cay[:, 1:-1] * ezy[:, 1:-1]
- cby[:, 1:-1] * s * np.diff(ux, axis=1))
ezx[m, m] += 0.5 * gauss(n)
ezy[m, m] += 0.5 * gauss(n)
e = ezx + ezy
else: # "pec" oder "rampe"
e, ux, uy = felder(nx, nx)
sigma = (smax * (prof_e[:, None] + prof_e[None, :])
if rand == "rampe" else np.zeros((nx, nx)))
ca, cb = koeff(sigma)
for n in range(schritte):
ux -= s * np.diff(e, axis=1)
uy += s * np.diff(e, axis=0)
wirbel = (np.diff(uy, axis=0)[:, 1:-1]
- np.diff(ux, axis=1)[1:-1, :])
e[1:-1, 1:-1] = (ca[1:-1, 1:-1] * e[1:-1, 1:-1]
+ cb[1:-1, 1:-1] * s * wirbel)
e[m, m] += gauss(n)
return np.abs(e[60:-60, 60:-60]).max()
reste = {rand: restfeld(rand) for rand in ("pec", "rampe", "pml")}
print("Restfeld im Innern, nachdem der Ring den Rand passiert hat:")
for rand, wert in reste.items():
print(f" {rand:6}: {wert:.2e} ({wert / reste['pec'] * 100:7.3f} % "
f"der PEC-Reflexion)")Restfeld im Innern, nachdem der Ring den Rand passiert hat:
pec : 1.13e-02 (100.000 % der PEC-Reflexion)
rampe : 1.68e-03 ( 14.913 % der PEC-Reflexion)
pml : 4.40e-05 ( 0.391 % der PEC-Reflexion)Die Messleiter: Die \(\sigma\)-Rampe drückt das Echo auf rund 15 % der PEC-Reflexion, die PML auf 0,4 % — noch einmal fast vierzigmal leiser, bei gleicher Rahmendicke. Wer mehr will, dreht an Dicke und \(\sigma\)-Profil; in den Werkzeugen von Teil III ist genau diese Schicht (in modernisierter Form) eingebaut, und du weißt jetzt, was sie tut und woran du eine zu dünn geratene erkennst: am Echo.
In dieser HTML-Fassung läuft die Messleiter als Film — derselbe Ring trifft in drei Welten auf den Rand. Die Farbskala ist dabei eine Lupe: Sie ist auf das PEC-Echo geeicht, der auslaufende Ring selbst übersteuert sie also bewusst. Sieh per Einzelschritt zu, was nach dem Auftreffen übrig bleibt: Im PEC-Panel schwappt der komplette Ring zurück und vermüllt die Szene, die \(\sigma\)-Rampe lässt nur ein blasses Echo zurück — und im PML-Panel bleibt die Bühne leer. Die Titel messen das größte Feld im Innenbereich live mit; am Ende stehen dort die drei Stufen der Messleiter.
Code der Animation (nur in der HTML-Fassung)
# von oben: np, plt, felder(), gauss(), S_2D — Schleifen wie in restfeld()
from matplotlib import animation
from IPython.display import HTML
def lauf_mit_bildern(rand, nx=241, dicke=20, schritte=400, s=S_2D,
smax=0.25, je=12):
"""Die restfeld()-Läufe, nur dass alle 12 Schritte fotografiert wird."""
m = nx // 2
r_e = np.arange(nx, dtype=float)
prof_e = (np.maximum(0, (dicke - r_e) / dicke)**2
+ np.maximum(0, (r_e - (nx - 1 - dicke)) / dicke)**2)
r_u = np.arange(nx - 1) + 0.5
prof_u = (np.maximum(0, (dicke - r_u) / dicke)**2
+ np.maximum(0, (r_u - (nx - 1 - dicke)) / dicke)**2)
def koeff(sigma):
return (1 - sigma) / (1 + sigma), 1 / (1 + sigma)
bilder = []
if rand == "pml":
ezx, ezy = np.zeros((nx, nx)), np.zeros((nx, nx))
ux, uy = np.zeros((nx, nx - 1)), np.zeros((nx - 1, nx))
cax, cbx = [k[:, None] for k in koeff(smax * prof_e)]
cay, cby = [k[None, :] for k in koeff(smax * prof_e)]
caux, cbux = [k[None, :] for k in koeff(smax * prof_u)]
cauy, cbuy = [k[:, None] for k in koeff(smax * prof_u)]
for n in range(schritte):
e = ezx + ezy
ux = caux * ux - cbux * s * np.diff(e, axis=1)
uy = cauy * uy + cbuy * s * np.diff(e, axis=0)
ezx[1:-1, :] = (cax[1:-1] * ezx[1:-1, :]
+ cbx[1:-1] * s * np.diff(uy, axis=0))
ezy[:, 1:-1] = (cay[:, 1:-1] * ezy[:, 1:-1]
- cby[:, 1:-1] * s * np.diff(ux, axis=1))
ezx[m, m] += 0.5 * gauss(n)
ezy[m, m] += 0.5 * gauss(n)
if n % je == 0:
bilder.append((ezx + ezy).copy())
else:
e, ux, uy = felder(nx, nx)
sigma = (smax * (prof_e[:, None] + prof_e[None, :])
if rand == "rampe" else np.zeros((nx, nx)))
ca, cb = koeff(sigma)
for n in range(schritte):
ux -= s * np.diff(e, axis=1)
uy += s * np.diff(e, axis=0)
wirbel = (np.diff(uy, axis=0)[:, 1:-1]
- np.diff(ux, axis=1)[1:-1, :])
e[1:-1, 1:-1] = (ca[1:-1, 1:-1] * e[1:-1, 1:-1]
+ cb[1:-1, 1:-1] * s * wirbel)
e[m, m] += gauss(n)
if n % je == 0:
bilder.append(e.copy())
return bilder
filme_r = {rand: lauf_mit_bildern(rand)
for rand in ("pec", "rampe", "pml")}
w_lupe = np.abs(filme_r["pec"][-1]).max() # Lupe: aufs PEC-Echo geeicht
fig_a, achsen = plt.subplots(1, 3, figsize=(7.4, 2.9))
bilder_a, titel_a = {}, {"pec": "PEC-Wand", "rampe": "σ-Rampe",
"pml": "Split-Feld-PML"}
for ax, rand in zip(achsen, filme_r):
bilder_a[rand] = ax.imshow(filme_r[rand][0].T, origin="lower",
cmap="RdBu_r", vmin=-w_lupe, vmax=w_lupe)
ax.set_xticks([]); ax.set_yticks([])
def zeichne(j):
for ax, rand in zip(achsen, filme_r):
feld = filme_r[rand][j]
bilder_a[rand].set_data(feld.T)
rest = np.abs(feld[60:-60, 60:-60]).max()
ax.set_title(f"{titel_a[rand]} — innen {rest:.1e}", fontsize=9)
fig_a.suptitle(f"Schritt {12 * j} von 400", fontsize=10, y=1.0)
return list(bilder_a.values())
anim = animation.FuncAnimation(fig_a, zeichne, frames=len(filme_r["pec"]),
interval=110)
plt.close(fig_a)
HTML(anim.to_jshtml(default_mode="loop"))9.5 Ausblick: die dritte Dimension
Der Schritt von 2D nach 3D bringt keine neue Idee mehr, nur mehr Buchhaltung: Alle sechs Komponenten leben, \(E\) wohnt auf den Kanten eines Würfelgitters, \(B\) auf den Flächenmitten — die berühmte Yee-Zelle, seit 1966 unverändert das Herz jedes FDTD-Programms. Drei Richtungen addieren drei Sinus-Terme in der Dispersionsrelation, die Grenze rückt auf \(S \le 1/\sqrt{3}\), und Dispersion gibt es entlang dreier Achsen-, Flächen- und Raumdiagonalen-Familien. Der wahre Unterschied ist der Preis: Ein Gitter mit \(700\) Punkten pro Kante hat in 2D eine halbe Million Zellen, in 3D 343 Millionen — und weil mit dem Gitter auch die Schrittzahl wächst, skaliert der Aufwand wie \(N^4\). Ab hier lohnt es sich, die Handarbeit an professionelle Löser zu übergeben. Genau dorthin sind wir unterwegs: Kapitel 10 baut zuvor noch die zwei Messinstrumente, ohne die man keinem Löser trauen sollte — Spektren und Konvergenztests.
9.6 Das Kapitel-Programm
programme/kap09/kap09_yee_2d.py bündelt alle vier Befunde eigenständig und prüft sie mit assert-Schranken: Kamm-Tempo der Ringwelle (\(|v/c - 1| < 0{,}01\)), Schachbrett-Explosion (Wachstum auf \(0{,}02\) an der Spurformel, Korrelation \(< -0{,}9\) in beiden Richtungen, Gegenprobe stabil), Richtungsmessung gegen beide Formeln (\(< 2\cdot 10^{-3}\)) und die Rand-Messleiter (Rampe \(< 25\,\%\) von PEC, PML \(< 5\,\%\) der Rampe).
Roter Faden
Fast alles hier war Wiederverwendung: Der Yee-Versatz aus Kapitel 5 wurde vom Stab zum Stadtplan; Spurkriterium und Dispersionsrelation aus Kapitel 8 lieferten Schachbrett-Faktor und Diagonal-Magie; die \(\sigma\)-Rampe aus Kapitel 6 bekam ihren versprochenen magnetischen Zwilling und wurde zur PML; und die Farb-Wahlregel aus Kapitel 4 hat ihre ersten echten Feldbilder eingelöst. Nach vorn: Kapitel 10 macht aus Zeitsignalen Spektren (der „Chor aus Sinuswellen” wird messbar) und aus Gitterverfeinerung einen Vertrauenstest. Kapitel 12 schickt Wellen schräg auf Materialgrenzen (Brechung), Kapitel 14 durch Doppelspalte — beides auf genau dem Gitter, das du in diesem Kapitel gebaut hast.
Übungen
Ü 9.1 (Verstehen). Begründe ohne Rechnung, warum die Stabilitätsgrenze in 3D bei \(S \le 1/\sqrt{3}\) liegt — einmal über die Dispersionsrelation, einmal über das Rauten-Argument. Was ist die gefährlichste Welle in 3D?
HinweisMusterlösung zu Ü 9.1
Über die Relation: In 3D addieren sich drei \(\sin^2\)-Terme; im schlimmsten Fall (alle drei Sinus gleich 1, also Zickzack in \(x\), \(y\) und \(z\) — ein dreidimensionales Schachbrett, bei dem jeder Wert das Negativ seiner sechs Nachbarn ist) verlangt die linke Seite \(\sin(\omega\Delta t/2) = S\sqrt{3} \le 1\). Über die Raute: Ein Zeitschritt erreicht in 3D nur den Oktaeder \(|\Delta i| + |\Delta j| + |\Delta k| \le 1\); die größte einbeschriebene Kugel hat den Radius \(\Delta x/\sqrt{3}\) (Abstand des Mittelpunkts zur schrägen Dreiecksfläche). Pro Dimension kommt ein Summand unter der Wurzel hinzu — die Grenze ist allgemein \(1/\sqrt{d}\), weil die Raumdiagonale mit jeder Dimension länger wird, der Ein-Schritt-Horizont des Schemas aber nicht.
Ü 9.2 (Verstehen). In 1D behielt der laufende Puls seine volle Höhe. Die Ringwelle wird dagegen nach außen leiser. Mit welchem Gesetz muss ihre Amplitude abfallen? (Denke an die Energie — sie verteilt sich auf den wachsenden Umfang.) Prüfe deine Antwort an den Kammhöhen der Ring-Simulation.
HinweisMusterlösung zu Ü 9.2
Die Energie des Rings verteilt sich auf einen Umfang \(2\pi r\), die Energiedichte fällt also wie \(1/r\) — und weil Energie quadratisch in der Amplitude steckt (Kapitel 4), fällt die Amplitude wie \(1/\sqrt{r}\). Dann muss \(E_\text{Kamm}\cdot\sqrt{r}\) konstant sein:
# von oben: felder(), schritt(), gauss(), S_2D
nx = 801
m = nx // 2
e, ux, uy = felder(nx, nx)
radien, hoehen = [], []
for n in range(521):
e, ux, uy = schritt(e, ux, uy, S_2D)
e[m, m] += gauss(n)
if n >= 180 and n % 60 == 0:
strahl = e[m, m:]
r = np.argmax(strahl)
radien.append(r); hoehen.append(strahl[r])
produkt = np.array(hoehen) * np.sqrt(np.array(radien))
print("Radius (Zellen): ", radien)
print("Kammhöhe · sqrt(Radius):", np.round(produkt, 4))Radius (Zellen): [np.int64(109), np.int64(151), np.int64(193), np.int64(235), np.int64(278), np.int64(320)]
Kammhöhe · sqrt(Radius): [0.0973 0.097 0.0961 0.0944 0.0942 0.0933]Das Produkt bleibt nahezu konstant, während der Radius um den Faktor drei wächst — das \(1/\sqrt{r}\)-Gesetz trägt. (Der leichte Restabfall von wenigen Prozent kommt von der numerischen Dispersion, die den Kamm unterwegs etwas verbreitert.) In 3D verteilt sich die Energie auf Kugelflächen \(4\pi r^2\), dort fällt die Amplitude wie \(1/r\) — das wird in Kapitel 21 zur Reichweitenformel des Funks.
Ü 9.3 (Verändern). Lass die Sinusquelle aus der Richtungsmessung mit nur 4 Punkten pro Wellenlänge laufen (bei \(S = 1/\sqrt{2}\)). Berechne zuerst aus den beiden Formeln \(v/c\) für Achse und Diagonale — und sieh dir dann die Form der Wellenfronten an.
HinweisMusterlösung zu Ü 9.3
# von oben: felder(), schritt(), lambda_aus_strahl(), DX, S_2D
n_lambda = 4
w_dt = 2 * np.pi * S_2D / n_lambda
k_achse = 2 * np.arcsin(np.sin(w_dt / 2) / S_2D)
k_diag = 2 * np.sqrt(2) * np.arcsin(np.sin(w_dt / 2) / (np.sqrt(2) * S_2D))
print(f"Theorie: Achse v/c = {w_dt / (S_2D * k_achse):.4f}, "
f"Diagonale v/c = {w_dt / (S_2D * k_diag):.4f}")
nx = 501
m = nx // 2
anlauf = int(4 * 2 * np.pi / w_dt)
e, ux, uy = felder(nx, nx)
for n in range(330):
e, ux, uy = schritt(e, ux, uy, S_2D)
e[m, m] += min(1.0, n / anlauf)**2 * np.sin(w_dt * n)
fig, ax = plt.subplots(figsize=(4.6, 4.0))
ii, jj = np.meshgrid(np.arange(nx), np.arange(nx), indexing="ij")
bild = e * np.sqrt(np.hypot(ii - m, jj - m) + 1)
w = np.abs(bild[m, m + 50:m + 200]).max()
ax.imshow(bild.T, origin="lower", cmap="RdBu_r", vmin=-w, vmax=w)
ax.set_xlabel("Zelle i"); ax.set_ylabel("Zelle j")
plt.tight_layout(); plt.show()Theorie: Achse v/c = 0.9333, Diagonale v/c = 1.0000
Die Theorie sagt \(v/c = 0{,}933\) entlang der Achse gegen exakt \(1{,}000\) diagonal — fast sieben Prozent Unterschied, und das Bild zeigt die Quittung: Aus dem Ring ist ein deutliches Quadrat auf der Spitze geworden. Vier Punkte pro Wellenlänge sind eben keine Auflösung, sondern eine Karikatur — die 10–20-Punkte-Faustregel aus Kapitel 8 gilt in 2D erst recht.
Ü 9.4 (Übertragen). Baue eine PEC-Wand quer durchs Gitter (eine Zeile, in der \(E_z\) nach jedem Schritt auf null gesetzt wird) und lasse darin einen Spalt von etwa einer halben Wellenlänge offen. Schicke die Ringwelle einer Sinusquelle von einer Seite darauf. Was erwartest du hinter der Wand — Schatten, Strahl oder etwas Drittes?
HinweisMusterlösung zu Ü 9.4
# von oben: felder(), schritt(), S_2D
nx = 401
m = nx // 2
n_lambda = 12
w_dt = 2 * np.pi * S_2D / n_lambda
anlauf = int(4 * 2 * np.pi / w_dt)
wand = 240 # Wandzeile (i-Index)
spalt = slice(m - 4, m + 4) # Öffnung: 8 Zellen ≈ 2/3 λ
e, ux, uy = felder(nx, nx)
for n in range(420):
e, ux, uy = schritt(e, ux, uy, S_2D)
e[140, m] += min(1.0, n / anlauf)**2 * np.sin(w_dt * n)
e[wand, :spalt.start] = 0.0 # PEC-Wand ...
e[wand, spalt.stop:] = 0.0 # ... mit offenem Spalt
fig, ax = plt.subplots(figsize=(4.6, 4.0))
w = 0.025 # harte Sättigung vor der Wand,
ax.imshow(e.T, origin="lower", cmap="RdBu_r", # damit das leise Feld dahinter
vmin=-w, vmax=w) # sichtbar wird
ax.axvline(wand, color="k", lw=1.0)
ax.set_xlabel("Zelle i"); ax.set_ylabel("Zelle j")
plt.tight_layout(); plt.show()
Etwas Drittes: Hinter dem Spalt entsteht eine neue Ringwelle, als säße im Spalt selbst eine Punktquelle — die Wellenfronten fächern in den gesamten Halbraum auf, auch in den „Schatten”, den ein Lichtstrahl-Bild vorhersagen würde. Ein Spalt, der klein gegen die Wellenlänge ist, ist aus Sicht der Welle eine Punktquelle. Genau diese Idee — jede Stelle einer Wellenfront als Quelle einer neuen Elementarwelle — heißt Huygens-Prinzip, und sie wird in Kapitel 14 zum Doppelspalt und zur Interferenz ausgebaut: Dort lassen wir zwei solche Spalte gegeneinander antreten.