ダイポールアンテナやホーンアンテナのように、形状が単純なアンテナは解析解（厳密解）が得られます。しかし、実用的なアンテナ — パッチアンテナの給電プローブ周辺、複雑な形状のスロットアレイ、車体に取り付けたアンテナなど — は、形状が複雑すぎて解析解を求めることができません。
このような問題を解決するのが数値電磁界シミュレーションです。中でもFDTD法（Finite-Difference Time-Domain method）は、マクスウェル方程式を直接、時間ステップごとに逐次計算する手法であり、その概念的なシンプルさと広い適用性から、アンテナ設計で最も広く使われる手法の一つです。
FDTD法の理論を理解することは、以下の場面で大きな力を発揮します。
- アンテナの設計と最適化: 複雑な形状のアンテナの放射パターン、インピーダンス、帯域幅をシミュレーションで予測し、試作前に設計を最適化できます
- EMC（電磁両立性）解析: 電子機器からの不要放射や、外部電磁波による影響を予測し、電磁シールドの設計に活用できます
- 電磁波伝搬解析: 建物内の電波伝搬、人体近傍の電界分布、レーダー断面積（RCS）の計算などに応用できます
本記事の内容
- FDTD法の基本原理（マクスウェル方程式の差分化）
- Yeeセルの空間配置の意味
- 時間ステップの安定条件（CFL条件）の導出
- 吸収境界条件（PML）の原理
- Pythonによる2次元FDTD実装
- 点波源とダイポールのシミュレーション
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
FDTD法の基本原理
マクスウェル方程式の時間発展形
FDTD法の出発点は、源（電流・電荷）を含むマクスウェル方程式の回転形です。
$$ \frac{\partial \bm{H}}{\partial t} = -\frac{1}{\mu}\nabla \times \bm{E} $$
$$ \frac{\partial \bm{E}}{\partial t} = \frac{1}{\epsilon}\nabla \times \bm{H} – \frac{\sigma}{\epsilon}\bm{E} – \frac{1}{\epsilon}\bm{J}_s $$
ここで $\epsilon$ は誘電率、$\mu$ は透磁率、$\sigma$ は導電率、$\bm{J}_s$ は外部電流源です。
これらの式は、$\bm{E}$ の時間変化が $\bm{H}$ の空間微分で決まり、$\bm{H}$ の時間変化が $\bm{E}$ の空間微分で決まるという「相互作用」を表しています。日常的なアナロジーで言えば、2人がキャッチボールをしているようなものです。一方（$\bm{E}$）の状態変化が他方（$\bm{H}$）に伝わり、その変化がまた最初の一方に返ってくる — このプロセスが時間とともに繰り返されることで電磁波が伝搬していきます。
中心差分近似
FDTD法の核心は、時間微分と空間微分を中心差分で近似することです。
関数 $f(x)$ の点 $x_0$ における1階微分を中心差分で近似すると
$$ \frac{\partial f}{\partial x}\bigg|_{x_0} \approx \frac{f(x_0 + \Delta x/2) – f(x_0 – \Delta x/2)}{\Delta x} $$
この近似は2次精度 $O(\Delta x^2)$ を持ちます。つまり、$\Delta x$ を半分にすれば誤差は4分の1になります。前進差分 $[f(x_0+\Delta x) – f(x_0)]/\Delta x$ の1次精度と比べて、格段に高い精度を持ちます。
2次元TM波の場合
2次元問題に話を限定しましょう。$z$ 方向に一様な構造を仮定すると、マクスウェル方程式は2つの独立な偏波モードに分離されます。ここではTMモード（$E_z$, $H_x$, $H_y$ のみ非ゼロ）を扱います。
TMモードのマクスウェル方程式は
$$ \mu \frac{\partial H_x}{\partial t} = -\frac{\partial E_z}{\partial y} $$
$$ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} $$
$$ \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} – \frac{\partial H_x}{\partial y} – \sigma E_z – J_z $$
これらの式に中心差分を適用すると、更新式が得られます。$E_z$ と $\bm{H}$ の更新を時間的に半ステップずらす（リープフロッグ法）ことで、安定かつ高精度な時間発展が実現されます。
中心差分による近似がどのような空間配置を要求するかを、次にYeeセルの概念で見ていきましょう。
Yeeセル — 空間配置の美しさ
Yeeの格子配置
1966年にKane S. Yeeは、マクスウェル方程式の差分化に最適な電磁界成分の空間配置を提案しました。この配置はYeeセル（Yee cell, Yee格子）と呼ばれます。
Yeeセルの核心的なアイデアは、$\bm{E}$ と $\bm{H}$ の各成分を空間的に半格子点（$\Delta x/2$）ずらして配置することです。
2次元TMモードの場合:
- $E_z$ は格子点 $(i, j)$ に配置
- $H_x$ は $y$ 方向に半格子ずれた点 $(i, j+1/2)$ に配置
- $H_y$ は $x$ 方向に半格子ずれた点 $(i+1/2, j)$ に配置
この配置の最大の利点は、回転（curl）演算の中心差分が自然に半格子ずれの位置で定義されることです。たとえば $\partial E_z/\partial x$ は $E_z(i+1,j)$ と $E_z(i,j)$ の差として $(i+1/2, j)$ の位置で評価されますが、これはまさに $H_y$ が定義された位置です。
時間方向のリープフロッグ配置
空間方向と同様に、時間方向でも $\bm{E}$ と $\bm{H}$ を半ステップずらして配置します。
- $E_z$ は時刻 $n\Delta t$ で定義: $E_z^n(i,j)$
- $H_x$, $H_y$ は時刻 $(n+1/2)\Delta t$ で定義: $H_x^{n+1/2}(i,j+1/2)$
この配置により、$\bm{H}$ を更新する式の右辺には $\bm{E}^n$ が、$\bm{E}$ を更新する式の右辺には $\bm{H}^{n+1/2}$ が現れ、すべての微分が中心差分で評価されます。
更新式
上記の配置に基づく具体的な更新式を書き下しましょう。
$H_x$ の更新 ($H_x$ は時刻 $n$ から $n+1/2$ へ、位置 $(i, j+1/2)$ で):
$$ H_x^{n+1/2}(i,j+\tfrac{1}{2}) = H_x^{n-1/2}(i,j+\tfrac{1}{2}) – \frac{\Delta t}{\mu \Delta y}\left[E_z^n(i,j+1) – E_z^n(i,j)\right] $$
$H_y$ の更新 (位置 $(i+1/2, j)$):
$$ H_y^{n+1/2}(i+\tfrac{1}{2},j) = H_y^{n-1/2}(i+\tfrac{1}{2},j) + \frac{\Delta t}{\mu \Delta x}\left[E_z^n(i+1,j) – E_z^n(i,j)\right] $$
$E_z$ の更新 (位置 $(i,j)$):
$$ E_z^{n+1}(i,j) = C_a \cdot E_z^n(i,j) + C_b \left[\frac{H_y^{n+1/2}(i+\tfrac{1}{2},j) – H_y^{n+1/2}(i-\tfrac{1}{2},j)}{\Delta x} – \frac{H_x^{n+1/2}(i,j+\tfrac{1}{2}) – H_x^{n+1/2}(i,j-\tfrac{1}{2})}{\Delta y}\right] $$
ここで導電率 $\sigma$ がある場合の係数は
$$ C_a = \frac{1 – \sigma\Delta t/(2\epsilon)}{1 + \sigma\Delta t/(2\epsilon)}, \quad C_b = \frac{\Delta t/\epsilon}{1 + \sigma\Delta t/(2\epsilon)} $$
自由空間（$\sigma = 0$）では $C_a = 1$, $C_b = \Delta t/\epsilon$ に簡略化されます。
Yeeセルの空間配置と更新式を理解したところで、次に安定条件を導出しましょう。安定条件を満たさないと、シミュレーションが発散してしまいます。
CFL安定条件
なぜ安定条件が必要なのか
FDTD法は陽的（explicit）な時間発展法です。つまり、$n+1$ 時刻の値は $n$ 時刻の値だけから直接計算されます（暗黙的な連立方程式を解く必要がない）。この特長は計算を非常にシンプルにしますが、時間ステップ $\Delta t$ を大きくしすぎると数値解が発散するという代価があります。
CFL条件の導出
安定条件はCFL条件（Courant-Friedrichs-Lewy condition）と呼ばれます。物理的な意味は「数値的な情報伝達速度が、物理的な情報伝達速度（光速）以上でなければならない」ということです。
1次元の場合を考えましょう。光速は $c = 1/\sqrt{\mu\epsilon}$ です。1時間ステップ $\Delta t$ の間に、電磁波は $c\Delta t$ だけ進みます。この距離が格子間隔 $\Delta x$ 以下でなければ、電磁波の伝搬を正しく追跡できません。
$$ c\Delta t \leq \Delta x $$
$$ \Delta t \leq \frac{\Delta x}{c} = \Delta x\sqrt{\mu\epsilon} $$
2次元の場合、$x$ 方向と $y$ 方向の両方を考慮する必要があります。斜め方向（45°）に伝搬する波は、$x$ と $y$ 方向の格子間隔の対角線を超えてはいけません。
$$ c\Delta t \leq \frac{1}{\sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}}} $$
$\Delta x = \Delta y = \Delta$ の場合
$$ \Delta t \leq \frac{\Delta}{c\sqrt{2}} $$
3次元の場合は
$$ \Delta t \leq \frac{\Delta}{c\sqrt{3}} $$
一般に $n$ 次元では
$$ \Delta t \leq \frac{\Delta}{c\sqrt{n}} $$
実用上の設定
実用上は、CFL条件の上限値に安全係数（通常0.9〜0.99）を掛けた値を $\Delta t$ として使います。
$$ \Delta t = \frac{\text{CFL}_n}{c\sqrt{n}\cdot(1/\Delta x^2 + 1/\Delta y^2 + \cdots)^{1/2}} $$
$\text{CFL}_n < 1$ が安全係数です。
空間格子の大きさ $\Delta x$ は、波長 $\lambda$ に対して十分小さく設定する必要があります。一般的な目安は
$$ \Delta x \leq \frac{\lambda_{\min}}{10} \sim \frac{\lambda_{\min}}{20} $$
ここで $\lambda_{\min}$ はシミュレーションで考慮する最短波長です。$\lambda/20$ を使えば、位相誤差は1%以下に抑えられます。
CFL条件を理解したところで、計算領域の境界処理について見ていきましょう。
吸収境界条件 — PML
境界での反射問題
FDTD法の計算領域は有限です。計算領域の端に到達した電磁波は、何も処理しなければ完全に反射されて計算領域内に戻ってきます。アンテナの放射は本来、無限の自由空間に放射されるべきであり、境界からの反射は物理的に存在しないアーティファクトです。
この問題を解決するのが吸収境界条件（ABC: Absorbing Boundary Condition）です。その中でも最も高性能なのがPML（Perfectly Matched Layer）です。
PMLの原理
1994年にBerenger が提案したPML は、「計算領域を取り囲む仮想的な吸収層」です。PMLに入射した電磁波は反射なく層内に進入し、損失によって減衰していきます。
PMLの核心的なアイデアは、界のインピーダンスが計算領域の内部と連続（整合）になるように吸収媒質を設計することです。通常の損失性媒質（$\sigma > 0$）では、電磁波は減衰しますが、媒質の界面で反射が生じます（インピーダンスの不連続）。PMLでは、座標を複素数に拡張（座標伸長）することで、反射のない完全な吸収を実現します。
PMLの実装
PMLの実装で広く使われるのがCPML（Convolutional PML）です。基本的な考え方は、PML層内で導電率を外側に向かって徐々に増加させることです。
PML層内の導電率プロファイルは、典型的には多項式で設定されます。
$$ \sigma(d) = \sigma_{\max}\left(\frac{d}{d_{\text{PML}}}\right)^p $$
ここで $d$ はPML層の内側表面からの距離、$d_{\text{PML}}$ はPML層の厚さ、$p$ は次数（通常 $p = 3$ 〜 $4$）、$\sigma_{\max}$ は最大導電率です。
$\sigma_{\max}$ の最適値は
$$ \sigma_{\max} = \frac{0.8(p+1)}{\eta \Delta} $$
で与えられます（$\eta = \sqrt{\mu/\epsilon}$ は波動インピーダンス、$\Delta$ は格子間隔）。
PML層の厚さは通常10〜20格子点で十分であり、反射レベルは-40 dB〜-80 dBに抑制できます。
PMLの原理を理解したところで、Pythonで2次元FDTD法を実装しましょう。
Pythonによる2次元FDTD実装
2次元TMモードのFDTDを実装し、点波源からの電磁波伝搬をシミュレーションします。
import numpy as np
import matplotlib.pyplot as plt
# --- 2D FDTD シミュレーション (TM mode) ---
# 物理定数
c0 = 3e8 # 光速 [m/s]
epsilon_0 = 8.854e-12 # 真空の誘電率
mu_0 = 4 * np.pi * 1e-7 # 真空の透磁率
# シミュレーションパラメータ
freq = 1e9 # 周波数 1 GHz
wavelength = c0 / freq # 波長 0.3 m
dx = wavelength / 20 # 格子間隔 (λ/20)
dy = dx
Nx = 200 # x方向の格子数
Ny = 200 # y方向の格子数
# CFL条件による時間ステップ
dt = 0.9 * dx / (c0 * np.sqrt(2)) # 安全係数0.9
Nt = 300 # 時間ステップ数
# 電磁界配列の初期化
Ez = np.zeros((Nx, Ny))
Hx = np.zeros((Nx, Ny))
Hy = np.zeros((Nx, Ny))
# 更新係数（自由空間）
Ca = 1.0
Cb = dt / (epsilon_0 * dx)
Da = 1.0
Db = dt / (mu_0 * dx)
# 波源の位置（中心）
src_x, src_y = Nx // 2, Ny // 2
# 簡易吸収境界（Mur 1次ABC）
# 境界の前ステップの値を保存
Ez_prev_x0 = np.zeros(Ny)
Ez_prev_xN = np.zeros(Ny)
Ez_prev_y0 = np.zeros(Nx)
Ez_prev_yN = np.zeros(Nx)
Ez_prev2_x0 = np.zeros(Ny)
Ez_prev2_xN = np.zeros(Ny)
Ez_prev2_y0 = np.zeros(Nx)
Ez_prev2_yN = np.zeros(Nx)
coeff_abc = (c0 * dt - dx) / (c0 * dt + dx)
# スナップショット保存用
snapshots = []
snapshot_times = [50, 100, 150, 250]
# --- メインループ ---
for n in range(Nt):
# 1. Hx の更新: Hx -= Db * dEz/dy
Hx[:, :-1] -= Db * (Ez[:, 1:] - Ez[:, :-1])
# 2. Hy の更新: Hy += Db * dEz/dx
Hy[:-1, :] += Db * (Ez[1:, :] - Ez[:-1, :])
# 3. Ez の更新: Ez += Cb * (dHy/dx - dHx/dy)
Ez[1:, 1:] += Cb * (
(Hy[1:, 1:] - Hy[:-1, 1:]) -
(Hx[1:, 1:] - Hx[1:, :-1])
)
# 4. 波源（ガウシアンパルスの微分）
t_now = n * dt
t0 = 3 / freq # パルス中心
spread = 0.5 / freq # パルス幅
source = -2 * (t_now - t0) / spread**2 * np.exp(
-((t_now - t0) / spread)**2
)
Ez[src_x, src_y] += source
# 5. 簡易吸収境界条件 (Mur 1st order)
Ez[0, :] = Ez_prev_x0 + coeff_abc * (Ez[1, :] - Ez[0, :])
Ez[-1, :] = Ez_prev_xN + coeff_abc * (Ez[-2, :] - Ez[-1, :])
Ez[:, 0] = Ez_prev_y0 + coeff_abc * (Ez[:, 1] - Ez[:, 0])
Ez[:, -1] = Ez_prev_yN + coeff_abc * (Ez[:, -1] - Ez[:, -2])
# 境界値の更新
Ez_prev_x0 = Ez[1, :].copy()
Ez_prev_xN = Ez[-2, :].copy()
Ez_prev_y0 = Ez[:, 1].copy()
Ez_prev_yN = Ez[:, -2].copy()
# スナップショット保存
if n in snapshot_times:
snapshots.append(Ez.copy())
# --- 可視化 ---
fig, axes = plt.subplots(2, 2, figsize=(12, 12))
vmax = max(np.max(np.abs(s)) for s in snapshots) * 0.5
for ax, snap, t_step in zip(axes.flat, snapshots, snapshot_times):
im = ax.imshow(snap.T, origin='lower', cmap='RdBu_r',
vmin=-vmax, vmax=vmax,
extent=[0, Nx*dx*100, 0, Ny*dy*100])
ax.set_title(f'Time step n = {t_step}\n(t = {t_step*dt*1e9:.2f} ns)',
fontsize=12)
ax.set_xlabel('x [cm]', fontsize=11)
ax.set_ylabel('y [cm]', fontsize=11)
ax.plot(src_x*dx*100, src_y*dy*100, 'k+', markersize=10, markeredgewidth=2)
plt.colorbar(im, ax=ax, label='$E_z$ [V/m]', shrink=0.8)
plt.suptitle('2D FDTD Simulation: Point Source Radiation',
fontsize=15, y=1.02)
plt.tight_layout()
plt.show()
FDTDシミュレーションの結果から、電磁波の伝搬特性が視覚的に確認できます。
- 円形の波面が中心から広がる様子が明確に観察できます。点波源からの放射は等方的であり、波面は同心円状に広がります。これは2次元における自由空間の電磁波伝搬の正しい振る舞いです。
- 波長の一致: 波面間の間隔が $\lambda = 30$ cm（0.3 m）に相当することが確認できます。これはFDTDの空間離散化が十分細かく（$\lambda/20$）、数値分散による波長の歪みが小さいことを示しています。
- 吸収境界の効果: 境界付近で波が大きな反射なく吸収されています。Mur 1次ABCでも基本的な反射抑制は機能していますが、高い精度が必要な場合はPMLの実装が推奨されます。
- 時間発展: $n = 50$ では波面がまだ計算領域内に収まっていますが、$n = 250$ では境界に到達し始めています。シミュレーション時間の設定は、対象とする物理現象の時間スケールに合わせて調整する必要があります。
Pythonによるダイポールアンテナのシミュレーション
点波源をz方向の短い電流（微小ダイポール）に拡張し、放射パターンの方向依存性をシミュレーションします。
import numpy as np
import matplotlib.pyplot as plt
# --- 2D FDTD: ダイポール放射 ---
c0 = 3e8
epsilon_0 = 8.854e-12
mu_0 = 4 * np.pi * 1e-7
freq = 1e9
wavelength = c0 / freq
dx = wavelength / 20
dy = dx
Nx = 250
Ny = 250
dt = 0.9 * dx / (c0 * np.sqrt(2))
Nt = 400
Ez = np.zeros((Nx, Ny))
Hx = np.zeros((Nx, Ny))
Hy = np.zeros((Nx, Ny))
Cb = dt / (epsilon_0 * dx)
Db = dt / (mu_0 * dx)
# ダイポールの位置と長さ
center_x, center_y = Nx // 2, Ny // 2
dipole_length = 5 # 格子点数（λ/4 ≈ 5格子点 at λ/20）
# 簡易吸収境界
coeff_abc = (c0 * dt - dx) / (c0 * dt + dx)
Ez_prev = {
'x0': np.zeros(Ny), 'xN': np.zeros(Ny),
'y0': np.zeros(Nx), 'yN': np.zeros(Nx)
}
for n in range(Nt):
Hx[:, :-1] -= Db * (Ez[:, 1:] - Ez[:, :-1])
Hy[:-1, :] += Db * (Ez[1:, :] - Ez[:-1, :])
Ez[1:, 1:] += Cb * (
(Hy[1:, 1:] - Hy[:-1, 1:]) -
(Hx[1:, 1:] - Hx[1:, :-1])
)
# ダイポール励振（正弦波、y方向に沿って配置）
t_now = n * dt
envelope = np.sin(2 * np.pi * freq * t_now)
# ランプアップ（急な開始による過渡応答を抑制）
if t_now < 3 / freq:
envelope *= (1 - np.cos(np.pi * t_now * freq / 3)) / 2
for dy_offset in range(-dipole_length//2, dipole_length//2 + 1):
Ez[center_x, center_y + dy_offset] += envelope * 0.5
# 吸収境界
Ez[0, :] = Ez_prev['x0'] + coeff_abc * (Ez[1, :] - Ez[0, :])
Ez[-1, :] = Ez_prev['xN'] + coeff_abc * (Ez[-2, :] - Ez[-1, :])
Ez[:, 0] = Ez_prev['y0'] + coeff_abc * (Ez[:, 1] - Ez[:, 0])
Ez[:, -1] = Ez_prev['yN'] + coeff_abc * (Ez[:, -1] - Ez[:, -2])
Ez_prev['x0'] = Ez[1, :].copy()
Ez_prev['xN'] = Ez[-2, :].copy()
Ez_prev['y0'] = Ez[:, 1].copy()
Ez_prev['yN'] = Ez[:, -2].copy()
# --- 放射パターンの抽出 ---
# 遠方界の近似: 十分遠い円上で電界を読み取る
r_obs = 80 # 観測半径 [格子点]
n_angles = 360
angles = np.linspace(0, 2*np.pi, n_angles, endpoint=False)
pattern = np.zeros(n_angles)
for i, ang in enumerate(angles):
ix = int(center_x + r_obs * np.cos(ang))
iy = int(center_y + r_obs * np.sin(ang))
if 0 <= ix < Nx and 0 <= iy < Ny:
pattern[i] = np.abs(Ez[ix, iy])
# 正規化
pattern /= np.max(pattern) + 1e-10
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# (a) 電界分布
ax = axes[0]
vmax = np.max(np.abs(Ez)) * 0.3
im = ax.imshow(Ez.T, origin='lower', cmap='RdBu_r',
vmin=-vmax, vmax=vmax,
extent=[0, Nx*dx*100, 0, Ny*dy*100])
ax.set_title('Dipole Radiation: $E_z$ Field', fontsize=14)
ax.set_xlabel('x [cm]', fontsize=12)
ax.set_ylabel('y [cm]', fontsize=12)
# ダイポール位置をマーク
ax.plot([center_x*dx*100], [center_y*dy*100], 'k+', markersize=15,
markeredgewidth=2)
plt.colorbar(im, ax=ax, label='$E_z$ [V/m]', shrink=0.8)
# (b) 放射パターン（極座標）
ax = plt.subplot(122, projection='polar')
ax.plot(angles, pattern, 'c-', linewidth=2)
ax.set_title('Radiation Pattern\n(from FDTD)', fontsize=14, pad=20)
ax.set_rlabel_position(45)
plt.tight_layout()
plt.show()
ダイポールのシミュレーション結果から、以下の重要な特徴が確認できます。
- 8の字型の放射パターン: ダイポールの軸（y方向）に直交する方向（x方向）で放射が最大となり、軸方向では放射がほぼゼロ（ヌル）になっています。これはダイポールアンテナの理論パターン $|\sin\theta|$ と一致します。
- 電界分布の空間構造: 電界のスナップショットから、ダイポール近傍の近傍界（急激に変化する電界）と遠方界（規則的な波面）の違いが視覚的に把握できます。
- 数値的な検証: FDTDで得られた放射パターンが解析解と定性的に一致することで、実装の正しさが検証されます。格子の粗さ（$\lambda/20$）やダイポールの有限長による誤差はありますが、基本的な物理を正しく捉えています。
FDTDの実用上の注意点
数値分散
FDTD法では、電磁波の伝搬速度が方向によってわずかに異なる数値分散が生じます。格子上の離散的な伝搬は、連続空間の等方的な伝搬とは完全には一致しません。
数値分散の大きさは $\Delta x / \lambda$ に依存し、$\Delta x = \lambda/20$ のとき位相速度の誤差は約0.1%です。高い精度が要求される場合は $\Delta x = \lambda/30$ 以上にする必要があります。
計算コスト
FDTD法の計算コストは、3次元の場合 $O(N^3 \times N_t)$ です（$N$ は各方向の格子数）。格子間隔を半分にすると、格子数は8倍、CFL条件により時間ステップ数も2倍になるため、計算量は16倍に増加します。
大規模な問題では、GPUによる並列計算やサブグリッディング（局所的に格子を細かくする手法）が必要になります。
階段近似
FDTD法の格子は直交格子（正方形/直方体）であるため、曲面や斜面は階段状（staircasing）に近似されます。この近似は、曲面上の電界分布に誤差を生じさせます。
階段近似の影響を低減する手法としては、格子を十分細かくする、等角格子（conformal meshing）を用いる、サブセル法を適用するなどがあります。
まとめ
本記事では、FDTD法による電磁界シミュレーションの基礎と実装について解説しました。
- FDTD法の基本原理: マクスウェル方程式の時間微分と空間微分を中心差分で近似し、時間ステップごとに電磁界を更新します
- Yeeセル: $\bm{E}$ と $\bm{H}$ を空間的・時間的に半格子ずらして配置する美しい構造が、中心差分の自然な適用を可能にします
- CFL安定条件: 時間ステップは $\Delta t \leq \Delta/(c\sqrt{n})$（$n$ は次元数）を満たす必要があり、これを超えるとシミュレーションが発散します
- PML（吸収境界条件）: 計算領域の境界に損失層を配置して外向き波を反射なく吸収します。座標伸長に基づく理論的基盤を持ちます
- Pythonによる2D FDTD実装: 点波源とダイポールのシミュレーションにより、電磁波の伝搬と放射パターンの方向依存性を確認しました
FDTD法は、その概念のシンプルさと幅広い適用性から、初学者がまず学ぶべき電磁界解析手法です。本記事で示した2次元実装を3次元に拡張し、現実のアンテナ構造に適用することが次のステップとなります。
次のステップとして、以下の記事も参考にしてください。