Wyprowadzenie metody FDTD

1. Wprowadzenie

Na tej stronie przedstawiono szczegółowe wyprowadzenie matematyczne podstawowych równań aktualizacji w metodzie różnic skończonych w dziedzinie czasu (FDTD), stosowanych do rozwiązywania równań Maxwella w dziedzinie czasu. Celem tej strony jest inny niż celem głównej strony poświęconej teorii. Nie chodzi o przedstawienie zwięzłego przeglądu, lecz o to, aby wyraźnie pokazać, w jaki sposób ciągłe równania rotacji są przekształcane w dyskretne wyrażenia algebraiczne, które są przetwarzane przez komputer. Jest to również ten element metody FDTD, w którym najłatwiej popełnić błędy w kodzie: błędy w znakach, nieprawidłowe indeksy lub niespójne przesunięcia w czasie mogą natychmiast prowadzić do niestabilnych lub niefizycznych wyników.

W związku z tym rozpoczynamy od prawa Ampère’a i prawa Faradaya w postaci różniczkowej, rozwijamy je na składowe kartezjańskie, a następnie ograniczamy do dwuwymiarowego problemu optycznego. Pochodne przestrzenne zastępujemy skończonymi różnicami centralnymi, pochodne czasowe dyskretyzujemy za pomocą standardowej metody „żaby”, a otrzymane wyrażenia przekształcamy, aby uzyskać jawne wzory aktualizujące składowe pola elektrycznego i magnetycznego.

Aby uprościć prezentację wyprowadzenia, ograniczamy się tutaj do przypadku dwuwymiarowego, w którym \(\partial/\partial y = 0\) i przewodność magnetyczna są ustawione na zero, \(\sigma_m = 0\). Dzięki temu algebra staje się łatwiejsza do zrozumienia, a ponadto wyraźnie widać, jak w praktyce przeprowadzana jest redukcja wymiarowości. W pełnym solverze OghmaNano FDTD równania elektromagnetyczne są jednak rozwiązywane w trzech wymiarach. Dlatego też przedstawione tutaj dwuwymiarowe wyprowadzenie ma służyć jako jasny przykład metody, a nie jako stwierdzenie, że podstawowy solver jest ograniczony do dwóch wymiarów.

Na całej tej stronie wartości składowych pola elektrycznego są obliczane w dyskretnych momentach czasu \(n\Delta t\), natomiast wartości składowych pola magnetycznego są obliczane w momentach czasu odpowiadających połowie jednostki \(\left(n+\tfrac{1}{2}\right)\Delta t\). Jest to standardowa konfiguracja siatki Yee, stosowana w metodzie FDTD. W przeciwieństwie do krótszej strony wprowadzającej, tutaj indeksowanie jest wyraźnie określone, aby można było traktować to jako zarówno instrukcję implementacji, jak i notatkę matematyczną.

2. Równania początkowe

Rozpoczynamy od równań Maxwella w postaci rotora, stosowanych do ośrodka liniowego o przewodności elektrycznej \(\sigma\), przenikalności elektrycznej \(\epsilon\), przenikalności magnetycznej \(\mu\) oraz, dla ogólności, przewodności magnetycznej \(\sigma_m\):

\[ \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}. \]

W postaci macierzowej, loki są

\[ \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}. \]

Oto równania dla przestrzeni i czasu ciągłego. Celem dalszej części strony jest przekształcenie tych równań w wyraźny schemat różnic skończonych, który można wykorzystać do symulacji w czasie.

3. Postać składowa równań opisujących zjawisko wirowania

Rozwinięcie prawa Ampère'a na składowe kartezjańskie daje:

\[ \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}. \]

Rozszerzenie prawa Faradaya prowadzi do:

\[ -\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}. \]

Na tym etapie nie przyjęto jeszcze żadnych założeń dotyczących wymiarowości. Powyższe wyrażenia stanowią jedynie składowe pełnych, trójwymiarowych równań rotacji.

4. Redukcja do dwóch wymiarów

Aby uprościć proces wyprowadzania wzorów, ograniczamy teraz problem do dwóch wymiarów, zakładając, że nie występuje żadna zmiana w kierunku \(y\), a zatem \(\partial / \partial y = 0\). Ten krok jest podejmowany wyłącznie w celu uproszczenia obliczeń i ułatwienia zrozumienia struktury równań aktualizujących. Tę samą procedurę można bezpośrednio rozszerzyć na trzy wymiary, a w praktyce OghmaNano rozwiązuje pełne trójwymiarowe równania Maxwella.

Przy takim założeniu równania pola elektrycznego upraszczają się do:

\[ \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}. \]

Jeśli dodatkowo ustawimy przewodność magnetyczną na zero, \(\sigma_m = 0\), równania opisujące pole magnetyczne przyjmą postać:

\[ \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}. \]

Oto ciągłe równania, które teraz zostaną zdyskretyzowane. Warto wyraźnie zaznaczyć, że jest to redukcja do dwóch wymiarów w odniesieniu do pełnego problemu w trzech wymiarach. Celem nie jest twierdzenie, że wszystkie metody FDTD są z natury dwuwymiarowe, lecz wyprowadzenie prostych równań aktualizacji dla typowego przypadku, w którym jeden kierunek jest pomijany.

5. Siatka Yee, indeksowanie i notacja.

Aby numerycznie rozwiązać powyższe równania, zarówno przestrzeń, jak i czas są dyskretyzowane. Niech odstępy siatki obliczeniowej w kierunkach \(x\) i \(z\) wynoszą odpowiednio \(\Delta x\) i \(\Delta z\), a krok czasowy niech wynosi \(\Delta t\). Do oznaczania kierunków \(x\) i \(z\) używamy całkowitych indeksów \(i\) i \(k\):

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

Czas jest oznaczany liczbą całkowitą \(n\), dzięki czemu

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

W schemacie Yee pola elektryczne i magnetyczne nie są przechowywane w tych samych punktach przestrzeni ani w tym samym momencie czasu. Zamiast tego są one przesunięte względem siebie. Składowe pola elektrycznego są obliczane w określonych momentach czasu \(n\Delta t\), natomiast składowe pola magnetycznego są obliczane w momentach czasu przesuniętych o połowę \(\left(n+\tfrac{1}{2}\right)\Delta t\). To właśnie stąd pochodzi charakterystyczna struktura „skoku” w metodzie FDTD.

Pod względem przestrzennym składowe pola są przesunięte względem siebie w taki sposób, że każda wycentrowana skończono-różnicowa pochodna wektorowa jest naturalnie zdefiniowana pomiędzy sąsiednimi punktami siatki. W pełnej implementacji kodu każda składowa pola ma własną parę indeksów przestrzennych, na przykład \(E_x^n(i,k)\), \(H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2})\) i tak dalej. Aby uprościć zapis, niektóre powtarzające się „poboczne” indeksy nie zawsze są zapisywane w kolejnych krokach. Jednak interpretacja jest zawsze taka sama:

\(i\) odnosi się do kierunku \(x\).
\(k\) odnosi się do kierunku \(z\).
\(n\) indeksuje poziomy czasowe pola elektrycznego wyrażone liczbami całkowitymi, a
\(n+\tfrac{1}{2}\) wskazuje wartości pola magnetycznego w określonych momentach czasu, z dokładnością do półtonu.

Na przykład, centralna aproksymacja różnic skończonych dla \(\partial H_y / \partial z\) ma postać:

\[ \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}, \]

a w podobny sposób przybliżenie \(\partial H_z / \partial x\) względem środka wygląda następująco:

\[ \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}. \]

Po wprowadzeniu powyższego zapisu, przedstawione poniżej wyprowadzenie można jednoznacznie interpretować jako właściwy schemat aktualizacji siatki Yee.

6. Wyprowadzenie równań aktualizacji pola elektrycznego.

Pole elektryczne jest obliczane w dyskretnych momentach czasu \(n\Delta t\), natomiast pole magnetyczne jest obliczane w momentach czasu odpowiadających połowie przedziału \(\left(n+\tfrac{1}{2}\right)\Delta t\). Aby wyprowadzić równania aktualizacji, pochodną czasową przybliża się za pomocą metody różnic skończonych, a człon przewodzący oblicza się, wykorzystując średnią wartości starego i nowego pola elektrycznego. W ten sposób uzyskuje się numerycznie stabilne przybliżenie przewodnictwa za pomocą metody trapezów.

Ten etap uśredniania jest istotny. Gdyby składnik przewodności był traktowany niespójnie, współczynniki aktualizacji nie miałyby już standardowej, tłumionej postaci, a związek z powszechnie stosowanymi równaniami FDTD zostałby zerwany.

6.1 Zaktualizuj równanie dla \(E_x\).

Zacznij od

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

Oblicz wartość współczynnika przewodności, wykorzystując średnią z wartości \(E_x^{n+1}\) i \(E_x^n\), a wartość pochodnej czasowej oblicz, stosując centralną różnicę skończoną:

\[ \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}. \]

Przenieś wyrazy \(E_x^{n+1}\) na lewą stronę:

\[ \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). \]

W związku z tym / Z tego powodu / Dlatego

\[ 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} }. \]

Pomnożenie licznika i mianownika przez \(2\Delta t\) daje nam znaną postać współczynnikową:

\[ 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 Zaktualizuj równanie dla \(E_y\).

Zacznij od

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

Dyskretuj oba wyrażenia zawierające pochodne na siatce Yee:

\[ \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}. \]

Przestawianie rzeczy daje...

\[ 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]. \]

Jest to jedyny składnik pola elektrycznego w dwuwymiarowej redukcji, który zawiera wkład z dwóch oddzielnych pochodnych pola magnetycznego, dlatego najłatwiej popełnić błąd podczas implementacji w kodzie. W szczególności, względny znak między członami pochodnymi \(x\) i \(z\) wynika bezpośrednio z rozwinięcia rotacji i nie powinien być zmieniany.

6.3 Zaktualizuj równanie dla \(E_z\).

Zacznij od

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

Dyskretizacja zapewnia.

\[ \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}. \]

Z tego powodu / Stąd / W związku z tym

\[ 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}. \]

Na tym etapie uzyskano zestawy równań opisujących zmiany pola elektrycznego w dyskretnej formie, uwzględniające tłumienie. Następnym krokiem jest wyprowadzenie odpowiadających im równań pola magnetycznego, obliczanych w punktach czasowych odpowiadających połowie przedziału czasowego.

7. Wyprowadzenie równań aktualizujących pole magnetyczne.

Pole magnetyczne jest aktualizowane w krokach czasowych odpowiadających połowicznym wartościom liczb całkowitych. W przypadku \(\sigma_m = 0\) równania pola magnetycznego nie zawierają członu tłumiącego, dlatego wzory aktualizacji przyjmują prostszą postać iteracyjną. Ta część wyprowadzenia jest algebraicznie krótsza niż w przypadku pola elektrycznego, ponieważ nie występuje w niej człon przewodności trapezoidalnej, który należałoby uwzględnić.

7.1 Zaktualizuj równanie dla \(H_x\).

Zacznij od

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

Dyskretizacja w czasie pomiędzy punktami \(n-\tfrac{1}{2}\) i \(n+\tfrac{1}{2}\) oraz w przestrzeni z wykorzystaniem centralnej różnicy w kierunku \(z\) daje:

\[ \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}. \]

W zwięzłym zapisie półindeksowym, często stosowanym w dowodach, ten sam wynik można zapisać jako:

\[ \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}. \]

W związku z tym / Z tego powodu / Dlatego

\[ 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 Zaktualizuj równanie dla \(H_y\).

Zacznij od

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

Dyskretizacja zapewnia.

\[ \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}. \]

Ponownie, w skróconej notacji indeksów.

\[ \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}. \]

Z tego powodu / Stąd / W związku z tym

\[ 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 Zaktualizuj równanie dla \(H_z\).

Zacznij od

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

Dyskretizacja daje w wyniku

\[ \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}. \]

Równoważnie, w bardziej zwartej notacji:

\[ \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}. \]

W związku z tym / Z tego powodu / Dlatego

\[ 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}. \]

Po wprowadzeniu tych zmian, standardowe równania aktualizacyjne dla metody FDTD, wykorzystywanej w symulacjach propagacji fal, można bezpośrednio zastosować.

8. Ostateczne równania aktualizujące

Po zebraniu wyników, dwuwymiarowe równania aktualizacyjne FDTD przyjmują postać:

\[ 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}. \]

W tej formie przedstawienia schemat „skokowej” metody obliczeniowej jego zasada działania jest jasna. Pole elektryczne w momencie \(n+1\) jest obliczane na podstawie pola magnetycznego w momencie \(n+\tfrac{1}{2}\), natomiast pole magnetyczne w momencie \(n+\tfrac{1}{2}\) jest obliczane na podstawie pola elektrycznego w momencie \(n\). W ten sposób pola zmieniają się naprzemiennie w czasie.

9. Uwagi dotyczące wdrażania.

W praktycznym kodzie współczynniki, przez które mnożone są stare wartości pola oraz wyrazy związane z wirowaniem, są często obliczane z góry i przechowywane, zwłaszcza gdy \(\epsilon\), \(\mu\) lub \(\sigma\) zmieniają się w różnych komórkach. Na przykład, aktualizacja pola elektrycznego zazwyczaj ma następującą postać:

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

gdzie lokalne współczynniki \(C_a\) i \(C_b\) zawierają informacje o parametrach materiału i kroku czasowym. Dzięki temu implementacja jest bardziej przejrzysta i nie trzeba za każdym razem ponownie obliczać tych samych czynników.

Warto również podkreślić, że przedstawione powyżej wyprowadzenie dotyczy jedynie podstawowego schematu aktualizacji FDTD. Aby przeprowadzić użyteczną symulację, potrzebne są dodatkowo źródła, warunki brzegowe, warstwy pochłaniające, takie jak idealnie dopasowane warstwy (PML), oraz krok czasowy spełniający warunek stabilności Couranta. Żadne z tych dodatkowych elementów nie zmienia podstawowej struktury „leapfrog”, którą wyprowadziliśmy, ale każdy z nich modyfikuje sposób, w jaki równania aktualizacji są stosowane w określonych punktach siatki.

Podsumowując, choć powyższy zapis został opracowany z wystarczającą starannością, aby był przydatny przy implementacji, dokładne rozmieszczenie każdego elementu w przestrzeni zależy od przyjętej konwencji siatki w kodzie. Kluczowe jest nie tyle samo w sobie użyte pismo, co spójność: operator zwijający, przesunięcie i znaki muszą odpowiadać temu samemu układowi komórek Yee w całym kodzie.