Chapter 06

FDTD の限界 — クーラン条件・数値分散・階段近似

この章のゴール.
クーラン条件 $c\Delta t \le \Delta x/\sqrt{3}$ を自分で導出できること。「なぜ小さいセルが 1 つあるだけで計算全体が遅くなるのか」を説明できること。そして数値分散と階段近似が、実務のどの症状として現れるかを知ること。
この章で使う既出の用語（定義は各リンク先）. FDTD（01 章 5 節）、FEM（01 章 5 節）、次元（01 章 5 節）、領域（01 章 5 節）、発散（02 章 2 節）、分散性（03 章 6 節）、中心差分（04 章 2 節）、数値分散（05 章 6 節）、階段近似（05 章 6 節）

1. 安定性 — 破ると一瞬で発散する

FDTD には時間ステップ $\Delta t$ の上限がある。超えると解が指数的に発散し、数十ステップで数値が無限大になる。この条件をクーラン条件（CFL 条件、Courant–Friedrichs–Lewy）と呼ぶ。

導出にはフォン・ノイマンの安定性解析を使う。手順は 3 段だけである。

Step 1: 解を平面波と仮定して代入する

1 次元の更新式（05 章）に、空間で波数 $k$、時間で増幅率 $q$ を持つ解を仮定する:

\[ E^n_i = E_0\, q^n e^{-jki\Delta x}, \qquad H^{n+1/2}_{i+1/2} = H_0\, q^{n+1/2} e^{-jk(i+1/2)\Delta x} \]

「$|q| > 1$ なら毎ステップ増幅されて発散、$|q| \le 1$ なら安定」——これが判定基準になる。

Step 2: 更新式に代入して $q$ の方程式を作る

磁界の更新式に代入すると、$e^{-jk(i+1/2)\Delta x}$ が共通因子として括れて

\[ H_0\left(q^{1/2} - q^{-1/2}\right) = -\frac{\Delta t}{\mu\Delta x}E_0\left(e^{-jk\Delta x/2} - e^{+jk\Delta x/2}\right) = \frac{\Delta t}{\mu\Delta x}E_0 \cdot 2j\sin\!\left(\frac{k\Delta x}{2}\right) \]

電界の更新式からも同様に

\[ E_0\left(q^{1/2} - q^{-1/2}\right) = \frac{\Delta t}{\varepsilon\Delta x}H_0 \cdot 2j\sin\!\left(\frac{k\Delta x}{2}\right) \]

（両式とも $q^{n+1/2}$ で割ってある。）

Step 3: 2 式を掛けて $E_0 H_0$ を消す

\[ \left(q^{1/2} - q^{-1/2}\right)^2 = -\frac{4\Delta t^2}{\mu\varepsilon\,\Delta x^2}\sin^2\!\left(\frac{k\Delta x}{2}\right) \]

$c^2 = 1/(\mu\varepsilon)$ を使い、右辺を $-4S^2\sin^2(k\Delta x/2)$ と書く。ここで

\[ S \equiv \frac{c\,\Delta t}{\Delta x} \qquad\text{（クーラン数）} \]

$\xi = q^{1/2}$ と置くと $\xi^2 - 2A\xi + 1 = 0$（$A = 1 - 2S^2\sin^2(k\Delta x/2)$）の形になり、解は

\[ \xi = A \pm \sqrt{A^2-1} \]

$|A| \le 1$ のとき: $\sqrt{A^2-1}$ は純虚数 → $|\xi|^2 = A^2 + (1-A^2) = 1$ → $|q|=1$、安定（振幅が変わらない）
$|A| > 1$ のとき: 2 解のうち一方が $|\xi|>1$ → 発散

したがって安定条件は $|A|\le 1$、すなわち $2S^2\sin^2(k\Delta x/2) \le 2$。最悪ケース（$\sin^2 = 1$）で

\[ \boxed{\; S = \frac{c\Delta t}{\Delta x} \le 1 \;} \qquad\text{（1 次元）} \]

3 次元では 3 方向の波数が同時に効くので、$\sin^2$ の項が 3 つ足し合わさり

\[ \boxed{\; c\,\Delta t \le \frac{1}{\sqrt{\dfrac{1}{\Delta x^2}+\dfrac{1}{\Delta y^2}+\dfrac{1}{\Delta z^2}}} \;} \]

立方体セル（$\Delta x=\Delta y=\Delta z$）なら $c\Delta t \le \Delta x/\sqrt{3}$ である。

導出は「平面波を代入して増幅率 q を求め、|q| ≤ 1 を課す」だけ。フォン・ノイマンの安定性解析である

物理的な意味は一行で言える.
「1 時間ステップの間に、波が 1 セル分より遠くへ進んではいけない」。
FDTD は隣のセルの値しか見ない。もし波が 1 ステップで 2 セル先まで進むなら、 情報が届いていない場所の値を使って計算することになり、因果関係が壊れる。クーラン数 $S$ は「1 ステップで進む距離 ÷ セル幅」そのものである。

クーラン条件を破ってみる

クーラン数を 1 より大きくすると何が起きるか、実際に走らせて確かめられる。 数十ステップで数値が爆発する様子は、一度見ておくと忘れない。

2. 「最小セル」の呪い

クーラン条件は格子の中で最も小さいセルで決まる。ここに実務上の大問題がある。

細い配線を解像するために 1 箇所だけ 10 µm のセルを入れると、計算領域全体の時間ステップがそのセルに合わせて縮む。

最小セル	$\Delta t$ の上限	10 ns を計算するのに必要なステップ数
1 mm	1.9 ps	約 5,200
100 µm	190 fs	約 52,000
10 µm	19 fs	約 520,000

「ちょっと細かくしただけ」で計算が数日になるのはこの仕組みによる。細分化には常に理由が要る

セルを 1/100 にすると、ステップ数は 100 倍。これに 3 次元のセル数増加（$\times 10^6$）が掛かる。「ちょっと細かくしただけ」で計算が数日になるのはこのためである。

対策は 3 つある。

対策	中身	代償
サブグリッド	細かい領域だけ別格子にして、そこだけ小さい $\Delta t$ で回す	実装が複雑。境界で反射・不安定が起きうる
非一様格子	必要な場所だけ細かい直交格子（急激な変化は避ける）	急に変わると反射が出る。隣接セル比 1.5 以内が目安
陰解法（ADI-FDTD 等。ADI＝交互方向陰解法）	クーラン条件から解放される定式化	1 ステップの計算が重い。数値分散が悪化

実務の鉄則: 「なぜこのセルが小さいのか」を常に問う.
細いスリット、薄い誘電体、小さいギャップ——本当にその寸法を解像する必要があるのか。集中素子で置き換えられないか、等価的な薄板モデルは使えないか。 モデルの簡略化（20 章）が、最も効果の大きい高速化である。

3. 数値分散 — 波の速度が狂う

04 章で「中心差分は実効波数を $\sin(k\Delta x)/\Delta x$ にしてしまう」と見た。 FDTD ではこれが波の伝搬速度の誤差として現れる。

安定性解析で得た関係式（Step 3 の式）を $q = e^{j\omega\Delta t}$ と置いて整理すると、 数値分散関係式が得られる:

\[ \left[\frac{1}{c\Delta t}\sin\!\left(\frac{\omega\Delta t}{2}\right)\right]^2 = \left[\frac{1}{\Delta x}\sin\!\left(\frac{k_x\Delta x}{2}\right)\right]^2 + \left[\frac{1}{\Delta y}\sin\!\left(\frac{k_y\Delta y}{2}\right)\right]^2 + \left[\frac{1}{\Delta z}\sin\!\left(\frac{k_z\Delta z}{2}\right)\right]^2 \]

真の分散関係 $(\omega/c)^2 = k_x^2+k_y^2+k_z^2$ と比べると、 すべての $k$ が $\frac{2}{\Delta}\sin(k\Delta/2)$ に置き換わっている。 $\Delta \to 0$ で真の式に戻るが、有限の $\Delta$ では位相速度がずれる。

対策は結局セル/波長を増やすこと。高次差分の FDTD もあるが境界処理が複雑で商用では一般的でない

数値分散には 2 つの厄介な性質がある。

性質	中身	症状
周波数依存	高い周波数（$k\Delta x$ が大きい）ほど遅くなる	パルスが伝わるうちに崩れて尾を引く
伝搬方向依存（異方性）	格子軸方向と対角方向で速度が違う	円形に広がるはずの波がわずかに歪む

軸方向より対角方向のほうが誤差が小さい.
直感的には「対角は 45° に階段を刻むから悪そう」だが、逆である。対角方向では 3 方向の $\sin$ 項が分散して効くため、誤差が打ち消し合う。 $S = 1/\sqrt{3}$（3D の限界クーラン数）ちょうどで対角方向の分散がゼロになるという性質もある。 クーラン数はギリギリまで大きく取るのが得（安定限界の 0.95〜0.99 倍）—— 計算が速くなるうえに分散も小さい。

対策は結局のところ セル/波長を増やすことに尽きる（04 章の表）。高次差分（4 次精度など）を使う FDTD もあるが、境界処理が複雑になるため商用では一般的でない。

4. 階段近似 — 直交格子の宿命

FDTD は直交格子なので、斜面・曲面は階段状にしか表現できない。

商用ソルバではほぼ標準搭載だが、有効になっているかは確認する価値がある

症状	中身
共振周波数のずれ	円形パッチの実効半径が階段化でずれる → 共振が数 % 動く
偽の散乱	階段の角に電界が集中し、実在しない散乱・損失が出る
収束の遅さ	メッシュを細かくしても誤差が $\Delta x^2$ でなく $\Delta x$ でしか減らない（1 次精度に落ちる）

最後の項目が深刻である。 せっかく中心差分で 2 次精度にしたのに、斜めの境界があるだけで全体の精度が 1 次に落ちてしまう。

対策には適合メッシュ（conformal FDTD）がある。セルが境界で切られている場合に、切られた面積・辺長に応じて更新式の係数を補正する手法で、商用ソルバではほぼ標準搭載である。曲面を扱うときは必ず有効にする。

曲面が主役なら、そもそも FEM を選ぶ.
四面体メッシュ（07 章）なら曲面をそのまま貼れる。「アンテナの放射器が円形」「筐体が丸い」といった問題では、階段近似と戦うより手法を変えるほうが速い——それが 09 章の判断である。

5. 高 Q 共振 — 待ち時間の問題

FDTD はエネルギーが十分減衰するまで計算を続ける必要がある（そうしないと FFT が正しく出ない、16 章）。 Q 値の高い共振器では、減衰に必要な時間が

\[ t_{\text{decay}} \approx \frac{Q}{\pi f} \]

程度になる。$Q = 10{,}000$、$f = 1$ GHz なら 3.2 µs—— $\Delta t = 1$ ps として 300 万ステップ。現実的でない。

対策	中身
信号処理で外挿	減衰の途中で打ち切り、Prony 法・行列ペンシル法（打ち切った波形から減衰する正弦波の重ね合わせを推定し、その先を計算で伸ばす手法）などで残りを予測する（商用ソルバに搭載）
周波数領域に切り替える	高 Q なら FEM（07 章）が圧倒的に有利。共振は周波数領域の得意分野

6. FDTD の得手不得手 — まとめ表

場面	FDTD の適性
広帯域（1 回で全周波数）	◎ 最大の強み
大規模（電気的に大きい）	◎ $O(N)$ でメモリも線形
時間応答・過渡現象（TDR＝時間領域反射測定、16 章／EMI＝電磁妨害、19 章のパルス）	◎ 自然
直方体・層構造が主体（基板、パッケージ）	○
曲面・斜め構造が主役	△ 適合メッシュが要る
高 Q 共振（フィルタ、共振器）	× 待ち時間が非現実的
分散性材料が多数	△ 実装依存。畳み込みが要る
微細構造と大領域が混在	△ 最小セルの呪い

7. この章のまとめ

ポイント	内容
クーラン条件	$c\Delta t \le (\sum 1/\Delta^2)^{-1/2}$。1 ステップで 1 セル以上進めない
導出	平面波を代入 → 増幅率 $q$ → $`	q	\le1`$ の条件
最小セルの呪い	1 箇所の細いセルが全体の $\Delta t$ を支配する
数値分散	高周波ほど遅い・方向で速度が違う。$S$ は限界ギリギリが得
階段近似	曲面で 2 次精度が 1 次に落ちる。適合メッシュを使う
高 Q	減衰待ちが非現実的 → 周波数領域へ