Zehn Kapitel lang war in diesem Buch jede Zeile Physik selbst geschrieben: das Gitter, der Leapfrog, die Quellen, die Ränder, der Frequenz-Detektor. In diesem Kapitel schreiben wir die Ringwelle aus Kapitel 9 noch einmal — diesmal in einer Handvoll Zeilen, denn den Löser stellt jetzt Meep, ein professionelles, frei verfügbares FDTD-Programm, das seit zwei Jahrzehnten in Forschung und Industrie im Einsatz ist.
Das ist kein Abschied vom Eigenbau, sondern seine Belohnung. Wer einmal selbst ein Yee-Gitter aufgesetzt hat, betritt Meep nicht als Blackbox, sondern als Werkstatt voller Bekannter: Die Auflösung zählt Punkte pro Wellenlänge (Kapitel 8), die Courant-Zahl steht voreingestellt auf \(0{,}5\) (Kapitel 8 und 9), am Rand wartet eine PML (Kapitel 9), und eingebaute DFT-Monitore rechnen exakt unseren Detektor aus Kapitel 10. Jedes dieser Wiedersehen werden wir nachmessen, nicht nur behaupten — und dabei die Königsdisziplin der Validierungs-Checkliste einlösen: das Querrechnen, dasselbe Problem mit zwei unabhängigen Programmen. Am Ende des Kapitels bekommt außerdem eine Geschichte ihren Schluss, die in Kapitel 10 nur angedeutet wurde: die Yagi-Antenne, die fast falsch gebaut worden wäre.
Nach diesem Kapitel kannst du …
Warum überhaupt fremde Software, wenn der eigene Löser läuft? Weil zwischen unserem Lehr-Code und dem Ernstfall ein weites Feld liegt: echte Materialien (Glas, Gold, dispersive Medien), drei Dimensionen, gekrümmte Geometrien, die ein Treppenstufen-Gitter nur grob trifft (Meep glättet sie mit Subpixel-Mittelung), Parallelisierung auf viele Prozessorkerne, und — nicht zu unterschätzen — zwanzig Jahre Fehlersuche durch tausende Nutzer. All das kauft man mit einem import ein. Der Preis: Man muss dem Programm vertrauen. Genau dafür haben wir in Kapitel 10 die Checkliste gebaut, und genau darum wird dieses Kapitel Meep nicht vorstellen, sondern vermessen.
Meep wird am bequemsten über conda installiert (die Python-Distribution Miniforge genügt); ein einziger Befehl legt eine eigene Umgebung samt NumPy und Matplotlib an:
conda create -n meep -c conda-forge pymeep
conda activate meep
pip install necpp # die NEC-Antennen-Engine (für den Yagi-Abschnitt)Alle Code-Zellen dieses Kapitels (und der folgenden) laufen in dieser Umgebung. Die erste Begegnung:
import os
# Umgebungsvariable setzen, BEVOR Meep lädt: unterdrückt harmlose
# Hardware-Warnungen mancher MPI-Installationen
os.environ["HWLOC_HIDE_ERRORS"] = "2"
import numpy as np
import matplotlib.pyplot as plt
import meep as mp
mp.verbosity(0) # Meeps Fortschrittsmeldungen aus — wir messen selbstUsing MPI version 4.1, 1 processes1Die Zeile, die Meep beim Import ausgibt, ist schon die erste Wiedererkennung: Es meldet seine MPI-Version — das Programm ist von Haus aus dafür gebaut, eine große Simulation auf viele Prozessorkerne zu verteilen. Für unsere 2D-Experimente reicht ein Kern; ab echten 3D-Rechnungen wird daraus ein entscheidender Vorteil.
Bevor die erste Simulation läuft, eine Eigenheit, über die jeder Meep-Neuling stolpert: Meep fragt nie nach Metern oder Hertz. Eine Zelle ist „8 mal 8” — acht was?
Die Antwort steckt in den Maxwell-Gleichungen selbst. Im Vakuum enthalten sie genau eine Naturkonstante, die Geschwindigkeit \(c\) — aber keine Länge, keine Zeit, keine Frequenz. Vergrößert man ein Experiment in allen Längen um einen Faktor und streckt die Zeit um denselben Faktor, erfüllen die gestreckten Felder wieder exakt die Maxwell-Gleichungen: Nichts in ihnen kann den Maßstab bemerkt haben. Das heißt aber: Eine einzige Simulation gilt für alle Maßstäbe gleichzeitig, solange die Geometrie in Wellenlängen gemessen dieselbe bleibt — eine Antenne für 2,4 GHz und ihre tausendfach geschrumpfte Kopie für Licht sind dieselbe Rechnung. Diese Skaleninvarianz ist uns übrigens längst vertraut: Schon unser eigener Update-Code kannte weder \(\Delta x\) noch \(\Delta t\) einzeln, sondern nur ihr Verhältnis in der Courant-Zahl \(S\) — und Kapitel 8 hat eingeschärft, dass für die Genauigkeit allein Punkte pro Wellenlänge zählen.
Meep macht daraus ein Einheitensystem: Es setzt \(c = 1\) und überlässt dir die Wahl einer Bezugslänge \(a\). Alle Längen sind dann Vielfache von \(a\), alle Zeiten Vielfache von \(a/c\) (die Zeit, in der Licht die Strecke \(a\) schafft), alle Frequenzen Vielfache von \(c/a\). Daraus folgt die wichtigste Umrechnungsformel der Meep-Welt — eine Frequenz in Meep-Einheiten ist nichts als ein Längenverhältnis:
\[f_\text{meep} \;=\; \frac{f_\text{phys}}{c/a} \;=\; \frac{a}{\lambda}.\]
Die Handrechnung dazu, an einem konkreten Zahlenpunkt: Eine Quelle schwinge mit \(f_\text{meep} = 0{,}5\). Ihre Wellenlänge ist dann \(\lambda = a / 0{,}5 = 2a\) — zwei Bezugslängen, fertig, mehr sagt die Simulation nicht. Erst deine Wahl von \(a\) macht Physik daraus. Wählst du \(a = 1\ \text{µm}\) (Photonik-Konvention), ist \(\lambda = 2\ \text{µm}\) und \(f_\text{phys} = c/\lambda = (3 \cdot 10^8\ \text{m/s}) / (2 \cdot 10^{-6}\ \text{m}) \approx 1{,}5 \cdot 10^{14}\ \text{Hz}\) — Infrarotlicht. Wählst du \(a = 1\) cm, ist \(\lambda = 2\) cm und \(f_\text{phys} \approx 15\) GHz — Radartechnik. Das Ergebnis ist in beiden Fällen derselbe Lauf; nur das Etikett wechselt:
# von oben: nichts (eigenständige Handrechnung); c als Zahl
c = 299792458.0 # m/s — als Variable, Meep kennt kein SI
f_meep = 0.5
for a, einheit in ((1e-6, "1 µm"), (1e-2, "1 cm")):
lam = a / f_meep # λ = a / f_meep, die Formel von oben
print(f"a = {einheit}: λ = {lam:.2e} m,"
f" f = {c / lam:.3e} Hz")a = 1 µm: λ = 2.00e-06 m, f = 1.499e+14 Hz
a = 1 cm: λ = 2.00e-02 m, f = 1.499e+10 HzEine Zahl, zwei Wirklichkeiten — 150 THz Lichttechnik und 15 GHz Radartechnik aus demselben Lauf. Im restlichen Kapitel arbeiten wir direkt in Meep-Einheiten und sparen uns das Etikett; die Übung Ü 11.2 rechnet zur Kontrolle einmal ein WLAN-Problem um.
Jetzt das erste Experiment — und zwar als Querrechnung im Sinne der Checkliste aus Kapitel 10: dasselbe Problem, zwei unabhängige Programme, und am Ende wird verglichen. Die Frage: Liefert Meep dasselbe Feld wie unser eigener 2D-Code aus Kapitel 9? Die Bühne: ein Quadrat von \(8 \times 8\) Meep-Einheiten, Auflösung 50 Punkte pro Einheit, also \(\Delta x = 0{,}02\) und ein Gitter von \(400 \times 400\) Ez-Punkten; gerechnet wird mit Meeps Standard-Courant-Zahl \(S = 0{,}5\) (sicher unter der 2D-Grenze \(1/\sqrt{2}\) aus Kapitel 9), also \(\Delta t = S \cdot \Delta x = 0{,}01\). Anfangsbedingung: alle Felder null; angestoßen wird mit einer Stromquelle — Meep addiert in jedem Schritt einen vorgegebenen Strom \(J(t)\) auf einen Gitterpunkt, das ist exakt unsere weiche Quelle aus Kapitel 6. Als Zeitverlauf geben wir einen Gauß-Puls vor (Mitte \(t_0 = 0{,}6\), Breite \(\tau = 0{,}18\)). Randbedingung: gar keine — lassen wir die PML weg, schließt Meep die Zelle mit PEC-Metallwänden ab, genau wie unser Code. Gemessen wird ein Schnappschuss des Ez-Feldes bei \(t = 3\), kurz bevor der Ring die Wände erreicht. Erfolgskriterium: Die normierten Felder beider Programme unterscheiden sich höchstens um wenige Prozent — bei zwei Programmen, die dasselbe Verfahren auf demselben Gitter rechnen, darf man streng sein.
Ein Aufbaudetail ist begründungspflichtig: Die Quelle sitzt nicht bei \((0, 0)\), sondern beim krummen \((0{,}01,\ 0{,}01)\). Das ist kein Tippfehler, sondern der Yee-Versatz: Meeps Ez-Punkte liegen um einen halben Gitterabstand neben den Zellkanten, der Ursprung \((0,0)\) liegt also zwischen vier Ez-Punkten. Wir setzen die Quelle exakt auf den nächstgelegenen Ez-Punkt — sonst verschmiert Meep den Strom auf mehrere Nachbarpunkte, und unser Eigenbau (der die Quelle auf genau einen Punkt setzt) würde Äpfel mit halbverschobenen Birnen vergleichen.
Zuerst die gemeinsame Bühne und unser eigenes Programm — die drei np.diff-Zeilen aus Kapitel 9, wörtlich übernommen:
# von oben: nichts (Bühne wird hier definiert)
GROESSE = 8.0 # Zellgröße (Meep-Einheiten)
AUFL = 50 # Meeps "resolution": Gitterpunkte pro Einheit
DX = 1.0 / AUFL # Δx = 0,02
S = 0.5 # Meeps Standard-Courant-Zahl
DT = S * DX # Δt = 0,01 (c = 1)
T0, TAU = 0.6, 0.18 # Gauß-Puls der Stromquelle
T_BILD = 3.0 # Schnappschuss-Zeit
NX = int(GROESSE * AUFL) # 400 Ez-Punkte je Richtung
IQ = NX // 2 # Quellindex 200 …
XQ = (IQ + 0.5) * DX - GROESSE / 2 # … liegt bei x = +0,01
def quellpuls(t):
"""Der Strom-Zeitverlauf — für beide Programme derselbe."""
return np.exp(-((t - T0) / TAU) ** 2)
def eigener_lauf(versatz=1.0):
"""Der 2D-Yee-Code aus Kapitel 9 (PEC-Ränder, weiche Quelle).
versatz: zu welcher Zeit (n + versatz)·Δt der Quellstrom im
n-ten Schritt ausgewertet wird — Voreinstellung: die Zeit, auf
der die Felder nach dem Schritt gerade angekommen sind.
"""
e = np.zeros((NX, NX))
ux = np.zeros((NX, NX - 1))
uy = np.zeros((NX - 1, NX))
for n in range(int(round(T_BILD / DT))):
ux -= S * np.diff(e, axis=1)
uy += S * np.diff(e, axis=0)
e[1:-1, 1:-1] += S * (np.diff(uy, axis=0)[:, 1:-1]
- np.diff(ux, axis=1)[1:-1, :])
e[IQ, IQ] -= quellpuls((n + versatz) * DT) * DT # Strom J·Δt
return e
e_eigen = eigener_lauf()
print(f"eigener Lauf fertig: {int(round(T_BILD / DT))} Schritte,"
f" max|Ez| = {np.abs(e_eigen).max():.2e}")eigener Lauf fertig: 300 Schritte, max|Ez| = 1.07e-04(Das Minuszeichen vor der Quelle übernimmt Meeps Vorzeichenwahl: In der Ampère-Gleichung steht der Strom als \(-J\) auf der Seite von \(\partial E/\partial t\). Für den Vergleich ist das nur eine Konvention — aber wer querrechnet, muss Konventionen angleichen.)
Und jetzt dasselbe Experiment in Meep. Die gesamte Physik — Gitter, Update-Gleichungen, Quelle, Ränder — steckt in einem Simulation-Objekt:
# von oben: GROESSE, AUFL, XQ, T0, T_BILD, quellpuls()
quelle = [mp.Source( # mp.Source: eine Feldquelle
mp.CustomSource(quellpuls, # … mit frei wählbarem Zeitverlauf
end_time=10 * T0), # (danach gilt die Quelle als aus)
component=mp.Ez, # … speist die Ez-Komponente (TMz!)
center=mp.Vector3(XQ, XQ))] # mp.Vector3: ein Ort (x, y[, z])
sim = mp.Simulation( # die ganze Simulation als Objekt
cell_size=mp.Vector3(GROESSE, GROESSE), # 8 × 8; kein z → Meep rechnet 2D
resolution=AUFL, # Gitterpunkte pro Längeneinheit
sources=quelle,
boundary_layers=[]) # keine Randschichten → PEC-Wände
sim.run(until=T_BILD) # laufen lassen bis t = 3
# get_array: das Feld als NumPy-Array herausholen (hier: die ganze Zelle)
ez_meep = sim.get_array(component=mp.Ez, center=mp.Vector3(),
size=mp.Vector3(GROESSE, GROESSE))
print(f"Meep fertig: Feld {ez_meep.shape}, max|Ez| = {np.abs(ez_meep).max():.2e}")Meep fertig: Feld (400, 400), max|Ez| = 2.65e-01Kein eigener Zeitschritt, keine diff-Zeilen — alles, was wir in Kapitel 9 gebaut haben, verbirgt sich hinter sim.run(). Aber die Stellschrauben sind dieselben geblieben: resolution ist unser „Punkte pro Wellenlänge”-Regler aus Kapitel 8 (nur auf die Bezugslänge statt auf \(\lambda\) bezogen), die Courant-Zahl steht intern auf \(0{,}5\) und ließe sich mit Courant= verstellen (Ü 11.4 tut das — mit Ansage), und boundary_layers entscheidet über das Schicksal der Wellen am Rand.
Festgelegt? Dann lassen wir die beiden Programme nebeneinander antreten — zuerst entscheiden die Bilder, danach die Zahlen:
# von oben: e_eigen, ez_meep
fig, achsen = plt.subplots(1, 2, figsize=(7.0, 3.6))
for ax, feld, titel in ((achsen[0], e_eigen, "eigener Kapitel-9-Code"),
(achsen[1], ez_meep, "Meep")):
w = np.abs(feld).max() / 3
ax.imshow(feld.T, origin="lower", cmap="RdBu_r", vmin=-w, vmax=w,
extent=(-4, 4, -4, 4))
ax.set_title(titel, fontsize=10)
ax.set_xlabel("x (Meep-Einheiten)")
achsen[0].set_ylabel("y (Meep-Einheiten)")
plt.tight_layout(); plt.show()
Die Bilder sind nicht zu unterscheiden — und auch die erste Kontrollzahl passt: Der Wellenkamm sitzt in beiden Läufen am selben Ort.
# von oben: e_eigen, ez_meep, NX, IQ, DX, XQ, T_BILD, T0
x_achse = (np.arange(NX) + 0.5) * DX - GROESSE / 2 # Meeps Ez-Gitterachse
schnitt_eigen = e_eigen[IQ, :] / np.abs(e_eigen[IQ, :]).max()
schnitt_meep = ez_meep[IQ, :] / np.abs(ez_meep[IQ, :]).max()
rechts = x_achse > XQ + 1.0 # Kammsuche: rechts, außerhalb
r_eigen = x_achse[rechts][np.argmax(np.abs(schnitt_eigen[rechts]))] - XQ
r_meep = x_achse[rechts][np.argmax(np.abs(schnitt_meep[rechts]))] - XQ
print(f"Kammradius: eigener Code {r_eigen:.2f}, Meep {r_meep:.2f}"
f" (zum Vergleich: c·(t − t₀) = {T_BILD - T0:.2f})")Kammradius: eigener Code 2.46, Meep 2.46 (zum Vergleich: c·(t − t₀) = 2.40)„Sieht gleich aus” ist beim Querrechnen wertlos — das hat uns die Refutation in Kapitel 10 gelehrt. Also Zahlen: Wir normieren beide Felder auf ihr Maximum und nehmen die größte punktweise Differenz im gesamten \(400 \times 400\)-Feld.
# von oben: e_eigen, ez_meep
feld_eigen = e_eigen / np.abs(e_eigen).max()
feld_meep = ez_meep / np.abs(ez_meep).max()
print(f"max |Differenz der normierten Felder| ="
f" {np.abs(feld_eigen - feld_meep).max() * 100:.2f} %")max |Differenz der normierten Felder| = 4.73 %Fast fünf Prozent. Für zwei verschiedene Verfahren wäre das ein gutes Ergebnis — für zwei Programme, die dasselbe Verfahren auf demselben Gitter rechnen, ist es viel zu viel. Genau hier trennt sich gutes Querrechnen von Augenwischerei: Eine unerklärte Differenz ist kein Schönheitsfehler, sondern eine Botschaft. Irgendetwas machen die beiden Programme verschieden, und wir sollten wissen, was.
Gehen wir die Verdächtigen durch. Das Gitter? Ist mit dem krummen Quellort absichtlich deckungsgleich gemacht. Die Normierung? Beide Felder sind auf ihr Maximum skaliert. Die Quellform? Punkt für Punkt dieselbe Funktion. Bleibt — die Zeit. In unserem Code steckt eine stillschweigende Entscheidung, die wir in Kapitel 9 nie hinterfragt haben: Zu welchem Zeitpunkt werten wir den Quellstrom aus? Der \(n\)-te Schritt befördert das E-Feld von \(n\,\Delta t\) nach \((n+1)\,\Delta t\); unsere Voreinstellung füttert die Quelle mit der Ankunftszeit \((n+1)\,\Delta t\), genauso gut hätten wir die Startzeit \(n\,\Delta t\) nehmen können. Aber Kapitel 5 hat uns beigebracht, wie der Leapfrog über solche Fragen denkt: Das E-Update integriert über das Intervall von \(n\) bis \(n+1\), und was dabei während des Intervalls wirkt — wie das B-Feld, das auf den halben Schritten lebt —, gehört in die Mitte, auf \((n + \tfrac{1}{2})\,\Delta t\). Wenn Meep den Strom wie ein sauberes Leapfrog-Schema behandelt, müsste genau diese Wahl die Differenz zusammenfallen lassen. Der Test — alle drei Lesarten gegen denselben Meep-Lauf:
# von oben: eigener_lauf(), feld_meep
for versatz, name in ((0.0, "n·Δt (Startzeit)"),
(1.0, "(n+1)·Δt (Ankunftszeit)"),
(0.5, "(n+½)·Δt (Leapfrog-Mitte)")):
feld_v = eigener_lauf(versatz)
feld_v = feld_v / np.abs(feld_v).max()
print(f"Quellzeit {name}: "
f"max|Differenz| = {np.abs(feld_v - feld_meep).max() * 100:.2f} %")Quellzeit n·Δt (Startzeit): max|Differenz| = 4.14 %
Quellzeit (n+1)·Δt (Ankunftszeit): max|Differenz| = 4.73 %
Quellzeit (n+½)·Δt (Leapfrog-Mitte): max|Differenz| = 0.62 %Beide „naheliegenden” Wahlen liegen über vier Prozent daneben — die Leapfrog-Mitte trifft auf \(0{,}62\,\%\). Meep wertet den Quellstrom also am halben Zeitschritt aus, dort, wo er nach der Logik des Yee-Schemas hingehört. Das ist die schönste Sorte Querrechnen-Ergebnis: Die Differenz war kein Fehler, sondern eine Lektion — sie hat uns eine Konvention im Inneren von Meep verraten und nebenbei gezeigt, dass unser eigener Kapitel-9-Code die Quelle einen halben Schritt zu spät gefüttert hat (für alle Messungen dort war das egal: ein gemeinsamer Zeitversatz der ganzen Welle ändert weder Tempo- noch Stabilitätsbefunde).
# von oben: eigener_lauf(), schnitt_meep, x_achse, IQ
e_halb = eigener_lauf(0.5)
schnitt_halb = e_halb[IQ, :] / np.abs(e_halb[IQ, :]).max()
fig, ax = plt.subplots(figsize=(6.6, 3.2))
ax.plot(x_achse, schnitt_halb, lw=1.6, label="eigener Code, Quelle bei (n+½)Δt")
ax.plot(x_achse, schnitt_meep, "--", lw=1.2, label="Meep")
ax.plot(x_achse, (schnitt_halb - schnitt_meep) * 10 - 1.5, lw=0.8,
color="tab:red", label="Differenz × 10 (versetzt)")
ax.axhline(-1.5, color="gray", lw=0.5)
ax.set_xlabel("x (Meep-Einheiten)"); ax.set_ylabel("Ez (normiert)")
ax.legend(fontsize=8, loc="upper right")
plt.tight_layout(); plt.show()
Woher kommt der letzte halbe Prozentpunkt? Ein weiterer Zeitversatz ist es nicht — verschiebt man die Quellzeit testweise in kleinen Schritten, ist \(\tfrac{1}{2}\) exakt das Minimum. Übrig bleiben feinere Unterschiede in der Quellbehandlung (Meep glättet und integriert Quellen etwas anders als unser nackter Ein-Punkt-Strom), und die schlagen genau dort zu Buche, wo das Feld am steilsten ist. Für die Praxis ist das ein Spitzenwert: Zur Einordnung — als später NEC und openEMS dieselbe Antenne rechneten (Kapitel 11.8), galten zehn Prozent Unterschied im Gewinn als Bestätigung.
Jetzt drehen wir den einen Schalter um, der in Kapitel 9 viel Arbeit war: offene Ränder. Bei uns hieß das Split-Feld-PML — zwei Teilfelder, vier Verlustprofile, sorgfältig abgestimmte Koeffizienten. In Meep heißt es boundary_layers=[mp.PML(1.0)]: eine absorbierende Schicht von einer Längeneinheit Dicke, innen an alle Wände gelegt.
Der Versuchsaufbau für die Messung: Frage — wie viel reflektiert Meeps PML? Bühne, Quelle: wie im Ringexperiment. Ablauf: Wir lassen länger laufen, bis \(t = 7\); der Ring (Radius \(c \cdot (t - t_0) = 6{,}4\)) ist dann längst in den Rändern (bei \(4{,}0\)) verschwunden. Gemessen wird das Restfeld im zentralen \(3 \times 3\)-Quadrat, relativ zur Ringamplitude aus dem ersten Lauf. Erfolgskriterium: deutlich unter den \(\sim\!15\,\%\) der σ-Rampe aus Kapitel 9 — die PML soll ihren Namen verdienen. Und weil wir inzwischen misstrauische Messtechniker sind, bauen wir eine Referenz ein: denselben Lauf in einer doppelt so großen Zelle (\(16 \times 16\), ohne PML). Dort sind die Wände so weit weg, dass bis \(t = 7\) keine Reflexion ins Zentrum zurückkehrt — dieses Feld zeigt die reine Physik, und die punktweise Differenz zum PML-Lauf isoliert exakt das, was der Rand verschuldet.
# von oben: GROESSE, AUFL, XQ, T0, NX, quellpuls()
def meep_lauf(groesse, pml, courant=0.5, geometrie=(), bis=3.0):
"""Eine Meep-Simulation wie oben; gibt das Ez-Feld bei t=bis zurück."""
q = [mp.Source(mp.CustomSource(quellpuls, end_time=10 * T0),
component=mp.Ez, center=mp.Vector3(XQ, XQ))]
sim = mp.Simulation(cell_size=mp.Vector3(groesse, groesse),
resolution=AUFL, sources=q, Courant=courant,
geometry=list(geometrie), # Objekte in der Zelle
boundary_layers=[mp.PML(1.0)] if pml else [])
sim.run(until=bis)
return sim.get_array(component=mp.Ez, center=mp.Vector3(),
size=mp.Vector3(groesse, groesse))
ez_pml = meep_lauf(GROESSE, pml=True, bis=7.0)
ez_pec7 = meep_lauf(GROESSE, pml=False, bis=7.0) # PEC zum Kontrast
ez_gross = meep_lauf(2 * GROESSE, pml=False, bis=7.0) # die Referenzbox
ez_ref = ez_gross[NX // 2:NX // 2 + NX, NX // 2:NX // 2 + NX] # innerer Teil
ring_amp = np.abs(ez_meep).max() # Ringamplitude bei t = 3 (von oben)
innen = slice(NX // 2 - 75, NX // 2 + 75) # zentrale 3 × 3
kern = slice(75, NX - 75) # alles außerhalb der PML-Schicht
rest_pml = np.abs(ez_pml[innen, innen]).max()
nachleuchten = np.abs(ez_ref[innen, innen]).max()
reflexion = np.abs(ez_pml[kern, kern] - ez_ref[kern, kern]).max()
print("nach dem Ringdurchgang, relativ zur Ringamplitude:")
print(f" Restfeld im Zentrum (PEC-Wände): "
f"{np.abs(ez_pec7[innen, innen]).max() / ring_amp * 100:5.1f} %")
print(f" Restfeld im Zentrum (Meep-PML) : {rest_pml / ring_amp * 100:.2f} %")
print(f" reine Physik (Referenzbox) : "
f"{nachleuchten / ring_amp * 100:.2f} %")
print(f" echte PML-Reflexion (Differenz): {reflexion / ring_amp:.1e}")nach dem Ringdurchgang, relativ zur Ringamplitude:
Restfeld im Zentrum (PEC-Wände): 55.9 %
Restfeld im Zentrum (Meep-PML) : 0.56 %
reine Physik (Referenzbox) : 0.56 %
echte PML-Reflexion (Differenz): 7.0e-07Drei Lehren auf einmal. Erstens: Mit PEC-Wänden bleibt erwartungsgemäß alles im Kasten — gut die Hälfte der Ringamplitude schwappt als Echo durcheinander. Zweitens: Hinter Meeps PML bleiben \(0{,}56\,\%\) übrig — aber die Referenzbox entlarvt, was da übrig bleibt: Es ist haargenau das 2D-Nachleuchten aus Kapitel 9, die ganz normale Physik einer Linienquelle, kein Randeffekt. Wer dem Vorhersage-Punkt aufgesessen ist und es der PML angekreidet hätte, hat soeben am eigenen Leib erfahren, warum Messmethoden geprüft gehören. Drittens, und das ist die eigentliche Zahl: Die echte Reflexion von Meeps PML — die Differenz zwischen PML-Lauf und Referenz — liegt bei \(7 \cdot 10^{-7}\) der Ringamplitude. Unter einem Millionstel. Unsere selbstgebaute Split-Feld-PML aus Kapitel 9 war stolz auf ihre \(0{,}4\,\%\); Meeps PML (eine weiterentwickelte Variante mit gestreckten Koordinaten und optimiertem Profil) spielt mehrere Ligen darüber. Genau für solche Unterschiede bezahlt man gern mit einem import.
# von oben: ez_pec7, ez_pml, ring_amp, kern
fig, achsen = plt.subplots(1, 3, figsize=(9.6, 3.4))
felder = ((ez_pec7, ring_amp / 3, "PEC-Wände, t = 7"),
(ez_pml, ring_amp / 3, "Meep-PML, t = 7"),
(ez_pml, np.abs(ez_pml[kern, kern]).max(),
"Meep-PML, Farbskala aufs Restfeld"))
for ax, (feld, w, titel) in zip(achsen, felder):
ax.imshow(feld.T, origin="lower", cmap="RdBu_r", vmin=-w, vmax=w,
extent=(-4, 4, -4, 4))
ax.set_title(titel, fontsize=9)
ax.set_xlabel("x (Meep-Einheiten)")
achsen[0].set_ylabel("y (Meep-Einheiten)")
plt.tight_layout(); plt.show()
In Kapitel 10 haben wir Spektren so gemessen: Zeitreihe an einer Sonde aufzeichnen, am Ende rfft. Meep kann das eleganter — ein DFT-Monitor rechnet das Multiplizieren-und-Mitteln unseres Detektors während des Laufs mit, für jede gewünschte Frequenz (deshalb braucht er die Zeitreihe nie zu speichern; bei großen 3D-Feldern ist das der einzige gangbare Weg). Die Behauptung „das ist exakt unser Detektor” ist natürlich messpflichtig.
Der Aufbau: Ringexperiment mit PML, eine Sonde eine Längeneinheit rechts der Quelle. Dort lassen wir beide Messgeräte gleichzeitig laufen — Meeps DFT-Monitor für 14 Testfrequenzen zwischen \(0{,}3\) und \(1{,}6\), und parallel zeichnen wir die Zeitreihe selbst auf und halten hinterher den Kapitel-10-Detektor (\(2\langle s \cdot \sin\rangle\), \(2\langle s \cdot \cos\rangle\), np.hypot) an dieselben Frequenzen. Aufnahmedauer \(T = 12\). Erfolgskriterium: Beide Spektren stimmen bis auf einen konstanten Eichfaktor überein — konstant heißt: dasselbe Verhältnis bei allen 14 Frequenzen, sonst wäre es keine Eichung, sondern eine Abweichung.
# von oben: GROESSE, AUFL, XQ, T0, DT, quellpuls()
T_AUF = 12.0
frequenzen = np.linspace(0.3, 1.6, 14)
q = [mp.Source(mp.CustomSource(quellpuls, end_time=10 * T0),
component=mp.Ez, center=mp.Vector3(XQ, XQ))]
sim = mp.Simulation(cell_size=mp.Vector3(GROESSE, GROESSE), resolution=AUFL,
sources=q, boundary_layers=[mp.PML(1.0)])
sonde_ort = mp.Vector3(XQ + 1.0, XQ)
# add_dft_fields: der eingebaute Detektor — sammelt Ez bei den
# Testfrequenzen auf einem kleinen Fleck um die Sonde
dft = sim.add_dft_fields([mp.Ez], frequenzen.tolist(), center=sonde_ort,
size=mp.Vector3(2 * DX, 2 * DX))
reihe = [] # unsere eigene Zeitreihe
# mp.at_every(Δt, f): ruft f in jedem Zeitschritt auf;
# get_field_point: ein Feldwert an einem Ort
sim.run(mp.at_every(DT, lambda s: reihe.append(
s.get_field_point(mp.Ez, sonde_ort).real)), until=T_AUF)
def dft_wert(k):
"""|Ez(f_k)| am Sondenort (Mitte des Monitor-Flecks)."""
# np.atleast_2d: macht aus einem Skalar/Vektor notfalls eine
# 2D-Matrix — der Monitor-Fleck kann je nach Größe beides liefern
a = np.atleast_2d(sim.get_dft_array(dft, mp.Ez, k))
return np.abs(a[a.shape[0] // 2, a.shape[1] // 2])
spektrum_meep = np.array([dft_wert(k) for k in range(len(frequenzen))])
reihe = np.array(reihe)
t = np.arange(len(reihe)) * DT
spektrum_hand = np.array( # der Detektor aus Kapitel 10
[np.hypot(2 * np.mean(reihe * np.sin(2 * np.pi * f * t)),
2 * np.mean(reihe * np.cos(2 * np.pi * f * t)))
for f in frequenzen])
verhaeltnis = spektrum_meep / spektrum_hand
print(f"Verhältnis Meep-Monitor / Hand-Detektor über 14 Frequenzen:")
print(f" {verhaeltnis.mean():.4f} ± {verhaeltnis.std():.1e}")
print(f" zum Vergleich, T/(2·√(2π)) = "
f"{T_AUF / (2 * np.sqrt(2 * np.pi)):.4f}")Verhältnis Meep-Monitor / Hand-Detektor über 14 Frequenzen:
2.3937 ± 2.7e-15
zum Vergleich, T/(2·√(2π)) = 2.3937Das Verhältnis ist konstant bis auf \(3 \cdot 10^{-15}\) — Maschinengenauigkeit, wie beim rfft-Vergleich in Kapitel 10: Meeps Monitor ist unser Detektor. Und selbst der Eichfaktor lässt sich entzaubern: Meep summiert \(s(t)\,e^{i\omega t}\,\Delta t
/\sqrt{2\pi}\) statt zu mitteln und mit 2 zu multiplizieren — das Verhältnis der beiden Konventionen ist \(T/(2\sqrt{2\pi}) =
2{,}3937\), auf alle vier Stellen genau der gemessene Wert. Wo in Kapitel 10 der Faktor \(2/N\) stand, steht hier eben \(\Delta t /
\sqrt{2\pi}\); Eichungen sind Verabredungen, keine Physik.
Damit ist die Werkstatt-Tour komplett: Quelle, Courant, Gitter, PML, Detektor — wir haben jedes Bauteil von Meep an unserem Eigenbau gemessen und verstehen, was es tut. Zeit für den Blick über Meep hinaus.
„Es muss doch ein bestes Simulationsprogramm geben — das nimmt man dann für alles.”
Warum sie naheliegt: Bei Alltagssoftware funktioniert das ja: Ein guter Browser zeigt jede Webseite, eine gute Tabellenkalkulation rechnet jede Tabelle. Warum sollte ein „guter Feldlöser” nicht jedes Feld lösen?
Was stattdessen stimmt: Hinter jedem Löser steht eine Entscheidung, was diskretisiert wird — das Raumvolumen, nur die Metalloberflächen oder gleich nur Strahlwege. Diese Entscheidung macht ihn für eine Klasse von Problemen brillant und für andere prinzipiell ungeeignet; kein Programmieraufwand der Welt ändert das. Meep löst ein Volumengitter und scheitert darum ehrlich an einem 1 mm dünnen Antennendraht (das Gitter müsste den Draht auflösen — bei 12 cm Wellenlänge ein Auflösungs-Albtraum). NEC diskretisiert nur die Drähte selbst und erledigt dieselbe Antenne in Millisekunden — kann dafür aber mit einem Klotz aus Glas praktisch nichts anfangen. Werkzeugwahl ist Teil der Physik, nicht der Geschmackssache.
Die großen Familien, sortiert nach dem, was sie diskretisieren:
| Familie | rechnet auf | ideal für | schwach bei | freie Software |
|---|---|---|---|---|
| FDTD, Zeitbereich | Volumen + Zeit | Wellen unterwegs; Breitband (ein Puls, alle Frequenzen — Kap. 10); Photonik; dispersive Stoffe | dünne Drähte; el. riesige Objekte; viel Leerraum | Meep; openEMS (FIT, mit Dünndraht-Modell) |
| MoM, Momenten-Methode | Drähte und Leiter-Flächen | Drahtantennen; Abstrahlung im Freiraum; Impedanzen | große Volumen aus Dielektrikum | NEC2 (necpp, 4nec2, xnec2c) |
| FEM, meist Frequenz-Bereich | Volumen, freies Netz | komplexe 3D-Geometrie; Resonatoren; eine Frequenz je Lauf | sehr breitbandig; sehr groß | FreeFEM, Elmer (komm.: HFSS, COMSOL) |
| Strahlen (PO/GO/UTD) | Strahlwege | elektrisch Riesiges (Radarecho, Funkfeld einer Stadt) | Details in Wellenlängen-Größe | (in kommerziellen Hybriden, z. B. FEKO) |
Die Faustregel zum Mitnehmen knüpft an das Leitmotiv aus Kapitel 8 an — es zählen nur Verhältnisse zur Wellenlänge: Strukturen um \(\lambda\) herum (Photonik, Resonatoren, Metamaterialien) → Volumengitter, also FDTD/FEM. Dünne Leiter (\(\text{Radius} \ll \lambda\)) in viel Freiraum → MoM. Objekte \(\gg \lambda\) → Strahlverfahren. Und quer durch alle Familien gilt: scipy ist das Werkzeug zwischen den Werkzeugen — seine dünn besetzten Matrizen und Löser (Kapitel 7), die FFT (Kapitel 10) und seine ODE-Integratoren stecken als Bausteine in halb der Landschaft, und für ein maßgeschneidertes Problem (unser nächstes Kapitel rechnet Brechung wieder zu Fuß) bleiben sie die erste Wahl.
Noch eine Grenze der Landkarte: Sie gilt vom Mikrowellen- bis ins ultraviolette Regime. Bei 50 Hz (Wellenlänge 6000 km!) simuliert niemand Wellenausbreitung — dort regieren Schaltkreis- und Magnetostatik-Werkzeuge; und jenseits des UV übernimmt die Quantenphysik (Kapitel 30 kommt darauf zurück).
Kapitel 10 hat sie als Warnung zitiert, jetzt bekommt sie ihren Schluss — die Geschichte, der dieses Buch seine Querrechnen-Moral verdankt. Sie handelt von einer Yagi-Uda-Antenne (der klassischen Dachantennen-Bauform: ein gespeister Dipol, dahinter ein Reflektorstab, davor mehrere „Direktoren”, die die Abstrahlung bündeln), die für 2,4 GHz entworfen und wirklich gebaut werden sollte.
Akt 1 — das falsche Werkzeug. Der erste Versuch lief in Meep, versteht sich: ein 3D-Modell aus dünnen Metallstäben. Es scheiterte nicht an einem Bug, sondern an der Refutation von oben — 1 mm Drahtradius bei 122 mm Wellenlänge sprengt jedes bezahlbare Volumengitter; die Resonanzlängen kamen je nach Auflösung anders heraus. Der Rückzug auf ein 2D-Modell (Metallstreifen statt Drähte) brachte wunderschöne Bilder: eine saubere Keule, +11 dB Vor/Rück-Verhältnis. Nur: In 2D liegt die Resonanzlänge der Elemente bei \(\sim\!0{,}32\,\lambda\) — ein Artefakt der Flächenwelt. Wer diese Maße in Kupfer gesägt hätte, hätte Schrott gebaut. Das 2D-Modell war qualitativ richtig (so funktioniert eine Yagi!) und quantitativ unbaubar — im Rückblick die wertvollste Einsicht der ganzen Saga.
Akt 2 — die beschämend späte Recherche. Dem 2D-Modell ging wochenlanges Herumprobieren am 3D-Modell voraus — Maße tunen, Auflösung erhöhen, wieder tunen; die Keule zeigte hartnäckig nach hinten. Die Lösung stand längst in der Fachliteratur: Yagis rechnet man in der Lehre als 2D-Modell mit eben jener kurzen Resonanzlänge, und baubare Entwürfe macht man mit der Momentenmethode. Seither gilt im Buch wie im Labor die Regel: Bei etablierten Problemen zuerst die Literatur durchsuchen, dann rechnen. Raten ist die teuerste aller Methoden.
Akt 3 — das richtige Werkzeug. Für dünne Drähte im Freiraum ist die Momentenmethode gemacht: NEC zerlegt jeden Draht in Stromsegmente und löst direkt nach den Strömen — kein Volumengitter, kein Drahtradius-Albtraum (der Radius geht als Parameter in die Gleichungen ein, nicht ins Gitter). Der Versuchsaufbau des folgenden Laufs: fünf parallele Drahtdipole (Reflektor, gespeistes Element, drei Direktoren) mit den Maßen des finalen Entwurfs — Längen um \(0{,}43\)–\(0{,}50\,\lambda\), nicht \(0{,}32\)! —, Drahtradius 1 mm, frei im Raum bei 2,45 GHz. Gemessen werden das Richtdiagramm in der Horizontalebene, der Gewinn vorn und hinten und die Eingangsimpedanz am Speisepunkt (die entscheidet, ob sich die Antenne überhaupt an ein 50-Ω-Kabel anschließen lässt — Kapitel 22 widmet sich genau dieser Frage). Erfolgskriterium: ein Vor/Rück-Verhältnis deutlich über 10 dB und eine Impedanz ohne großen Blindanteil (Resonanz).
# von oben: c (der Zahlenwert aus dem Einheiten-Abschnitt)
import necpp
F_MHZ = 2450.0
LAMBDA = c / (F_MHZ * 1e6) # ≈ 0,122 m
# (Name, Position auf dem Träger in λ, Elementlänge in λ)
YAGI = [("Reflektor", 0.00, 0.500),
("Driven", 0.20, 0.473), # das gespeiste Element
("Direktor 1", 0.35, 0.440),
("Direktor 2", 0.55, 0.440),
("Direktor 3", 0.75, 0.430)]
ctx = necpp.nec_create() # eine NEC-Rechnung anlegen
for i, (_, xpos, laenge) in enumerate(YAGI):
x, halb = xpos * LAMBDA, laenge * LAMBDA / 2
# nec_wire: ein gerader Draht, in 11 Stromsegmente zerlegt
# (Anfangs- und Endpunkt in m, letzter Geometrieparameter: Radius 1 mm)
necpp.nec_wire(ctx, i + 1, 11, x, -halb, 0, x, halb, 0, 0.001, 1.0, 1.0)
necpp.nec_geometry_complete(ctx, 0) # 0 = Freiraum, kein Erdboden
necpp.nec_fr_card(ctx, 0, 1, F_MHZ, 0) # die Rechenfrequenz
necpp.nec_ex_card(ctx, 0, 2, 6, 0, 1.0, 0, 0, 0, 0, 0) # 1 V auf die
# ↑ Draht 2 (Driven), Segment 6 (die Mitte)
necpp.nec_rp_card(ctx, 0, 1, 73, 0, 5, 0, 0, # Richtdiagramm:
90.0, 0.0, 0.0, 5.0, 0.0, 0.0) # Horizont, 5°-Raster
azimut = np.array([necpp.nec_gain(ctx, 0, 0, j) for j in range(73)])
z_re = necpp.nec_impedance_real(ctx, 0)
z_im = necpp.nec_impedance_imag(ctx, 0)
necpp.nec_delete(ctx)
print(f"Gewinn vorwärts : {azimut[0]:5.1f} dBi")
print(f"Gewinn rückwärts: {azimut[36]:5.1f} dBi"
f" → Vor/Rück {azimut[0] - azimut[36]:.1f} dB")
print(f"Eingangsimpedanz: {z_re:.1f} {z_im:+.1f}j Ω")Gewinn vorwärts : 8.9 dBi
Gewinn rückwärts: -5.8 dBi → Vor/Rück 14.6 dB
Eingangsimpedanz: 19.0 +1.2j ΩMillisekunden statt Minuten — und alle Erfolgskriterien erfüllt: fast 9 dBi Gewinn, 14,6 dB Vor/Rück, und der Blindanteil der Impedanz ist praktisch null, die Antenne ist resonant. (Die 19 Ω Wirkanteil sind für Yagis dieser Bauart typisch und verlangen ein Anpassglied ans 50-Ω-Kabel — das Thema von Kapitel 22.)
# von oben: azimut
phi = np.radians(np.arange(73) * 5.0)
fig = plt.figure(figsize=(4.6, 4.6))
ax = fig.add_subplot(projection="polar") # Polardiagramm: Winkel + Radius
ax.plot(phi, np.maximum(azimut, -25)) # Nebenzipfel bei −25 dBi deckeln
ax.set_rmin(-25); ax.set_rmax(10) # radiale Achse: −25 … +10 dBi
ax.set_title("Yagi-Richtdiagramm (dBi), NEC", fontsize=10)
plt.tight_layout(); plt.show()
Akt 4 — das Querrechnen vor dem Bauen. Einem einzigen Programm glaubt man keinen Bauplan. Die Gegenrechnung übernahm openEMS — wieder ein Zeitbereichs-Volumenverfahren wie Meep, aber mit einem Dünndraht-Zusatzmodell, das Drähte unterhalb der Gitterauflösung korrekt behandelt (es läuft in einer eigenen conda-Umgebung und braucht Minuten statt Millisekunden, deshalb zitieren wir hier sein Ergebnis): Gewinn \(\approx 10\) dBi, Vor/Rück \(\approx 12\) dB, Hauptkeule nach vorn. Momentenmethode und Volumen-FDTD — zwei Verfahren, die mathematisch nichts voneinander wissen — stimmen auf etwa ein Dezibel überein. Erst damit war der Entwurf baubar; die Antenne wurde anschließend tatsächlich aus Kupferstäben gebaut.
Die Moral, in drei Zeilen Checkliste: Das Werkzeug muss zur Geometrie passen (Verhältnisse zur Wellenlänge!). Ein vereinfachtes Modell kann qualitativ glänzen und quantitativ unbaubar sein. Und gebaut wird erst, wenn zwei unabhängige Verfahren dasselbe sagen — Querrechnen ist kein Misstrauen gegen Programme, sondern Respekt vor dem Lötkolben.
programme/kap11/kap11_loeser_landschaft.py bündelt die vier Messungen des Kapitels eigenständig und sichert sie mit assert-Schranken: Ringwellen-Querrechnung (naive Quellzeiten \(> 3\,\%\), Leapfrog-Mitte \(< 1\,\%\), Kammradien gleich), PML (Restfeld = Nachleuchten auf \(10\,\%\), echte Reflexion \(< 2 \cdot 10^{-6}\)), DFT-Monitor (Verhältnis konstant auf \(10^{-9}\), Eichfaktor = \(T/(2\sqrt{2\pi})\) auf \(10^{-6}\)) und NEC-Yagi (Gewinn, Vor/Rück und Impedanz in engen Fenstern).
Dieses Kapitel war die Ernte von Teil II: Jedes Bauteil, das wir in den Kapiteln 5–10 selbst geschraubt haben — Yee-Gitter und Leapfrog (Kap. 5), weiche Quellen (Kap. 6), Courant-Zahl und Punkte-pro-\(\lambda\) (Kap. 8), PML und Nachleuchten (Kap. 9), Spektral-Detektor und Querrechnen-Checkliste (Kap. 10) — hängt in Meep am vertrauten Haken, und wir haben jeden Haken einzeln nachgewogen. Sogar die Yagi-Warnung aus Kapitel 10 hat jetzt Anfang, Mitte und versöhnliches Ende. Nach vorn öffnet sich Teil IV: Ab Kapitel 12 (Brechung und Totalreflexion) ist Meep eingearbeitetes Werkzeug statt Untersuchungsobjekt — und die Wahlregel aus Kapitel 4 („Farbe für räumliche Verteilung”) bekommt mit get_array ihren Dauerauftrag. Kapitel 22 löst das 19-Ω-Versprechen der Yagi ein (Leitungen und Anpassung), und Kapitel 24 baut sie — mit Stückliste.
Ü 11.1 (Verstehen). Wähle für jedes Problem das Verfahren (FDTD, MoM, FEM oder Strahlverfahren) und begründe mit einer Zeile: (a) eine UKW-Drahtantenne (145 MHz, 2 mm Drahtdurchmesser) auf dem Hausdach; (b) die Antireflex-Beschichtung einer Solarzelle über das ganze Sonnenspektrum; (c) ein Mikrowellen-Keramikfilter im Metallgehäuse, das bei genau 5,8 GHz arbeiten soll; (d) das Radarecho eines Flugzeugs bei 10 GHz.
(a) MoM (NEC): dünner Draht (\(2\ \text{mm} \ll \lambda \approx 2\ \text{m}\)) in viel Freiraum — das Bilderbuchproblem der Momentenmethode; ein Volumengitter müsste Millimeter auflösen und Meter umfassen. (b) FDTD (Meep): breitbandig (ein Puls, alle Frequenzen, Kapitel 10), dispersive dünne Schichten, Strukturgrößen um die Wellenlänge — Meeps Heimspiel. (c) FEM: eine einzige Arbeitsfrequenz, komplexe 3D-Geometrie aus Dielektrikum und Metall, gesucht ist eine scharfe Resonanz — das ist der Frequenzbereichs-Fall. (openEMS/FDTD ginge auch, müsste aber das lange Ausklingen des Resonators aussitzen.) (d) Strahlverfahren (PO/UTD): das Flugzeug misst hunderte Wellenlängen (\(\lambda = 3\) cm) — jedes Gitter wäre astronomisch; asymptotische Verfahren sind hier nicht Notlösung, sondern die einzige Lösung.
Ü 11.2 (Verstehen). Eine WLAN-Simulation bei \(f = 2{,}45\) GHz soll in Meep laufen; als Bezugslänge wird \(a = 1\) cm gewählt. Berechne von Hand (und prüfe mit drei Zeilen Code): Welche Meep-Frequenz \(f_\text{meep}\) gehört dazu? Und welche resolution brauchst du mindestens, damit das Gitter 20 Punkte pro Wellenlänge hat (Faustregel aus Kapitel 8)?
Die Wellenlänge ist \(\lambda = c/f = 0{,}1224\) m \(= 12{,}24\) cm \(= 12{,}24\,a\). Also \(f_\text{meep} = a/\lambda = 1/12{,}24 =
0{,}0817\). Die Auflösung zählt Punkte pro Bezugslänge; 20 Punkte pro Wellenlänge verlangen \(\text{resolution} = 20/\lambda_\text{meep}
= 20/12{,}24 \approx 1{,}6\), aufgerundet also resolution=2 (mehr schadet nie — Kapitel 8). Dass hier eine Zwei reicht, wo unser Ringexperiment 50 brauchte, liegt nur an der Wahl von \(a\): Die Wellenlänge ist zwölf Bezugslängen groß.
# von oben: c
a = 0.01
lam = c / 2.45e9 # Wellenlänge in m
f_meep = a / lam
print(f"λ = {lam / a:.2f}·a → f_meep = {f_meep:.4f}")
print(f"resolution ≥ {20 * f_meep:.2f} → gewählt: 2")λ = 12.24·a → f_meep = 0.0817
resolution ≥ 1.63 → gewählt: 2Ü 11.3 (Verändern). Stelle der Ringwelle ein Hindernis in den Weg: Fülle mit geometry die rechte Zellenhälfte mit Glas (\(\varepsilon = 4\), also Brechungsindex \(n = 2\) — Kapitel 4). Miss im Schnappschuss bei \(t = 2{,}6\) den Kammradius links (Vakuum) und rechts (Glas). Sage vorher: Welches Verhältnis erwartest du?
Im Glas läuft die Welle mit \(c/n = c/2\) — der rechte Kamm sollte also bei der halben Strecke stehen, Verhältnis 2. Die Messung:
# von oben: meep_lauf(), GROESSE, NX, IQ, DX, x_achse, XQ
# mp.Block: ein Quader; mp.Medium: ein Material (hier ε = 4)
glas = mp.Block(size=mp.Vector3(GROESSE / 2, GROESSE, mp.inf),
center=mp.Vector3(GROESSE / 4, 0),
material=mp.Medium(epsilon=4))
ez_glas = meep_lauf(GROESSE, pml=True, geometrie=[glas], bis=2.6)
schnitt = ez_glas[:, IQ] # Schnitt entlang x durch die Quelle
nach_rechts = schnitt[IQ + 30:] # Kammsuche jenseits der Quelle
nach_links = schnitt[:IQ - 29][::-1]
r_glas = (np.argmax(np.abs(nach_rechts)) + 30) * DX
r_vakuum = (np.argmax(np.abs(nach_links)) + 30) * DX
print(f"Kammradius Vakuum {r_vakuum:.2f}, Glas {r_glas:.2f}"
f" → Verhältnis {r_vakuum / r_glas:.2f} (Theorie: n = 2)")Kammradius Vakuum 2.06, Glas 1.02 → Verhältnis 2.02 (Theorie: n = 2)Gemessen: \(2{,}0\) — und ganz nebenbei steckt im Glas-Lauf schon das nächste Kapitel: An der Grenzfläche wird ein Teil des Rings reflektiert, und die Wellenlänge im Glas ist halbiert. Genau dort machen wir in Kapitel 12 weiter.
Ü 11.4 (Übertragen). Meep lässt dich die Courant-Zahl verstellen. Setze sie auf \(0{,}8\) — das liegt über der 2D-Grenze \(1/\sqrt{2} \approx 0{,}707\) aus Kapitel 9. Sage vorher, was passieren wird und wie das explodierte Feld aussehen muss; prüfe beides nach (die Nachbar-Korrelation aus Kapitel 8/9 überführt den Täter).
Vorhersage: Instabilität — und zwar getragen vom Schachbrett, der Zickzack-Welle in beiden Richtungen zugleich (Kapitel 9); seine Nachbar-Korrelation muss in x und y nahe \(-1\) liegen.
# von oben: meep_lauf(), GROESSE, NX
ez_boom = meep_lauf(GROESSE, pml=False, courant=0.8, bis=3.0)
print(f"max|Ez| nach t = 3: {np.abs(ez_boom).max():.1e}")
b = ez_boom[NX // 2 - 50:NX // 2 + 50, NX // 2 - 50:NX // 2 + 50]
korr_x = np.mean(b[1:, :] * b[:-1, :]) / np.mean(b ** 2)
korr_y = np.mean(b[:, 1:] * b[:, :-1]) / np.mean(b ** 2)
print(f"Nachbar-Korrelation: x {korr_x:+.3f}, y {korr_y:+.3f}")max|Ez| nach t = 3: 1.7e+71
Nachbar-Korrelation: x -0.998, y -0.998Siebzig Größenordnungen in drei Zeiteinheiten, Korrelation \(-0{,}998\) in beiden Richtungen: dasselbe Schachbrett wie in unserem Eigenbau — natürlich, denn es ist dasselbe Verfahren. Profi-Software schützt nicht vor der CFL-Bedingung; sie wählt nur die Voreinstellung (\(0{,}5\)) so, dass man nicht aus Versehen darüber stolpert. (Wer mag, halbiere stattdessen resolution im Ringexperiment und beobachte die Dispersions-Schleppe aus Kapitel 8 — auch sie wandert selbstverständlich mit um.)
Warum FEM-Löser Gleichungssysteme lösen — und Meep einfach marschiert. Tabelle 11.1 sortiert die Verfahren danach, was sie diskretisieren. Quer dazu verläuft die explizit/implizit-Trennlinie aus Kapitel 7, und dass sie verläuft, wie sie verläuft, hat einen benennbaren Grund. Finite Differenzen sind eine Kollokations-Methode: Die Gleichung soll punktweise an den Gitterpunkten gelten, vor der Unbekannten steht nichts — der Zeitschritt ist eine Auswertung, FDTD (und das verwandte FIT von openEMS) marschiert. FEM dagegen baut die Lösung aus überlappenden Ansatzfunktionen über den Netzelementen auf und verlangt die Gleichung nur im gewichteten Mittel (das Prinzip heißt Galerkin). Der Preis steht nach der Raum-Diskretisierung links vor der Zeitableitung: eine sogenannte Massenmatrix \(M\) — die Buchhaltung darüber, wie stark sich benachbarte Ansatzfunktionen überlappen —, sodass die halbdiskrete Gleichung \(M\,\dot u = K\,u\) lautet. Weil \(M\) nicht diagonal ist, kostet schon der explizit gemeinte Zeitschritt pro Schritt die Lösung von \(M\,x = b\): Das „alle Punkte gleichzeitig, Reihenfolge egal” aus Kapitel 7 ist konstruktionsbedingt verbaut. In der Praxis kommt ein zweiter Grund dazu: Die großen kommerziellen FEM-Feldlöser (HFSS, COMSOL) arbeiten meist gleich im Frequenzbereich — eine Frequenz, ein stationäres Problem, und das ist per Definition ein Gleichungssystem; dort gibt es gar keinen Zeitschritt, der explizit sein könnte. Und wer im Zeitbereich ohnehin pro Schritt löst, nimmt meist gleich das unbedingt stabile implizite Verfahren samt großem \(\Delta t\) mit — exakt die Kapitel-7-Abwägung.
Das Schlupfloch: explizites FEM gibt es doch. Zwingend ist das Gleichungssystem nämlich nicht. Beim Mass Lumping nähert man \(M\) durch eine Diagonalmatrix (geschickte Quadratur kann das sogar exakt einrichten) — dann ist \(M^{-1}\) geschenkt und der Schritt wieder eine Auswertung. Und die elegante moderne Lösung heißt Discontinuous Galerkin (DGTD): Die Ansatzfunktionen leben strikt pro Element, ohne Stetigkeitszwang zum Nachbarn; gekoppelt wird über Flüsse an den Elementgrenzen, wie bei Finite-Volumen-Verfahren. Dadurch zerfällt \(M\) in kleine, unabhängige Blöcke je Element, die man einmal vorab invertiert — danach marschiert das Verfahren voll explizit, auf unstrukturierten Netzen und mit hoher Ordnung: FEM-Geometrieflexibilität mit FDTD-Marschverhalten. Der Merksatz: Explizit geht überall dort, wo nach der Raum-Diskretisierung nichts Nichttriviales vor der Unbekannten steht — finite Differenzen haben das gratis, FIT/finite Volumen auch, FEM muss es sich mit Lumping oder DG erarbeiten. Und umgekehrt gibt es implizites FDTD (ADI, Kapitel 7): Die Trennlinie explizit/implizit ist eine Eigenschaft der Matrizen, nicht der Verfahrensfamilie.