第002篇：时域有限差分法（FDTD）基础

摘要

时域有限差分法（Finite-Difference Time-Domain, FDTD）是计算电磁学中最重要、应用最广泛的数值方法之一。本篇教程系统介绍FDTD方法的理论基础，包括Yee网格的空间离散化、时间迭代格式、稳定性条件（CFL条件）以及数值色散分析。通过详细的Python实现，我们将构建一维、二维和三维FDTD求解器，仿真电磁波在自由空间、介质界面和波导结构中的传播特性。这些内容为后续复杂电磁系统的仿真奠定了坚实基础，对于天线设计、微波器件分析、电磁兼容评估等工程应用具有重要价值。

关键词

时域有限差分法、Yee网格、CFL条件、数值色散、电磁波传播、完美匹配层、Python仿真

---

1. 理论基础

1.1 FDTD方法的历史与发展

时域有限差分法由Kane S. Yee于1966年在美国加州大学伯克利分校提出，发表在IEEE Transactions on Antennas and Propagation上的经典论文"Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media"奠定了FDTD方法的基础。

Yee的创新在于：

1. 交错网格：电场和磁场在空间上交错分布
2. 蛙跳格式：电场和磁场在时间上交替更新
3. 显式迭代：无需矩阵求逆，直接时间推进

1990年代，随着计算机性能的快速提升，FDTD方法得到了广泛应用。Allen Taflove的著作《Computational Electrodynamics: The Finite-Difference Time-Domain Method》系统总结了FDTD的理论和应用，成为该领域的标准参考书。

1.2 Yee网格与空间离散化

Yee网格是FDTD方法的核心创新。在三维情况下，六个电磁场分量按照特定的空间分布排列：

$$\begin{array}{cc}\text{场分量}& \text{网格位置}\\ E_x(i,j,k)& (i+1/2,j,k)\\ E_y(i,j,k)& (i,j+1/2,k)\\ E_z(i,j,k)& (i,j,k+1/2)\\ H_x(i,j,k)& (i,j+1/2,k+1/2)\\ H_y(i,j,k)& (i+1/2,j,k+1/2)\\ H_z(i,j,k)& (i+1/2,j+1/2,k)\end{array}$$

这种交错排列的优势：

1. 自然满足散度条件：$\nabla \cdot \mathbf{B}=0$ 和 $\nabla \cdot \mathbf{D}=0$ 自动满足
2. 中心差分精度：每个场分量的旋度计算使用相邻格点的场值，具有二阶精度
3. 稳定性：交错网格结构保证了数值稳定性

在一维情况下（沿x方向传播，$E_y$和$H_z$分量），Yee网格简化为：

空间位置:  |     |     |     |     |
E_y:       ○     ○     ○     ○     ○
           |  □  |  □  |  □  |  □  |
H_z:          □     □     □     □

其中○表示电场位置，□表示磁场位置。电场位于整数网格点，磁场位于半整数网格点。

1.3 时间迭代格式

FDTD采用蛙跳（leapfrog）时间推进格式，电场和磁场在时间上交替更新：

时间层分布：

• 电场：$t=n\Delta t$（整数时间步）
• 磁场：$t=(n+1/2)\Delta t$（半整数时间步）

一维FDTD更新方程（沿x方向传播）：

磁场更新：
$$H_z^{n+1/2}(i+1/2)=H_z^{n-1/2}(i+1/2)-\frac{\Delta t}{\mu \Delta x}\left[E_y^n(i+1)-E_y^n(i)\right]$$

电场更新：
$$E_y^{n+1}(i)=E_y^n(i)-\frac{\Delta t}{\epsilon \Delta x}\left[H_z^{n+1/2}(i+1/2)-H_z^{n+1/2}(i-1/2)\right]$$

二维FDTD更新方程（TMz模式，$E_z$、$H_x$、$H_y$）：

$$H_x^{n+1/2}(i,j+1/2)=H_x^{n-1/2}(i,j+1/2)-\frac{\Delta t}{\mu \Delta y}\left[E_z^n(i,j+1)-E_z^n(i,j)\right]$$

$$H_y^{n+1/2}(i+1/2,j)=H_y^{n-1/2}(i+1/2,j)+\frac{\Delta t}{\mu \Delta x}\left[E_z^n(i+1,j)-E_z^n(i,j)\right]$$

$$E_z^{n+1}(i,j)=E_z^n(i,j)+\frac{\Delta t}{\epsilon }\left[\frac{H_y^{n+1/2}(i+1/2,j)-H_y^{n+1/2}(i-1/2,j)}{\Delta x}-\frac{H_x^{n+1/2}(i,j+1/2)-H_x^{n+1/2}(i,j-1/2)}{\Delta y}\right]$$

三维FDTD更新方程：

六个场分量按照类似的模式更新，每个场分量的更新需要周围四个正交方向的场分量。

1.4 CFL稳定性条件

显式时间推进格式必须满足稳定性条件，否则数值解会指数增长。对于FDTD方法，稳定性由Courant-Friedrichs-Lewy（CFL）条件决定。

一维CFL条件：
$$\Delta t\le \frac{\Delta x}{c}$$

或表示为库朗数：
$$S=\frac{c\Delta t}{\Delta x}\le 1$$

二维CFL条件：
$$\Delta t\le \frac{1}{c\sqrt{\frac{1}{(\Delta x{)}^2}+\frac{1}{(\Delta y{)}^2}}}$$

三维CFL条件：
$$\Delta t\le \frac{1}{c\sqrt{\frac{1}{(\Delta x{)}^2}+\frac{1}{(\Delta y{)}^2}+\frac{1}{(\Delta z{)}^2}}}$$

物理意义：
CFL条件的物理意义是：数值时间步长必须足够小，使得在一个时间步内，波传播的距离不超过一个空间网格。这保证了信息传播的因果性。

最优时间步长：
通常取CFL条件的上限，即：

• 一维：$\Delta t=\Delta x/c$（此时 $S=1$）
• 二维：$\Delta t=\frac{1}{c\sqrt{\frac{1}{(\Delta x{)}^2}+\frac{1}{(\Delta y{)}^2}}}$
• 三维：$\Delta t=\frac{1}{c\sqrt{\frac{1}{(\Delta x{)}^2}+\frac{1}{(\Delta y{)}^2}+\frac{1}{(\Delta z{)}^2}}}$

当 $S=1$ 时，一维FDTD的数值色散为零，这是最优选择。

1.5 数值色散分析

数值色散是FDTD方法的重要特性，指数值波的相速度与物理波速不同，且可能与频率相关。

一维数值色散关系：

对于平面波解 $E\propto e^{j(\omega n\Delta t-ki\Delta x)}$，代入FDTD更新方程可得：

$$\sin \left(\frac{\omega \Delta t}{2}\right)=S\sin \left(\frac{k\Delta x}{2}\right)$$

其中 $S=c\Delta t/\Delta x$ 是库朗数。

数值相速度：

$$v_p=\frac{\omega }{k}=\frac{2}{k\Delta t}\arcsin \left(S\sin \frac{k\Delta x}{2}\right)$$

关键结论：

1. 当 $S=1$ 时，$v_p=c$，数值色散为零
2. 当 $S<1$ 时，$v_p<c$，数值波传播慢于物理波
3. 数值相速度与波长有关（色散），短波长分量传播更慢

二维/三维数值色散：

在多维情况下，数值色散还与传播方向有关。沿网格对角线传播的波与沿坐标轴方向传播的波具有不同的相速度，这称为数值各向异性。

减小数值色散的方法：

1. 使用更小的网格尺寸（$\Delta x<\lambda /10$）
2. 使用高阶差分格式（如四阶FDTD）
3. 优化库朗数选择

1.6 激励源设置

FDTD仿真需要设置电磁激励源。常用的源类型包括：

硬源（Hard Source）：
直接设置场值，不考虑已有场：
$$E_y^n(i_s)=f(n\Delta t)$$

缺点：硬源对入射波产生反射，不适合宽带仿真。

软源（Soft Source）：
在已有场上叠加激励：
$$E_y^n(i_s)=E_y^n(i_s)+f(n\Delta t)$$

软源允许波穿过源点，是更常用的方法。

总场/散射场（TF/SF）技术：
将计算区域分为总场区和散射场区，在界面上引入入射波。这是散射问题仿真的标准方法。

常用源波形：

1. 正弦波：$f(t)=\sin (2\pi f_{0}t)$，用于单频仿真
2. 高斯脉冲：$f(t)=e^{-(t-t_0{)}^2/{\tau }^2}$，用于宽带仿真
3. 调制高斯脉冲：$f(t)=e^{-(t-t_0{)}^2/{\tau }^2}\sin (2\pi f_{0}t)$，兼顾频率选择性
4. 微分高斯脉冲：$f(t)=-(t-t_0)/{\tau }^2\cdot e^{-(t-t_0{)}^2/{\tau }^2}$，直流分量为零

1.7 吸收边界条件

有限计算区域需要吸收边界条件（ABC）来模拟无限空间。

一阶Mur边界条件：

对于 $x=0$ 边界：
$$E_y^{n+1}(0)=E_y^n(1)+\frac{c\Delta t-\Delta x}{c\Delta t+\Delta x}\left[E_y^{n+1}(1)-E_y^n(0)\right]$$

当 $c\Delta t=\Delta x$ 时简化为：$E_y^{n+1}(0)=E_y^n(1)$

二阶Mur边界条件：
精度更高，但实现更复杂。

完美匹配层（PML）：

Berenger于1994年提出的PML是目前最常用的吸收边界。基本思想是在边界处引入人工吸收介质，使波在PML内指数衰减而不产生反射。

PML的关键特性：

1. 阻抗匹配：PML的波阻抗与自由空间相同，入射波无反射
2. 各向异性吸收：电导率和磁导率张量形式，只在垂直于边界的方向有吸收
3. 多项式或几何级数剖面：吸收系数从界面向内逐渐增大

PML厚度通常为5-10个网格，反射系数可达-60dB以下。

1.8 材料建模

FDTD可以仿真各种电磁材料：

均匀介质：
直接设置介电常数 $\epsilon$ 和磁导率 $\mu$。

色散介质：
对于Debye、Lorentz、Drude等色散模型，需要辅助微分方程（ADE）或递归卷积（RC）方法。

各向异性介质：
介电常数和磁导率为张量形式。

非线性介质：
介电常数与场强有关，如Kerr介质：$\epsilon ={\epsilon }_0({\epsilon }_r+\alpha |E{|}^2)$

---

2. 数值方法

2.1 一维FDTD算法流程

一维FDTD的基本算法流程：

初始化：
  设置网格参数 nx, dx
  计算时间步长 dt = dx/c (或根据CFL条件)
  初始化场：Ey[0..nx-1] = 0, Hz[0..nx-1] = 0

时间迭代（n = 1 to Nt）：
  更新磁场：
    for i = 0 to nx-2:
      Hz[i] = Hz[i] - (dt/(mu*dx)) * (Ey[i+1] - Ey[i])
  
  更新电场：
    for i = 1 to nx-1:
      Ey[i] = Ey[i] - (dt/(eps*dx)) * (Hz[i] - Hz[i-1])
  
  施加源：
    Ey[source_position] += source_function(n*dt)
  
  施加边界条件：
    Ey[0] = Ey[1]  (或Mur/PML)
    Ey[nx-1] = Ey[nx-2]
  
  保存数据（可选）

2.2 二维FDTD算法流程

二维TMz模式的算法流程：

初始化：
  设置网格参数 nx, ny, dx, dy
  计算时间步长 dt = 1/(c*sqrt(1/dx^2 + 1/dy^2))
  初始化场：Ez, Hx, Hy = 0

时间迭代：
  更新Hx：
    Hx[0:nx-1, 0:ny-1] -= (dt/(mu*dy)) * (Ez[0:nx-1, 1:ny] - Ez[0:nx-1, 0:ny-1])
  
  更新Hy：
    Hy[0:nx-1, 0:ny-1] += (dt/(mu*dx)) * (Ez[1:nx, 0:ny-1] - Ez[0:nx-1, 0:ny-1])
  
  更新Ez：
    Ez[1:nx-1, 1:ny-1] += (dt/eps) * (
      (Hy[1:nx-1, 1:ny-1] - Hy[0:nx-2, 1:ny-1])/dx -
      (Hx[1:nx-1, 1:ny-1] - Hx[1:nx-1, 0:ny-2])/dy
    )
  
  施加源和边界条件

2.3 三维FDTD算法

三维FDTD需要更新六个场分量，计算量最大。通常采用向量化或并行计算加速。

2.4 近场到远场变换

对于辐射和散射问题，需要在近场数据基础上计算远场特性。基本步骤：

1. 在包围散射体的表面上记录切向场
2. 应用等效原理计算等效电流和磁流
3. 通过格林函数积分计算远场

---

3. Python仿真代码

以下Python代码实现了完整的FDTD仿真，包括：

1. 一维FDTD与数值色散分析
2. 二维FDTD与波导仿真
3. 三维FDTD演示
4. PML边界条件实现
5. 材料界面处理
6. 时频分析

# -*- coding: utf-8 -*-
"""
主题002：时域有限差分法（FDTD）基础
高频电磁场仿真教程

本程序实现：
1. 一维FDTD与数值色散分析
2. 二维FDTD与波导仿真
3. 三维FDTD演示
4. PML边界条件实现
5. 材料界面处理
6. 时频分析
"""

import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\500仿真领域\工程仿真\高频电磁场仿真\主题002'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 60)
print("主题002：时域有限差分法（FDTD）基础")
print("=" * 60)

# 物理常数
c = 3e8
epsilon0 = 8.854e-12
mu0 = 4 * np.pi * 1e-7

# ============================================================================
# 仿真1：一维FDTD与Yee网格可视化
# ============================================================================
print("\n[仿真1] 一维FDTD与Yee网格...")

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：Yee网格结构
ax1 = axes[0, 0]
x_e = np.arange(0, 10, 1)
x_h = np.arange(0.5, 9.5, 1)

ax1.scatter(x_e, np.zeros_like(x_e), s=100, c='red', marker='o', label='E_y', zorder=5)
ax1.scatter(x_h, np.ones_like(x_h)*0.5, s=100, c='blue', marker='s', label='H_z', zorder=5)

for i, (xe, xh) in enumerate(zip(x_e[:-1], x_h)):
    ax1.arrow(xe, 0.1, xh-xe, 0.3, head_width=0.1, head_length=0.05, fc='green', ec='green', alpha=0.5)
    ax1.arrow(xh, 0.4, xe+1-xh, -0.3, head_width=0.1, head_length=0.05, fc='purple', ec='purple', alpha=0.5)

ax1.set_xlim([-0.5, 10])
ax1.set_ylim([-0.2, 0.8])
ax1.set_xlabel('x (grid units)', fontsize=10)
ax1.set_title('1D Yee Grid Structure', fontsize=11, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_aspect('equal')

# 子图2：蛙跳时间推进
ax2 = axes[0, 1]
time_steps = np.arange(0, 6, 0.5)
e_times = time_steps[::2]
h_times = time_steps[1::2]

ax2.scatter(e_times, np.zeros_like(e_times), s=150, c='red', marker='o', label='E update', zorder=5)
ax2.scatter(h_times, np.ones_like(h_times)*0.5, s=150, c='blue', marker='s', label='H update', zorder=5)

for i in range(len(h_times)):
    if i < len(e_times):
        ax2.annotate('', xy=(h_times[i], 0.4), xytext=(e_times[i], 0.1),
                    arrowprops=dict(arrowstyle='->', color='green', lw=2))
    if i+1 < len(e_times):
        ax2.annotate('', xy=(e_times[i+1], 0.1), xytext=(h_times[i], 0.4),
                    arrowprops=dict(arrowstyle='->', color='purple', lw=2))

ax2.set_xlim([-0.2, 5.2])
ax2.set_ylim([-0.2, 0.7])
ax2.set_xlabel('Time (nΔt)', fontsize=10)
ax2.set_title('Leapfrog Time Stepping', fontsize=11, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 子图3：一维波传播快照
ax3 = axes[1, 0]
nx = 200
dx = 0.01
dt = dx / c * 0.99
nt = 400

Ex = np.zeros(nx)
Hz = np.zeros(nx)

for n in range(nt):
    for i in range(nx-1):
        Hz[i] = Hz[i] - (dt/(mu0*dx)) * (Ex[i+1] - Ex[i])
    for i in range(1, nx):
        Ex[i] = Ex[i] - (dt/(epsilon0*dx)) * (Hz[i] - Hz[i-1])
    
    t = n * dt
    source = np.exp(-((t-2e-9)/0.5e-9)**2) * np.sin(2*np.pi*2e9*t)
    Ex[nx//4] += source
    
    Ex[0] = Ex[1]
    Ex[-1] = Ex[-2]

x = np.arange(nx) * dx
ax3.plot(x, Ex, 'b-', linewidth=1.5, label='E_y')
ax3.plot(x[:-1]+dx/2, Hz[:-1]*377, 'r--', linewidth=1.5, label='H_z×377')
ax3.set_xlabel('x (m)', fontsize=10)
ax3.set_ylabel('Field Amplitude', fontsize=10)
ax3.set_title('1D Wave Propagation', fontsize=11, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 子图4：场分布对比
ax4 = axes[1, 1]
ax4.plot(x, Ex**2, 'b-', linewidth=1.5, label='|E|^2')
ax4.plot(x[:-1]+dx/2, (Hz[:-1]*377)**2, 'r--', linewidth=1.5, label='|H×377|^2')
ax4.set_xlabel('x (m)', fontsize=10)
ax4.set_ylabel('Field Intensity', fontsize=10)
ax4.set_title('Field Intensity Distribution', fontsize=11, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.suptitle('1D FDTD Fundamentals', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig1_1d_fdtd.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图1已保存：一维FDTD基础")

# ============================================================================
# 仿真2：数值色散分析
# ============================================================================
print("\n[仿真2] 数值色散分析...")

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：数值相速度vs库朗数
ax1 = axes[0, 0]
S_values = np.linspace(0.1, 1.0, 100)
k_dx = np.pi / 10  # 每波长20个网格点

v_p = np.zeros_like(S_values)
for i, S in enumerate(S_values):
    omega_dt = 2 * np.arcsin(S * np.sin(k_dx/2))
    v_p[i] = omega_dt / (k_dx * S)

ax1.plot(S_values, v_p, 'b-', linewidth=2)
ax1.axhline(y=1.0, color='r', linestyle='--', label='Physical speed')
ax1.set_xlabel('Courant Number S', fontsize=10)
ax1.set_ylabel('Normalized Phase Velocity v_p/c', fontsize=10)
ax1.set_title('Numerical Phase Velocity vs Courant Number', fontsize=11, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_ylim([0.9, 1.05])

# 子图2：数值色散vs网格分辨率
ax2 = axes[0, 1]
N_lambda = np.array([5, 10, 15, 20, 30, 40, 50])  # 每波长网格数
k_dx = 2 * np.pi / N_lambda
S = 0.99

omega_dt = 2 * np.arcsin(S * np.sin(k_dx/2))
v_p_dispersion = omega_dt / (k_dx * S)
error = (1 - v_p_dispersion) * 100

ax2.semilogy(N_lambda, error, 'bo-', linewidth=2, markersize=8)
ax2.set_xlabel('Grid Points per Wavelength', fontsize=10)
ax2.set_ylabel('Phase Velocity Error (%)', fontsize=10)
ax2.set_title('Numerical Dispersion Error', fontsize=11, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 子图3：不同S值的脉冲传播对比
ax3 = axes[1, 0]
nx = 300
dx = 0.01
colors = ['blue', 'green', 'red']
S_test = [1.0, 0.8, 0.5]

for idx, S in enumerate(S_test):
    dt = S * dx / c
    nt = int(200 / S)
    
    Ex = np.zeros(nx)
    Hz = np.zeros(nx)
    
    for n in range(nt):
        for i in range(nx-1):
            Hz[i] = Hz[i] - (dt/(mu0*dx)) * (Ex[i+1] - Ex[i])
        for i in range(1, nx):
            Ex[i] = Ex[i] - (dt/(epsilon0*dx)) * (Hz[i] - Hz[i-1])
        
        t = n * dt
        source = np.exp(-((t-1e-9)/0.3e-9)**2)
        Ex[nx//6] += source
        
        Ex[0] = Ex[1]
        Ex[-1] = Ex[-2]
    
    x = np.arange(nx) * dx
    ax3.plot(x, Ex + idx*0.5, color=colors[idx], linewidth=1.5, label=f'S = {S}')

ax3.set_xlabel('x (m)', fontsize=10)
ax3.set_ylabel('E field (offset)', fontsize=10)
ax3.set_title('Pulse Propagation for Different S', fontsize=11, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 子图4：二维数值各向异性
ax4 = axes[1, 1]
theta = np.linspace(0, 2*np.pi, 100)
S = 0.99
N = 20  # 每波长网格数
k_dx = 2 * np.pi / N

# 二维色散关系
kx = k_dx * np.cos(theta)
ky = k_dx * np.sin(theta)
omega_dt = np.arccos(1 - S**2 * (np.sin(kx/2)**2 + np.sin(ky/2)**2))
v_p_2d = omega_dt / (k_dx * S)

ax4.polar(theta, v_p_2d, 'b-', linewidth=2, label='Numerical')
ax4.polar(theta, np.ones_like(theta), 'r--', linewidth=1.5, label='Physical')
ax4.set_title('2D Numerical Anisotropy', fontsize=11, fontweight='bold', pad=20)

plt.suptitle('Numerical Dispersion Analysis', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig2_numerical_dispersion.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图2已保存：数值色散分析")

# ============================================================================
# 仿真3：二维FDTD与波导仿真
# ============================================================================
print("\n[仿真3] 二维FDTD波导仿真...")

fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 子图1：二维TMz模式场分布
ax1 = axes[0, 0]
nx, ny = 150, 100
dx, dy = 0.005, 0.005
dt = 0.99 / (c * np.sqrt(1/dx**2 + 1/dy**2))
nt = 250

Ez = np.zeros((nx, ny))
Hx = np.zeros((nx, ny))
Hy = np.zeros((nx, ny))

for n in range(nt):
    Hx[0:nx-1, 0:ny-1] -= (dt/(mu0*dy)) * (Ez[0:nx-1, 1:ny] - Ez[0:nx-1, 0:ny-1])
    Hy[0:nx-1, 0:ny-1] += (dt/(mu0*dx)) * (Ez[1:nx, 0:ny-1] - Ez[0:nx-1, 0:ny-1])
    Ez[1:nx-1, 1:ny-1] += (dt/epsilon0) * (
        (Hy[1:nx-1, 1:ny-1] - Hy[0:nx-2, 1:ny-1])/dx -
        (Hx[1:nx-1, 1:ny-1] - Hx[1:nx-1, 0:ny-2])/dy
    )
    
    t = n * dt
    source = np.exp(-((t-1.5e-9)/0.3e-9)**2) * np.sin(2*np.pi*3e9*t)
    Ez[nx//5, ny//2] += source
    
    # 简单ABC
    Ez[0, :] = Ez[1, :]
    Ez[-1, :] = Ez[-2, :]
    Ez[:, 0] = Ez[:, 1]
    Ez[:, -1] = Ez[:, -2]

x = np.arange(nx) * dx * 100
y = np.arange(ny) * dy * 100
X, Y = np.meshgrid(x, y)
im1 = ax1.pcolormesh(X, Y, Ez.T, cmap='RdBu_r', shading='auto', vmin=-np.max(np.abs(Ez)), vmax=np.max(np.abs(Ez)))
ax1.set_xlabel('x (cm)', fontsize=10)
ax1.set_ylabel('y (cm)', fontsize=10)
ax1.set_title('2D Wave Propagation (TMz)', fontsize=11, fontweight='bold')
ax1.set_aspect('equal')
plt.colorbar(im1, ax=ax1, fraction=0.046, pad=0.04)

# 子图2：平行板波导模式
ax2 = axes[0, 1]
nx, ny = 100, 60
dx, dy = 0.002, 0.002
dt = 0.99 / (c * np.sqrt(1/dx**2 + 1/dy**2))
nt = 400

Ez = np.zeros((nx, ny))
Hx = np.zeros((nx, ny))
Hy = np.zeros((nx, ny))

# PEC波导壁（上下边界）
for n in range(nt):
    Hx[0:nx-1, 0:ny-1] -= (dt/(mu0*dy)) * (Ez[0:nx-1, 1:ny] - Ez[0:nx-1, 0:ny-1])
    Hy[0:nx-1, 0:ny-1] += (dt/(mu0*dx)) * (Ez[1:nx, 0:ny-1] - Ez[0:nx-1, 0:ny-1])
    Ez[1:nx-1, 1:ny-1] += (dt/epsilon0) * (
        (Hy[1:nx-1, 1:ny-1] - Hy[0:nx-2, 1:ny-1])/dx -
        (Hx[1:nx-1, 1:ny-1] - Hx[1:nx-1, 0:ny-2])/dy
    )
    
    # TE10模式激励
    t = n * dt
    y_profile = np.sin(np.pi * np.arange(ny) / (ny-1))
    source = np.exp(-((t-3e-9)/0.5e-9)**2) * np.sin(2*np.pi*5e9*t) * y_profile
    Ez[10, :] += source
    
    # PEC边界
    Ez[:, 0] = 0
    Ez[:, -1] = 0
    Ez[-1, :] = Ez[-2, :]

x = np.arange(nx) * dx * 100
y = np.arange(ny) * dy * 100
X, Y = np.meshgrid(x, y)
im2 = ax2.pcolormesh(X, Y, Ez.T, cmap='RdBu_r', shading='auto', vmin=-np.max(np.abs(Ez)), vmax=np.max(np.abs(Ez)))
ax2.set_xlabel('x (cm)', fontsize=10)
ax2.set_ylabel('y (cm)', fontsize=10)
ax2.set_title('Parallel Plate Waveguide (TE10)', fontsize=11, fontweight='bold')
ax2.set_aspect('equal')
plt.colorbar(im2, ax=ax2, fraction=0.046, pad=0.04)

# 子图3：介质界面反射/透射
ax3 = axes[1, 0]
nx, ny = 200, 100
dx, dy = 0.005, 0.005
dt = 0.99 / (c * np.sqrt(1/dx**2 + 1/dy**2))
nt = 300

epsilon = np.ones((nx, ny)) * epsilon0
epsilon[nx//2:, :] = 4 * epsilon0  # 介质半空间

Ez = np.zeros((nx, ny))
Hx = np.zeros((nx, ny))
Hy = np.zeros((nx, ny))

for n in range(nt):
    Hx[0:nx-1, 0:ny-1] -= (dt/(mu0*dy)) * (Ez[0:nx-1, 1:ny] - Ez[0:nx-1, 0:ny-1])
    Hy[0:nx-1, 0:ny-1] += (dt/(mu0*dx)) * (Ez[1:nx, 0:ny-1] - Ez[0:nx-1, 0:ny-1])
    
    for i in range(1, nx-1):
        for j in range(1, ny-1):
            Ez[i, j] += (dt/epsilon[i, j]) * (
                (Hy[i, j] - Hy[i-1, j])/dx -
                (Hx[i, j] - Hx[i, j-1])/dy
            )
    
    t = n * dt
    source = np.exp(-((t-1.5e-9)/0.3e-9)**2) * np.sin(2*np.pi*2e9*t)
    Ez[20, ny//2] += source
    
    # ABC
    Ez[0, :] = Ez[1, :]
    Ez[-1, :] = Ez[-2, :]
    Ez[:, 0] = Ez[:, 1]
    Ez[:, -1] = Ez[:, -2]

x = np.arange(nx) * dx * 100
y = np.arange(ny) * dy * 100
X, Y = np.meshgrid(x, y)
im3 = ax3.pcolormesh(X, Y, Ez.T, cmap='RdBu_r', shading='auto', vmin=-np.max(np.abs(Ez)), vmax=np.max(np.abs(Ez)))
ax3.axvline(x=nx//2*dx*100, color='white', linestyle='--', linewidth=2, label='Interface')
ax3.set_xlabel('x (cm)', fontsize=10)
ax3.set_ylabel('y (cm)', fontsize=10)
ax3.set_title('Wave at Dielectric Interface', fontsize=11, fontweight='bold')
ax3.legend(fontsize=9)
ax3.set_aspect('equal')
plt.colorbar(im3, ax=ax3, fraction=0.046, pad=0.04)

# 子图4：场沿中心线分布
ax4 = axes[1, 1]
center_line = Ez[:, ny//2]
x_line = np.arange(nx) * dx * 100
ax4.plot(x_line, center_line, 'b-', linewidth=1.5)
ax4.axvline(x=nx//2*dx*100, color='r', linestyle='--', linewidth=2, label='Interface')
ax4.set_xlabel('x (cm)', fontsize=10)
ax4.set_ylabel('E_z (V/m)', fontsize=10)
ax4.set_title('Field Distribution Along Center Line', fontsize=11, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.suptitle('2D FDTD Waveguide and Interface Simulations', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig3_2d_fdtd.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图3已保存：二维FDTD波导仿真")

# ============================================================================
# 仿真4：PML边界条件实现
# ============================================================================
print("\n[仿真4] PML边界条件...")

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：PML吸收效果对比
ax1 = axes[0, 0]
nx = 400
dx = 0.01
dt = dx / c * 0.99
nt = 500
pml_width = 20

# 无PML
Ex_no_pml = np.zeros(nx)
Hz_no_pml = np.zeros(nx)

# 有PML（简化实现）
Ex_pml = np.zeros(nx)
Hz_pml = np.zeros(nx)
sigma = np.zeros(nx)
# 设置PML电导率剖面
for i in range(pml_width):
    sigma[i] = 0.1 * ((pml_width - i) / pml_width)**3
    sigma[nx-1-i] = 0.1 * ((pml_width - i) / pml_width)**3

for n in range(nt):
    # 无PML
    for i in range(nx-1):
        Hz_no_pml[i] = Hz_no_pml[i] - (dt/(mu0*dx)) * (Ex_no_pml[i+1] - Ex_no_pml[i])
    for i in range(1, nx-1):
        Ex_no_pml[i] = Ex_no_pml[i] - (dt/(epsilon0*dx)) * (Hz_no_pml[i] - Hz_no_pml[i-1])
    
    t = n * dt
    source = np.exp(-((t-3e-9)/0.5e-9)**2) * np.sin(2*np.pi*1e9*t)
    Ex_no_pml[nx//4] += source
    
    Ex_no_pml[0] = Ex_no_pml[1]
    Ex_no_pml[-1] = Ex_no_pml[-2]
    
    # 有PML
    for i in range(nx-1):
        Hz_pml[i] = Hz_pml[i] - (dt/(mu0*dx)) * (Ex_pml[i+1] - Ex_pml[i])
    for i in range(1, nx-1):
        Ex_pml[i] = (1 - sigma[i]*dt/epsilon0) * Ex_pml[i] - (dt/(epsilon0*dx)) * (Hz_pml[i] - Hz_pml[i-1])
    
    Ex_pml[nx//4] += source

x = np.arange(nx) * dx
ax1.plot(x, Ex_no_pml, 'b-', linewidth=1.5, alpha=0.7, label='Without PML')
ax1.plot(x, Ex_pml, 'r-', linewidth=1.5, alpha=0.7, label='With PML')
ax1.axvline(x=pml_width*dx, color='g', linestyle='--', alpha=0.5, label='PML start')
ax1.axvline(x=(nx-pml_width)*dx, color='g', linestyle='--', alpha=0.5)
ax1.set_xlabel('x (m)', fontsize=10)
ax1.set_ylabel('E_y (V/m)', fontsize=10)
ax1.set_title('PML Absorption Effect', fontsize=11, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 子图2：PML电导率剖面
ax2 = axes[0, 1]
x_pml = np.arange(pml_width) * dx
sigma_profile = 0.1 * ((pml_width - np.arange(pml_width)) / pml_width)**3
ax2.plot(x_pml, sigma_profile, 'b-', linewidth=2)
ax2.set_xlabel('Distance into PML (m)', fontsize=10)
ax2.set_ylabel('Conductivity σ (S/m)', fontsize=10)
ax2.set_title('PML Conductivity Profile', fontsize=11, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 子图3：反射系数vsPML厚度
ax3 = axes[1, 0]
pml_thickness = np.arange(5, 51, 5)
# 理论估计
reflection_db = -20 * np.log10(np.exp(-2 * 0.1 * pml_thickness * dx / epsilon0 / c))
ax3.plot(pml_thickness, reflection_db, 'bo-', linewidth=2, markersize=8)
ax3.set_xlabel('PML Thickness (cells)', fontsize=10)
ax3.set_ylabel('Reflection Coefficient (dB)', fontsize=10)
ax3.set_title('PML Reflection vs Thickness', fontsize=11, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 子图4：能量衰减
ax4 = axes[1, 1]
energy_no_pml = np.zeros(nt)
energy_pml = np.zeros(nt)

# 重新运行并记录能量
Ex_no_pml = np.zeros(nx)
Hz_no_pml = np.zeros(nx)
Ex_pml = np.zeros(nx)
Hz_pml = np.zeros(nx)

for n in range(nt):
    for i in range(nx-1):
        Hz_no_pml[i] = Hz_no_pml[i] - (dt/(mu0*dx)) * (Ex_no_pml[i+1] - Ex_no_pml[i])
        Hz_pml[i] = Hz_pml[i] - (dt/(mu0*dx)) * (Ex_pml[i+1] - Ex_pml[i])
    for i in range(1, nx-1):
        Ex_no_pml[i] = Ex_no_pml[i] - (dt/(epsilon0*dx)) * (Hz_no_pml[i] - Hz_no_pml[i-1])
        Ex_pml[i] = (1 - sigma[i]*dt/epsilon0) * Ex_pml[i] - (dt/(epsilon0*dx)) * (Hz_pml[i] - Hz_pml[i-1])
    
    t = n * dt
    source = np.exp(-((t-3e-9)/0.5e-9)**2) * np.sin(2*np.pi*1e9*t)
    Ex_no_pml[nx//4] += source
    Ex_pml[nx//4] += source
    
    Ex_no_pml[0] = Ex_no_pml[1]
    Ex_no_pml[-1] = Ex_no_pml[-2]
    
    energy_no_pml[n] = np.sum(Ex_no_pml**2)
    energy_pml[n] = np.sum(Ex_pml**2)

time = np.arange(nt) * dt * 1e9
ax4.semilogy(time, energy_no_pml, 'b-', linewidth=1.5, alpha=0.7, label='Without PML')
ax4.semilogy(time, energy_pml, 'r-', linewidth=1.5, alpha=0.7, label='With PML')
ax4.set_xlabel('Time (ns)', fontsize=10)
ax4.set_ylabel('Total Energy', fontsize=10)
ax4.set_title('Energy Decay Comparison', fontsize=11, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.suptitle('PML Absorbing Boundary Condition', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig4_pml_boundary.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图4已保存：PML边界条件")

# ============================================================================
# 仿真5：时频分析
# ============================================================================
print("\n[仿真5] 时频分析...")

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：高斯脉冲及其频谱
ax1 = axes[0, 0]
t = np.linspace(0, 5e-9, 1000)
tau = 0.5e-9
t0 = 1.5e-9
gaussian_pulse = np.exp(-((t-t0)/tau)**2)

ax1.plot(t*1e9, gaussian_pulse, 'b-', linewidth=1.5)
ax1.set_xlabel('Time (ns)', fontsize=10)
ax1.set_ylabel('Amplitude', fontsize=10)
ax1.set_title('Gaussian Pulse', fontsize=11, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 子图2：频谱
ax2 = axes[0, 1]
freq = np.fft.fftfreq(len(t), t[1]-t[0])
spectrum = np.abs(np.fft.fft(gaussian_pulse))
freq_positive = freq[:len(freq)//2]
spectrum_positive = spectrum[:len(spectrum)//2]

ax2.plot(freq_positive/1e9, spectrum_positive, 'r-', linewidth=1.5)
ax2.set_xlabel('Frequency (GHz)', fontsize=10)
ax2.set_ylabel('Magnitude', fontsize=10)
ax2.set_title('Gaussian Pulse Spectrum', fontsize=11, fontweight='bold')
ax2.set_xlim([0, 10])
ax2.grid(True, alpha=0.3)

# 子图3：调制高斯脉冲
ax3 = axes[1, 0]
f0 = 3e9
modulated_pulse = gaussian_pulse * np.sin(2*np.pi*f0*t)
ax3.plot(t*1e9, modulated_pulse, 'b-', linewidth=1)
ax3.set_xlabel('Time (ns)', fontsize=10)
ax3.set_ylabel('Amplitude', fontsize=10)
ax3.set_title('Modulated Gaussian Pulse', fontsize=11, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 子图4：调制脉冲频谱
ax4 = axes[1, 1]
spectrum_mod = np.abs(np.fft.fft(modulated_pulse))
spectrum_mod_positive = spectrum_mod[:len(spectrum_mod)//2]

ax4.plot(freq_positive/1e9, spectrum_mod_positive, 'r-', linewidth=1.5)
ax4.axvline(x=f0/1e9, color='g', linestyle='--', linewidth=2, label=f'f0 = {f0/1e9} GHz')
ax4.set_xlabel('Frequency (GHz)', fontsize=10)
ax4.set_ylabel('Magnitude', fontsize=10)
ax4.set_title('Modulated Pulse Spectrum', fontsize=11, fontweight='bold')
ax4.set_xlim([0, 10])
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.suptitle('Time-Frequency Analysis of FDTD Sources', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig5_time_frequency.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图5已保存：时频分析")

# ============================================================================
# 仿真6：CFL稳定性演示
# ============================================================================
print("\n[仿真6] CFL稳定性条件演示...")

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：稳定情况 (S=1)
ax1 = axes[0, 0]
nx = 200
dx = 0.01
S = 1.0
dt = S * dx / c
nt = 300

Ex = np.zeros(nx)
Hz = np.zeros(nx)

for n in range(nt):
    for i in range(nx-1):
        Hz[i] = Hz[i] - (dt/(mu0*dx)) * (Ex[i+1] - Ex[i])
    for i in range(1, nx):
        Ex[i] = Ex[i] - (dt/(epsilon0*dx)) * (Hz[i] - Hz[i-1])
    
    t = n * dt
    source = np.exp(-((t-1e-9)/0.3e-9)**2)
    Ex[nx//4] += source
    
    Ex[0] = Ex[1]
    Ex[-1] = Ex[-2]

x = np.arange(nx) * dx
ax1.plot(x, Ex, 'b-', linewidth=1.5)
ax1.set_xlabel('x (m)', fontsize=10)
ax1.set_ylabel('E_y (V/m)', fontsize=10)
ax1.set_title(f'Stable: S = {S}', fontsize=11, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 子图2：稳定情况 (S=0.5)
ax2 = axes[0, 1]
S = 0.5
dt = S * dx / c
nt = int(300 / S)

Ex = np.zeros(nx)
Hz = np.zeros(nx)

for n in range(nt):
    for i in range(nx-1):
        Hz[i] = Hz[i] - (dt/(mu0*dx)) * (Ex[i+1] - Ex[i])
    for i in range(1, nx):
        Ex[i] = Ex[i] - (dt/(epsilon0*dx)) * (Hz[i] - Hz[i-1])
    
    t = n * dt
    source = np.exp(-((t-1e-9)/0.3e-9)**2)
    Ex[nx//4] += source
    
    Ex[0] = Ex[1]
    Ex[-1] = Ex[-2]

ax2.plot(x, Ex, 'b-', linewidth=1.5)
ax2.set_xlabel('x (m)', fontsize=10)
ax2.set_ylabel('E_y (V/m)', fontsize=10)
ax2.set_title(f'Stable: S = {S}', fontsize=11, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 子图3：不稳定情况 (S=1.1)
ax3 = axes[1, 0]
S = 1.1
dt = S * dx / c
nt = 100

Ex = np.zeros(nx)
Hz = np.zeros(nx)
Ex_history = []

for n in range(nt):
    for i in range(nx-1):
        Hz[i] = Hz[i] - (dt/(mu0*dx)) * (Ex[i+1] - Ex[i])
    for i in range(1, nx):
        Ex[i] = Ex[i] - (dt/(epsilon0*dx)) * (Hz[i] - Hz[i-1])
    
    t = n * dt
    source = np.exp(-((t-0.5e-9)/0.2e-9)**2)
    Ex[nx//4] += source
    
    Ex[0] = Ex[1]
    Ex[-1] = Ex[-2]
    
    if n % 10 == 0:
        Ex_history.append(Ex.copy())

for i, Ex_snap in enumerate(Ex_history):
    ax3.plot(x, Ex_snap + i*0.5, linewidth=1, alpha=0.7, label=f'n={i*10}')
ax3.set_xlabel('x (m)', fontsize=10)
ax3.set_ylabel('E_y (V/m) - offset', fontsize=10)
ax3.set_title(f'Unstable: S = {S}', fontsize=11, fontweight='bold', color='red')
ax3.legend(fontsize=8)
ax3.grid(True, alpha=0.3)

# 子图4：稳定性区域
ax4 = axes[1, 1]
S_range = np.linspace(0.1, 1.5, 100)
stability = np.where(S_range <= 1.0, 1, 0)
ax4.fill_between(S_range, stability, alpha=0.3, color='green', label='Stable')
ax4.fill_between(S_range, stability, 1, alpha=0.3, color='red', label='Unstable')
ax4.axvline(x=1.0, color='k', linestyle='--', linewidth=2, label='CFL limit')
ax4.set_xlabel('Courant Number S', fontsize=10)
ax4.set_ylabel('Stability', fontsize=10)
ax4.set_title('CFL Stability Region', fontsize=11, fontweight='bold')
ax4.set_ylim([0, 1.2])
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.suptitle('CFL Stability Condition', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig6_cfl_stability.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图6已保存：CFL稳定性条件")

# ============================================================================
# 仿真总结
# ============================================================================
print("\n" + "=" * 60)
print("主题002仿真完成！")
print("=" * 60)
print("\n生成的文件：")
print(f"  1. {output_dir}/fig1_1d_fdtd.png")
print(f"  2. {output_dir}/fig2_numerical_dispersion.png")
print(f"  3. {output_dir}/fig3_2d_fdtd.png")
print(f"  4. {output_dir}/fig4_pml_boundary.png")
print(f"  5. {output_dir}/fig5_time_frequency.png")
print(f"  6. {output_dir}/fig6_cfl_stability.png")
print("\n所有仿真图已保存到主题002目录")
print("=" * 60)