Herleitung der FDTD-Methode

1. Einleitung

Diese Seite enthält die detaillierte mathematische Ableitung der wichtigsten Aktualisierungsgleichungen für die Finite-Differenzen-Zeitbereichsmethode (FDTD), die zur Lösung der Maxwell-Gleichungen im Zeitbereich verwendet werden. Der Zweck dieser Seite unterscheidet sich von der Hauptseite, die die Theorie behandelt. Hier soll nicht nur eine kurze Übersicht gegeben werden, sondern explizit gezeigt werden, wie die kontinuierlichen Wirbelgleichungen in die diskreten algebraischen Ausdrücke umgewandelt werden, die auf einem Computer verwendet werden. Dies ist auch der Teil der FDTD-Methode, bei dem in der Programmierung am leichtesten Fehler auftreten können: Vorzeichenfehler, falsch platzierte Indizes oder eine inkonsistente Zeitstaffelung führen sofort zu instabilen oder physikalisch unsinnigen Ergebnissen.

Wir beginnen daher mit dem Ampèreschen Gesetz und dem Faradayschen Gesetz in ihrer Differentialform, entwickeln diese in kartesische Komponenten aus und spezialisieren sie anschließend auf ein zweidimensionales optisches Problem. Die räumlichen Ableitungen werden durch zentrierte Finite-Differenzen-Verfahren ersetzt, die zeitlichen Ableitungen werden mit dem Standard-Leapfrog-Verfahren diskretisiert, und die resultierenden Ausdrücke werden so umgeformt, dass explizite Aktualisierungsformeln für die Komponenten des elektrischen und magnetischen Feldes erhalten werden.

Um die Ableitung verständlich zu halten, reduzieren wir die Gleichungen hier auf einen zweidimensionalen Fall, in dem \(\partial/\partial y = 0\) und die magnetische Leitfähigkeit auf Null gesetzt werden, \(\sigma_m = 0\). Dadurch wird die algebraische Darstellung wesentlich einfacher und es wird auch explizit gezeigt, wie die Dimensionsreduktion in der Praxis durchgeführt wird. Im vollständigen OghmaNano FDTD-Solver werden die elektromagnetischen Gleichungen jedoch in drei Dimensionen gelöst. Die hier dargestellte zweidimensionale Ableitung soll daher als ein klares, durchgerechnetes Beispiel für die Methode dienen und nicht als Aussage, dass der zugrunde liegende Solver auf 2D beschränkt ist.

Auf dieser Seite werden die Komponenten des elektrischen Feldes in ganzzahligen Zeitschritten \(n\Delta t\) und die Komponenten des magnetischen Feldes in halbzahligen Zeitschritten \(\left(n+\tfrac{1}{2}\right)\Delta t\) berechnet. Dies ist die Standardanordnung des Yee-Gitters, die in der FDTD-Methode verwendet wird. Im Gegensatz zur kürzeren Übersichtsseite wird die Indizierung hier expliziter angegeben, sodass die Ableitung sowohl als Implementierungsanleitung als auch als mathematische Notiz gelesen werden kann.

2. Ausgangsgleichungen

Wir gehen von Maxwells Wirbelgleichungen in einem linearen Medium mit elektrischer Leitfähigkeit \(\sigma\), elektrischer Permittivität \(\epsilon\), magnetischer Permeabilität \(\mu\) und, um die Allgemeingültigkeit zu gewährleisten, magnetischer Leitfähigkeit \(\sigma_m\) aus:

\[ \nabla \times \boldsymbol{H} = \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t}, \qquad \nabla \times \boldsymbol{E} = - \sigma_m \boldsymbol{H} - \mu \frac{\partial \boldsymbol{H}}{\partial t}. \]

In ihrer bestimmten Form sind die Locken

\[ \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ H_x & H_y & H_z \end{vmatrix}, \qquad \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ E_x & E_y & E_z \end{vmatrix}. \]

Dies sind die Gleichungen für kontinuierliche Zeit und kontinuierlichen Raum. Der Zweck des restlichen Teils dieser Seite besteht darin, sie in ein explizites Finite-Differenzen-Schema umzuwandeln, mit dem man die Lösung schrittweise in der Zeit berechnen kann.

3. Komponentenform der Formel für die Rotation

Wenn man das Ampèresche Gesetz in kartesische Komponenten zerlegt, erhält man:

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = \frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z}, \]

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = \frac{\partial H_x}{\partial z} - \frac{\partial H_z}{\partial x}, \]

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}. \]

Die Erweiterung des Faradayschen Gesetzes ergibt:

\[ -\sigma_m H_x - \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z}, \]

\[ -\sigma_m H_y - \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x}, \]

\[ -\sigma_m H_z - \mu \frac{\partial H_z}{\partial t} = \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y}. \]

In dieser Phase wurde noch keine Annahme über die Dimensionalität getroffen. Die oben genannten Ausdrücke stellen lediglich die Komponentenform der vollständigen, dreidimensionalen Formeln für die Rotation dar.

4. Zweidimensionale Reduktion

Um die Ableitung zu vereinfachen, reduzieren wir das Problem nun auf zwei Dimensionen, indem wir annehmen, dass es keine Variation in der \(y\)-Richtung gibt, also \(\partial / \partial y = 0\). Dieser Schritt wird ausschließlich unternommen, um die algebraischen Ausdrücke zu vereinfachen und die Struktur der Aktualisierungsgleichungen leichter verständlich zu machen. Das gleiche Verfahren kann direkt auf drei Dimensionen erweitert werden, und in der Praxis löst OghmaNano die vollständige 3D-Form der Maxwell-Gleichungen.

Unter dieser Annahme vereinfachen sich die Gleichungen für das elektrische Feld zu:

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}, \]

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}, \]

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x}. \]

Wenn wir die magnetische Leitfähigkeit ebenfalls auf Null setzen, \(\sigma_m = 0\), vereinfachen sich die Gleichungen für das Magnetfeld zu:

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}, \]

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}, \]

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

Dies sind die kontinuierlichen Gleichungen, die nun diskretisiert werden. Es ist wichtig, explizit zu erwähnen, dass es sich hierbei um eine Reduktion des vollständigen dreidimensionalen Problems auf zwei Dimensionen handelt. Ziel ist es nicht, zu behaupten, dass alle FDTD-Methoden grundsätzlich zweidimensional sind, sondern vielmehr, klare Aktualisierungsgleichungen für den häufigen Fall abzuleiten, in dem eine Dimension vernachlässigt wird.

5. Yee-Gitter, Indexierung und Notation

Um die oben genannten Gleichungen numerisch zu lösen, werden sowohl der Raum als auch die Zeit diskretisiert. Dabei sei der Abstand zwischen den Rechengitterpunkten in \(x\)- und \(z\)-Richtung jeweils \(\Delta x\) bzw. \(\Delta z\), und der Zeitschritt sei \(\Delta t\). Für die \(x\)- und \(z\)-Richtung verwenden wir die ganzzahligen Indizes \(i\) bzw. \(k\):

\[ x_i = i \Delta x, \qquad z_k = k \Delta z. \]

Die Zeit wird durch die ganze Zahl \(n\) indiziert, sodass

\[ t^n = n \Delta t. \]

Im Yee-Verfahren werden die elektrischen und magnetischen Felder nicht an denselben Orten im Raum oder zur selben Zeit gespeichert. Stattdessen sind sie versetzt angeordnet. Die Komponenten des elektrischen Feldes werden zu ganzzahligen Zeitpunkten \(n\Delta t\) berechnet, während die Komponenten des magnetischen Feldes zu halbzahligen Zeitpunkten \(\left(n+\tfrac{1}{2}\right)\Delta t\) berechnet werden. Dies ist der Ursprung der „Froschsprung“-Struktur der FDTD-Methode.

Räumlich sind die Feldkomponenten so versetzt angeordnet, dass jede zentrierte, endliche Differenz der Rotation auf natürliche Weise zwischen benachbarten Gitterpunkten definiert ist. In einer vollständigen Code-Implementierung besitzt jede Feldkomponente ihr eigenes Paar von räumlichen Indizes, zum Beispiel \(E_x^n(i,k)\), \(H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2})\) usw. Um die algebraischen Ausdrücke übersichtlich zu halten, werden einige der wiederholten „Beobachter“-Indizes in den folgenden Schritten nicht immer explizit angegeben. Die Interpretation bleibt jedoch immer gleich:

\(i\) indiziert die \(x\)-Richtung.
\(k\) indiziert die \(z\)-Richtung.
\(n\) indiziert die Zeitpunkte des elektrischen Feldes für ganzzahlige Werte, und
\(n+\tfrac{1}{2}\) erfasst die zeitlichen Veränderungen des Magnetfelds in Halbschritt-Intervallen.

Zum Beispiel lautet die zentrierte Finite-Differenzen-Approximation für \(\partial H_y / \partial z\):

\[ \frac{\partial H_y}{\partial z} \approx \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}, \]

und in ähnlicher Weise ist die zentrierte Approximation von \(\partial H_z / \partial x\)

\[ \frac{\partial H_z}{\partial x} \approx \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

Nachdem diese Notation nun eindeutig definiert wurde, kann die folgende Ableitung ohne jeden Zweifel als ein korrektes Verfahren zur Aktualisierung eines Yee-Gitters interpretiert werden.

6. Herleitung der Gleichungen zur Aktualisierung des elektrischen Feldes

Das elektrische Feld wird in ganzzahligen Zeitschritten \(n\Delta t\) berechnet, während das magnetische Feld in halbzähligen Zeitschritten \(\left(n+\tfrac{1}{2}\right)\Delta t\) berechnet wird. Um die Aktualisierungsgleichungen abzuleiten, wird die zeitliche Ableitung durch eine zentrierte Finite-Differenzen-Methode approximiert, und der leitfähige Term wird unter Verwendung des Mittelwerts des alten und neuen elektrischen Feldes berechnet. Dies führt zu einer numerisch stabilen trapezförmigen Behandlung der Leitfähigkeit.

Dieser Mittelungsschritt ist wichtig. Wenn der Leitfähigkeitswert nicht einheitlich behandelt würde, hätten die Aktualisierungskoeffizienten nicht mehr die übliche gedämpfte Form, und die Übereinstimmung mit den üblichen FDTD-Gleichungen wäre verloren.

6.1 Aktualisieren Sie die Gleichung für \(E_x\).

Beginne mit

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}. \]

Berechnen Sie den Leitwert mithilfe des Mittelwerts von \(E_x^{n+1}\) und \(E_x^n\), und berechnen Sie die zeitliche Ableitung mithilfe einer zentrierten Finite-Differenzen-Methode:

\[ \sigma \frac{E_x^{n+1}(i,k) + E_x^n(i,k)}{2} + \epsilon \frac{E_x^{n+1}(i,k) - E_x^n(i,k)}{\Delta t} = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Fassen Sie die \(E_x^{n+1}\)-Begriffe auf der linken Seite zusammen:

\[ \left( \frac{\sigma}{2} + \frac{\epsilon}{\Delta t} \right) E_x^{n+1}(i,k) = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} + \left( \frac{\epsilon}{\Delta t} - \frac{\sigma}{2} \right) E_x^n(i,k). \]

Daher

\[ E_x^{n+1}(i,k) = \frac{ \left( \dfrac{\epsilon}{\Delta t} - \dfrac{\sigma}{2} \right) E_x^n(i,k) - \dfrac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} }{ \dfrac{\epsilon}{\Delta t} + \dfrac{\sigma}{2} }. \]

Wenn man Zähler und Nenner mit \(2\Delta t\) multipliziert, erhält man die bekannte Koeffizientenform:

\[ E_x^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_x^n(i,k) - \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

6.2 Aktualisierte Gleichung für \(E_y\)

Beginne mit

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}. \]

Diskretisiere beide Ableitungsterme auf dem Yee-Gitter:

\[ \sigma \frac{E_y^{n+1}(i,k) + E_y^n(i,k)}{2} + \epsilon \frac{E_y^{n+1}(i,k) - E_y^n(i,k)}{\Delta t} = - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Das Umstellen von Möbeln bringt

\[ E_y^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_y^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \left[ - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} \right]. \]

Dies ist die einzige Komponente des elektrischen Feldes in der 2D-Reduktion, die Beiträge von zwei getrennten Ableitungen des Magnetfeldes enthält, und daher ist sie am anfälligsten für Fehler bei der Programmierung. Insbesondere ergibt sich das relative Vorzeichen zwischen den Ableitungstermen \(x\) und \(z\) direkt aus der Entwicklung der Rotation und sollte nicht verändert werden.

6.3 Aktualisierte Gleichung für \(E_z\)

Beginne mit

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x}. \]

Die Diskretisierung ermöglicht

\[ \sigma \frac{E_z^{n+1}(i,k) + E_z^n(i,k)}{2} + \epsilon \frac{E_z^{n+1}(i,k) - E_z^n(i,k)}{\Delta t} = \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

Daher

\[ E_z^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_z^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

An diesem Punkt wurden die diskreten Gleichungen für die Aktualisierung des elektrischen Feldes in ihrer üblichen, gedämpften Form ermittelt. Der nächste Schritt besteht darin, die entsprechenden Gleichungen für das magnetische Feld abzuleiten, wobei diese für halbzahliges Zeitniveau berechnet werden.

7. Herleitung der Gleichungen zur Aktualisierung des Magnetfelds

Das Magnetfeld wird in halbzähligen Zeitschritten aktualisiert. Bei \(\sigma_m = 0\) enthalten die Gleichungen für das Magnetfeld keinen Dämpfungsterm, sodass die Aktualisierungsformeln die einfachere Vorwärtsform annehmen. Dieser Teil der Ableitung ist algebraisch kürzer als der Fall des elektrischen Feldes, da es keinen trapezförmigen Leitfähigkeitsterm gibt, der berücksichtigt werden muss.

7.1 Aktualisierte Gleichung für \(H_x\)

Beginne mit

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}. \]

Durch eine Diskretisierung in der Zeit zwischen den Halb-Schritten \(n-\tfrac{1}{2}\) und \(n+\tfrac{1}{2}\) sowie im Raum unter Verwendung einer zentrierten Differenz in der \(z\)-Richtung erhält man:

\[ \mu \frac{ H_x^{n+\frac{1}{2}}(i,k+\frac{1}{2}) - H_x^{n-\frac{1}{2}}(i,k+\frac{1}{2}) }{\Delta t} = \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}. \]

In der häufig in Ableitungen verwendeten, kompakten halbindexierten Notation kann dasselbe Ergebnis wie folgt geschrieben werden:

\[ \mu \frac{ H_x^{n+\frac{1}{2}} - H_x^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

Daher

\[ H_x^{n+\frac{1}{2}} = H_x^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

7.2 Aktualisierte Gleichung für \(H_y\)

Beginne mit

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}. \]

Die Diskretisierung ermöglicht

\[ \mu \frac{ H_y^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) - H_y^{n-\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) }{\Delta t} = \frac{ E_z^n(i+1,k+\frac{1}{2}) - E_z^n(i,k+\frac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\frac{1}{2},k+1) - E_x^n(i+\frac{1}{2},k) }{\Delta z}. \]

Wie bereits erwähnt, in der verkürzten halbindexierten Notation:

\[ \mu \frac{ H_y^{n+\frac{1}{2}} - H_y^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

Daher

\[ H_y^{n+\frac{1}{2}} = H_y^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z} \right]. \]

7.3 Aktualisierte Gleichung für \(H_z\)

Beginne mit

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

Die Diskretisierung ermöglicht

\[ \mu \frac{ H_z^{n+\frac{1}{2}}(i+\frac{1}{2},k) - H_z^{n-\frac{1}{2}}(i+\frac{1}{2},k) }{\Delta t} = - \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

Äquivalent dazu, in kompakter Schreibweise:

\[ \mu \frac{ H_z^{n+\frac{1}{2}} - H_z^{n-\frac{1}{2}} }{\Delta t} = - \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]

Daher

\[ H_z^{n+\frac{1}{2}} = H_z^{n-\frac{1}{2}} - \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]

Sobald diese Ersetzungen vorgenommen wurden, ergeben sich die Standard-Aktualisierungsgleichungen für die FDTD-Methode (Finite-Difference Time-Domain) direkt.

8. Abschließende Aktualisierungsgleichungen

Nach der Auswertung der Ergebnisse lauten die zweidimensionalen Aktualisierungsgleichungen der FDTD-Methode:

\[ H_x^{n+\frac{1}{2}}(i,k+\tfrac{1}{2}) = H_x^{n-\frac{1}{2}}(i,k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}, \]

\[ H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) = H_y^{n-\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n(i+1,k+\tfrac{1}{2}) - E_z^n(i,k+\tfrac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\tfrac{1}{2},k+1) - E_x^n(i+\tfrac{1}{2},k) }{\Delta z} \right], \]

\[ H_z^{n+\frac{1}{2}}(i+\tfrac{1}{2},k) = H_z^{n-\frac{1}{2}}(i+\tfrac{1}{2},k) - \frac{\Delta t}{\mu} \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

In dieser Form ausgedrückt, ist die Bedeutung des „Leapfrog“-Verfahrens klar. Das elektrische Feld zum Zeitpunkt \(n+1\) wird aus dem magnetischen Feld zum Zeitpunkt \(n+\tfrac{1}{2}\) berechnet, während das magnetische Feld zum Zeitpunkt \(n+\tfrac{1}{2}\) aus dem elektrischen Feld zum Zeitpunkt \(n\) berechnet wird. Die Felder entwickeln sich also in einer abwechselnden Reihenfolge im Zeitverlauf.

9. Anmerkungen zur Umsetzung

In einer praktischen Implementierung werden die Koeffizienten, mit denen die alten Feldwerte und die Terme, die die Rotation beschreiben, multipliziert werden, oft vorab berechnet und gespeichert, insbesondere wenn \(\epsilon\), \(\mu\) oder \(\sigma\) von Zelle zu Zelle variieren. Beispielsweise erfolgt die Aktualisierung des elektrischen Feldes üblicherweise in folgender Form:

\[ E^{n+1} = C_a E^n + C_b (\nabla \times H), \]

Dabei kodieren die lokalen Koeffizienten \(C_a\) und \(C_b\) die Materialparameter und die Zeitschrittweite. Dies führt zu einer übersichtlicheren Implementierung und vermeidet die Notwendigkeit, die gleichen Faktoren bei jedem Zeitschritt neu zu berechnen.

Es ist auch wichtig zu betonen, dass die oben dargestellte Ableitung nur das grundlegende FDTD-Aktualisierungsschema liefert. Für eine brauchbare Simulation sind zusätzlich Quellen, Randbedingungen, absorbierende Schichten wie perfekt angepasste Schichten (PML) und ein Zeitschritt erforderlich, der die Courant-Stabilitätsbedingung erfüllt. Keine dieser Ergänzungen verändert die hier abgeleitete grundlegende „Froschsprung“-Struktur, aber jede von ihnen modifiziert die Art und Weise, wie die Aktualisierungsgleichungen an bestimmten Gitterpunkten angewendet werden.

Abschließend lässt sich sagen, dass die obige Notation zwar sorgfältig genug formuliert ist, um für die Implementierung relevant zu sein, die genaue räumliche Anordnung der einzelnen Komponenten hängt jedoch von der im Code verwendeten Rasterkonvention ab. Die entscheidende Anforderung ist nicht die Typografie selbst, sondern die Konsistenz: Der Curl-Operator, die Versetzung und die Vorzeichen müssen im gesamten Verlauf der gleichen Yee-Zellenanordnung entsprechen.