Chapter 05

FDTD — マクスウェル方程式を、そのまま時間発展させる

この章のゴール.
マクスウェル方程式から FDTD の更新式を自分の手で導くこと。 Yee 格子で電界と磁界が半セルずれている理由、時間も半ステップずれている理由を、「中心差分を使いたいから」の一言で説明できるようになること。
この章で使う既出の用語（定義は各リンク先）. FDTD（01 章 5 節）、次元（01 章 5 節）、領域（01 章 5 節）、回転（02 章 2 節）、分散（03 章 2 節）、垂直（03 章 5 節）、中心差分（04 章 2 節）

1. 発想 — 交互に更新するだけ

02 章で見た通り、マクスウェル方程式の回転 2 本は互いを餌にしている。

\[ \frac{\partial \mathbf{H}}{\partial t} = -\frac{1}{\mu}\nabla\times\mathbf{E}, \qquad \frac{\partial \mathbf{E}}{\partial t} = \frac{1}{\varepsilon}\left(\nabla\times\mathbf{H} - \sigma\mathbf{E}\right) \]

左辺は「時間変化」、右辺は「空間の傾き」である。つまり

いまの空間分布が分かれば、次の瞬間の値が計算できる。

これを繰り返すだけ。行列も連立方程式も出てこない。これが FDTD（Finite-Difference Time-Domain）である。 1966 年に Kane Yee が提案した。

マクスウェル方程式の回転 2 本が互いを餌にしている構造を、そのまま計算手順にしたのが FDTD である

E を更新 → H を更新 → E を更新 → … と交互に、蛙飛び（leapfrog）のように進む。各ステップの計算は「隣のセルの値を引き算する」だけなので、 巨大な行列を持つ必要がない——これが FDTD が大規模問題に強い根本的な理由である。

2. 1 次元で更新式を導く

まず 1 次元（\(x\) 方向に進む波、\(E_z\) と \(H_y\) だけ）で導出する。マクスウェル方程式の該当成分は

\[ \frac{\partial H_y}{\partial t} = -\frac{1}{\mu}\frac{\partial E_z}{\partial x}, \qquad \frac{\partial E_z}{\partial t} = -\frac{1}{\varepsilon}\frac{\partial H_y}{\partial x} \]

（符号は座標の取り方による。ここでは \(+x\) 方向へ進む波の標準的な取り方に従う。）

両辺を中心差分にしたい。 04 章で見た通り、中心差分は 2 次精度で得だからである。そのためには「微分したい点の両側に値がある」必要がある。そこで:

\(E_z\) を整数格子点 \(i\)、整数時刻 \(n\) に置く
\(H_y\) を半整数格子点 \(i+\frac12\)、半整数時刻 \(n+\frac12\) に置く

×印は紙面奥向きの磁界を表す。この 1 次元の配置を 3 次元に拡張したものが Yee セルである

こう配置すると、\(H_y\) の時間微分を評価したい点 \((i+\frac12,\ n)\) の両側に \(H^{n-1/2}\) と \(H^{n+1/2}\) があり、空間微分を評価したい点の両側に \(E^n_i\) と \(E^n_{i+1}\) がある。 すべてが中心差分になる。

\[ \frac{H^{n+1/2}_{i+1/2} - H^{n-1/2}_{i+1/2}}{\Delta t} = -\frac{1}{\mu}\cdot\frac{E^{n}_{i+1} - E^{n}_{i}}{\Delta x} \]

\(H\) について解けば磁界の更新式:

\[ \boxed{\; H^{n+1/2}_{i+1/2} = H^{n-1/2}_{i+1/2} - \frac{\Delta t}{\mu\,\Delta x}\left(E^{n}_{i+1} - E^{n}_{i}\right) \;} \]

同様に電界の式を点 \((i,\ n+\frac12)\) で中心差分すると電界の更新式:

\[ \boxed{\; E^{n+1}_{i} = E^{n}_{i} - \frac{\Delta t}{\varepsilon\,\Delta x}\left(H^{n+1/2}_{i+1/2} - H^{n+1/2}_{i-1/2}\right) \;} \]

これで FDTD は完成である。 あとはこの 2 行を交互にループで回すだけで、電磁波が計算機の中を伝わる。

式の読み方.
「隣り合う E の差が、その間の H を変化させる」「隣り合う H の差が、その間の E を変化させる」 ——マクスウェル方程式の「回転」が、そのまま「隣との差」になっている。 \(\Delta t/(\mu\Delta x)\) という係数の中に \(\mu, \varepsilon\) が入るので、 材料はセルごとにこの係数を変えるだけで表現できる。誘電体を置くのに特別な処理は要らない——これも FDTD の美点である。

損失（σ）を入れると係数が 2 つになるが、構造は同じ。時間平均を使って無条件安定にする

1D FDTD を走らせる

ブラウザの中で本物の FDTD が動く。 上の 2 行の更新式をそのまま実装したものである。誘電体スラブを置いて反射・透過を見たり、セル数を変えて数値分散（06 章）を体感したりできる。

3. 損失を入れる — 半陰的な扱い

導電率 \(\sigma\) があると電界の式に \(-\sigma E\) が加わる:

\[ \varepsilon\frac{\partial E_z}{\partial t} = -\frac{\partial H_y}{\partial x} - \sigma E_z \]

ここで問題が起きる。\(E_z\) の時間微分を \(n+\frac12\) で中心差分するなら、右辺の \(E_z\) も同じ時刻 \(n+\frac12\) の値でなければ精度が揃わない。しかし \(E\) は整数時刻にしか存在しない。

そこで時間平均で代用する（半陰的な扱い）:

\[ E_z^{n+1/2} \approx \frac{E^{n+1}_z + E^{n}_z}{2} \]

これを代入して \(E^{n+1}\) について解くと

\[ E^{n+1}_i = \underbrace{\frac{1 - \frac{\sigma\Delta t}{2\varepsilon}}{1 + \frac{\sigma\Delta t}{2\varepsilon}}}_{\text{減衰係数 } C_a} E^n_i \;-\; \underbrace{\frac{\frac{\Delta t}{\varepsilon\Delta x}}{1 + \frac{\sigma\Delta t}{2\varepsilon}}}_{C_b}\left(H^{n+1/2}_{i+1/2} - H^{n+1/2}_{i-1/2}\right) \]

なぜ単純に \(E^n\) を使わないのか.
\(-\sigma E^n\) をそのまま使う（陽的な扱い）と、\(\sigma\Delta t/\varepsilon > 2\) で発散する。良導体では \(\sigma\) が \(10^7\) オーダーなので、この条件は簡単に破れる。上の時間平均を使った形は \(\sigma\) がどれだけ大きくても \(|C_a| < 1\) を保つので、無条件に安定である。ソルバの中身にはこういう工夫が随所に入っている。

4. 3 次元へ — Yee セル

3 次元では 6 成分（\(E_x,E_y,E_z,H_x,H_y,H_z\)）を扱う。 Yee の配置は驚くほど自然である。

電界は立方体セルの「辺」の中点に、辺に沿った向きで置く
磁界はセルの「面」の中心に、面に垂直な向きで置く

電界は 1-形式（線に沿って積分する量）、磁界は 2-形式（面を貫く量）——微分幾何の構造が背後にある

こうすると、たとえば \(H_x\) を更新するのに必要な \(\nabla\times\mathbf{E}\) の \(x\) 成分

\[ (\nabla\times\mathbf{E})_x = \frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z} \]

に現れる 4 本の電界（\(E_z\) 2 本、\(E_y\) 2 本）が、 ちょうどその面を囲む 4 辺に載っている。面の周りを 1 周ぶん足すと、その面を貫く磁束の変化になる—— 02 章の積分形（ファラデーの法則）そのものである。

Yee 格子は「偶然うまくいく配置」ではない.
電界を辺に、磁界を面に置くのは、微分幾何の言葉で言えば「電界は 1-形式、磁界は 2-形式」という構造に対応している。だから \(\nabla\cdot\mathbf{B}=0\) が数値的にも厳密に保たれる（磁束が数値誤差で湧き出さない）。適当に差分を当てた手法ではこの保存が壊れ、長時間計算で解が汚れていく。 Yee 格子が 60 年経っても現役なのは、この構造的な正しさによる。

3 次元の更新式は 6 本になるが、形は 1 次元と同じである。\(H_x\) を例に挙げれば

\[ H_x^{n+1/2}\big|_{i,j+\frac12,k+\frac12} = H_x^{n-1/2}\big|_{i,j+\frac12,k+\frac12} + \frac{\Delta t}{\mu}\left[\frac{E_y^{n}\big|_{k+1} - E_y^{n}\big|_{k}}{\Delta z} - \frac{E_z^{n}\big|_{j+1} - E_z^{n}\big|_{j}}{\Delta y}\right] \]

（添字は変化する成分だけ書いた。）残り 5 本は添字を巡回させれば得られる。

5. アルゴリズムの全体像

終了条件はエネルギーの減衰（13 章）。十分減衰する前に切ると FFT の結果が汚れる

段	処理
1	格子・材料係数（\(C_a, C_b\) など）を用意する
2	ループ開始: \(\mathbf{H}\) を半ステップ更新
3	境界条件を適用（PML など、11 章）
4	\(\mathbf{E}\) を 1 ステップ更新
5	励振を加える（ソース、12 章）
6	出力点の値を記録する（時間波形）
7	2 へ戻る。終了条件（エネルギー減衰）まで繰り返す
8	記録した時間波形を FFT して周波数特性へ（16 章）

メモリは場の値だけ、計算は隣接セルとの引き算だけ。 だから

計算量・メモリはセル数に比例（\(O(N)\)）——行列を持つ手法（\(O(N^2)\) 以上）より圧倒的に軽い
隣としかやり取りしないので並列化が容易（GPU との相性が抜群）
1 回の計算で広帯域（16 章）

という強みが出る。大規模・広帯域なら FDTD、という原則の根拠がここにある。

6. FDTD の弱点（次章の予告）

万能ではない。次章で 1 つずつ扱う。

弱点	中身
時間ステップに上限がある	クーラン条件。細かいセルが 1 つでもあると全体が遅くなる
数値分散	波の速度が離散化のせいでわずかに狂う。長距離で効く
階段近似	直交格子なので、斜面や曲面がギザギザになる
鋭い共振が苦手	高 Q だと減衰を待つのに膨大なステップが要る
分散性材料が面倒	\(\varepsilon(\omega)\) を時間領域で扱うには畳み込みが要る

7. この章のまとめ

ポイント	内容
発想	回転 2 本を交互に時間更新するだけ。行列が要らない
Yee 格子	E と H を空間で半セル、時間で半ステップずらす。すべてを中心差分にするため
更新式	「隣の差が、間の場を変える」。材料は係数に入るだけ
損失	時間平均で半陰的に扱い、無条件安定にする
3D	電界は辺、磁界は面。\(\nabla\cdot\mathbf{B}=0\) が構造的に保たれる
強み	\(O(N)\)・並列化容易・1 回で広帯域