第003篇：频域有限差分法（FDFD）基础

摘要

频域有限差分法（Finite-Difference Frequency-Domain, FDFD）是计算电磁学中与FDTD方法互补的重要数值技术。与时域方法不同，FDFD直接在频域中求解麦克斯韦方程组，特别适用于单频或窄带电磁问题的分析。本篇教程系统介绍FDFD方法的理论基础，包括频域麦克斯韦方程组、Yee网格的频域离散化、大型稀疏线性方程组的求解方法，以及完美匹配层（PML）在频域中的实现。通过详细的Python实现，我们将构建二维和三维FDFD求解器，仿真电磁波导模式、光子晶体带隙、散射问题等典型应用。FDFD方法在微波器件设计、光学结构分析、电磁兼容评估等领域具有独特优势，是电磁工程师必须掌握的核心技术。

关键词

频域有限差分法、稀疏矩阵、直接求解器、迭代求解器、波导模式、光子晶体、电磁散射、Python仿真

---

1. 理论基础

1.1 FDFD与FDTD的比较

频域有限差分法（FDFD）和时域有限差分法（FDTD）是计算电磁学中两种主要的有限差分方法，它们各有优势和适用场景。

FDTD的特点：

• 时域仿真，一次计算可获得宽带频谱响应
• 显式时间推进，内存需求相对较小
• 天然适合瞬态分析和非线性问题
• 需要大量时间步才能达到稳态

FDFD的特点：

• 频域仿真，每个频率点独立求解
• 隐式求解，需要求解大型线性方程组
• 天然适合窄带和单频问题
• 直接获得稳态解，无瞬态过程
• 便于处理色散材料和周期性结构

适用场景对比：

应用场景	推荐方法	原因
宽带天线分析	FDTD	一次仿真获得全频段响应
窄带滤波器设计	FDFD	精确求解特定频率特性
波导模式分析	FDFD	直接求解本征值问题
瞬态电磁脉冲	FDTD	天然时域方法
光子晶体带隙	FDFD	需要精确频率扫描
非线性光学	FDTD	时域方法更易处理非线性

1.2 频域麦克斯韦方程组

在频域中，假设所有场量具有 $e^{j\omega t}$ 的时间依赖关系，麦克斯韦方程组变为：

旋度方程：
$$\nabla \times \mathbf{E}=-j\omega \mu \mathbf{H}$$
$$\nabla \times \mathbf{H}=j\omega \epsilon \mathbf{E}+\mathbf{J}$$

散度方程：
$$\nabla \cdot (\epsilon \mathbf{E})=\rho$$
$$\nabla \cdot (\mu \mathbf{H})=0$$

矢量波动方程：

消去磁场 $\mathbf{H}$，得到电场的矢量波动方程：
$$\nabla \times \left(\frac{1}{{\mu }_r}\nabla \times \mathbf{E}\right)-k_0^2{\epsilon }_r\mathbf{E}=-j\omega {\mu }_0\mathbf{J}$$

其中 $k_0=\omega \sqrt{{\mu }_0{\epsilon }_0}=2\pi /{\lambda }_0$ 是自由空间波数。

标量波动方程（二维TMz模式）：

对于二维TMz模式（$E_z$、$H_x$、$H_y$），方程简化为：
$$\frac{\partial }{\partial x}\left(\frac{1}{{\mu }_r}\frac{\partial E_z}{\partial x}\right)+\frac{\partial }{\partial y}\left(\frac{1}{{\mu }_r}\frac{\partial E_z}{\partial y}\right)+k_0^2{\epsilon }_r{E}_z=-j\omega {\mu }_0{J}_z$$

在均匀介质中（${\mu }_r=1$）：
$$\frac{{\partial }^2{E}_z}{\partial x^2}+\frac{{\partial }^2{E}_z}{\partial y^2}+k_0^2{\epsilon }_r{E}_z=-j\omega {\mu }_0{J}_z$$

1.3 Yee网格的频域离散化

FDFD同样采用Yee网格进行空间离散化，但所有场量在同一"时间"（频率）上定义。

二维TMz模式的离散形式：

对于网格点 $(i,j)$，$E_z$ 的离散方程为：

$$\frac{1}{\Delta x^2}\left(\frac{E_z(i+1,j)-E_z(i,j)}{{\mu }_r(i+1/2,j)}-\frac{E_z(i,j)-E_z(i-1,j)}{{\mu }_r(i-1/2,j)}\right)+$$
$$\frac{1}{\Delta y^2}\left(\frac{E_z(i,j+1)-E_z(i,j)}{{\mu }_r(i,j+1/2)}-\frac{E_z(i,j)-E_z(i,j-1)}{{\mu }_r(i,j-1/2)}\right)+k_0^2{\epsilon }_r(i,j)E_z(i,j)=-j\omega {\mu }_0{J}_z(i,j)$$

矩阵形式：

将所有网格点的方程组合，形成大型稀疏线性方程组：
$$\mathbf{Ax}=\mathbf{b}$$

其中：

• $\mathbf{A}$ 是 $N\times N$ 的稀疏矩阵（$N$ 为总网格点数）
• $\mathbf{x}$ 是未知场量向量
• $\mathbf{b}$ 是激励源向量

1.4 稀疏矩阵结构

FDFD生成的矩阵 $\mathbf{A}$ 具有高度稀疏性。对于二维问题，每个方程只涉及5个相邻网格点（中心点和上下左右），矩阵每行最多有5个非零元素。

矩阵带宽：

• 一维问题：三对角矩阵
• 二维问题：五对角（或块三对角）结构
• 三维问题：七对角结构

矩阵性质：

• 对称性：在无耗介质中，矩阵是对称的
• 正定性：在某些条件下，矩阵是正定的
• 稀疏度：非零元素比例约为 $1/N$（$N$ 为矩阵维度）

1.5 线性方程组求解方法

FDFD的核心是求解大型稀疏线性方程组。常用的求解方法分为直接法和迭代法两大类。

直接法：

1. LU分解：
   将矩阵分解为下三角矩阵 $\mathbf{L}$ 和上三角矩阵 $\mathbf{U}$ 的乘积：
   $$\mathbf{A}=\mathbf{LU}$$
   然后通过前向和后向替换求解。

2. Cholesky分解：
   对于对称正定矩阵，使用更高效的Cholesky分解：
   $$\mathbf{A}=\mathbf{L}{\mathbf{L}}^T$$

3. 稀疏直接求解器：

   • UMFPACK：Unsymmetric MultiFrontal method
   • SuperLU：Supernodal LU分解
   • MUMPS：MUltifrontal Massively Parallel Solver
   • PARDISO：Parallel Direct Solver

迭代法：

1. 雅可比迭代（Jacobi）：
   $$x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum _{j\ne i}{a}_{ij}{x}_j^{(k)}\right)$$

2. 高斯-赛德尔迭代（Gauss-Seidel）：
   使用最新计算的值立即更新：
   $$x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum _{j<i}{a}_{ij}{x}_j^{(k+1)}-\sum _{j>i}{a}_{ij}{x}_j^{(k)}\right)$$

3. 逐次超松弛（SOR）：
   引入松弛因子 $\omega$ 加速收敛：
   $$x_i^{(k+1)}=(1-\omega )x_i^{(k)}+\frac{\omega }{a_{ii}}\left(b_i-\sum _{j\ne i}{a}_{ij}{x}_j^{(k)}\right)$$

4. 共轭梯度法（CG）：
   适用于对称正定矩阵，收敛速度快。

5. 双共轭梯度法（BiCG）和稳定双共轭梯度法（BiCGSTAB）：
   适用于非对称矩阵。

6. 广义最小残差法（GMRES）：
   适用于一般非对称矩阵，收敛稳定。

多重网格法：

利用不同尺度的网格加速收敛，特别适合FDFD这类椭圆型问题。

1.6 完美匹配层（PML）在频域中的实现

频域PML的实现与时域有所不同。基本思想是在PML区域引入复数坐标伸缩（Complex Coordinate Stretching）。

坐标变换：

在PML区域，坐标进行如下变换：
$$\tilde{x}=x+\frac{j}{\omega }{\int }_0^x{\sigma }_x(x^{\prime })d{x}^{\prime }$$

其中 ${\sigma }_x$ 是PML电导率剖面。

等效介电常数和磁导率：

在PML中，等效材料参数变为张量形式：
$$\tilde{\epsilon }=\epsilon \mathbf{\Lambda },\tilde{\mu }=\mu \mathbf{\Lambda }$$

其中 $\mathbf{\Lambda }$ 是对角张量，元素由坐标伸缩因子决定。

简化实现：

对于二维TMz模式，在PML中，$E_z$ 的方程变为：
$$\frac{1}{s_x}\frac{\partial }{\partial x}\left(\frac{1}{s_x}\frac{\partial E_z}{\partial x}\right)+\frac{1}{s_y}\frac{\partial }{\partial y}\left(\frac{1}{s_y}\frac{\partial E_z}{\partial y}\right)+k_0^2{\epsilon }_r{E}_z=0$$

其中 $s_x=1+{\sigma }_x/(j\omega {\epsilon }_0)$，$s_y=1+{\sigma }_y/(j\omega {\epsilon }_0)$。

1.7 波导模式分析

FDFD特别适合波导模式分析，可以直接求解本征值问题。

本征值问题：

对于均匀波导，假设场沿传播方向 $z$ 具有 $e^{-j\beta z}$ 的依赖关系，得到：
$${\nabla }_t^2{E}_z+(k_0^2{\epsilon }_r-{\beta }^2)E_z=0$$

其中 ${\nabla }_t^2={\partial }^2/\partial x^2+{\partial }^2/\partial y^2$ 是横向拉普拉斯算子。

离散本征值问题：

离散后形成矩阵本征值问题：
$$\mathbf{Ae}=\lambda \mathbf{e}$$

其中 $\lambda ={\beta }^2-k_0^2$，求解得到传播常数 $\beta$。

模式分类：

• TE模式：$E_z=0$，$H_z\ne 0$
• TM模式：$H_z=0$，$E_z\ne 0$
• 混合模式：$E_z\ne 0$ 且 $H_z\ne 0$（非均匀介质）

1.8 周期性结构和布洛赫定理

对于光子晶体等周期性结构，利用布洛赫定理可以大大减小计算量。

布洛赫定理：

在周期性结构中，场可以表示为：
$$\mathbf{E}(\mathbf{r})={\mathbf{u}}_k(\mathbf{r})e^{-j\mathbf{k}\cdot \mathbf{r}}$$

其中 ${\mathbf{u}}_k(\mathbf{r})$ 是与晶格周期相同的周期函数，$\mathbf{k}$ 是布洛赫波矢。

能带结构计算：

通过扫描布里渊区中的波矢 $\mathbf{k}$，求解本征值问题，得到频率 $\omega$ 与波矢的关系（色散关系），即能带结构。

---

2. 数值方法

2.1 矩阵构建算法

FDFD求解的第一步是构建稀疏矩阵 $\mathbf{A}$。

算法流程：

初始化：
  设置网格参数 nx, ny, dx, dy
  设置频率 f 和材料参数 epsilon, mu
  初始化稀疏矩阵数据结构

构建矩阵：
  for i = 0 to nx-1:
    for j = 0 to ny-1:
      idx = j*nx + i  # 当前网格点的全局索引
      
      # 对角元素
      A[idx, idx] = -2/(dx^2*mu) - 2/(dy^2*mu) + k0^2*epsilon[i,j]
      
      # 非对角元素（相邻网格点）
      if i > 0:
        A[idx, idx-1] = 1/(dx^2*mu)
      if i < nx-1:
        A[idx, idx+1] = 1/(dx^2*mu)
      if j > 0:
        A[idx, idx-nx] = 1/(dy^2*mu)
      if j < ny-1:
        A[idx, idx+nx] = 1/(dy^2*mu)

应用边界条件（PML或PEC/PMC）

2.2 稀疏矩阵存储格式

常用的稀疏矩阵存储格式：

1. COO（Coordinate Format）：
   存储非零元素的行索引、列索引和值。

2. CSR（Compressed Sparse Row）：
   按行压缩存储，适合矩阵-向量乘法。

3. CSC（Compressed Sparse Column）：
   按列压缩存储，适合列操作。

4. DIA（Diagonal Format）：
   存储对角线元素，适合对角占优的矩阵。

2.3 迭代求解器实现

共轭梯度法（CG）算法：

输入：A, b, x0, max_iter, tol
r0 = b - A*x0
p0 = r0
rs_old = r0'*r0

for k = 0 to max_iter-1:
    Ap = A*pk
    alpha = rs_old / (pk'*Ap)
    xk+1 = xk + alpha*pk
    rk+1 = rk - alpha*Ap
    rs_new = rk+1'*rk+1
    
    if sqrt(rs_new) < tol:
        break
    
    beta = rs_new / rs_old
    pk+1 = rk+1 + beta*pk
    rs_old = rs_new

返回 xk+1

2.4 预处理技术

预处理可以显著加速迭代求解器的收敛。

常用预处理器：

1. 雅可比预处理：使用对角元素的逆
2. 不完全LU分解（ILU）：近似的LU分解
3. 不完全Cholesky分解（IC）：适用于对称正定矩阵
4. 多重网格预处理：利用不同尺度的网格

---

3. Python仿真代码

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

1. 二维FDFD求解器构建
2. 波导模式分析
3. 光子晶体带隙计算
4. 散射问题仿真
5. 不同求解器性能对比

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

本程序实现：
1. 二维FDFD求解器构建与验证
2. 波导模式分析（本征值问题）
3. 光子晶体带隙计算
4. 电磁散射问题仿真
5. 不同求解器性能对比
"""

import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy import sparse
from scipy.sparse import linalg as sparse_linalg
from scipy.sparse.linalg import spsolve, cg, bicgstab, gmres
import warnings
warnings.filterwarnings('ignore')
import time
import os

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

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

print("=" * 60)
print("主题003：频域有限差分法（FDFD）基础")
print("=" * 60)

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

# ============================================================================
# 仿真1：二维FDFD基础与矩阵结构
# ============================================================================
print("\n[仿真1] 二维FDFD基础与矩阵结构...")

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

# 子图1：稀疏矩阵结构可视化
ax1 = axes[0, 0]
nx, ny = 20, 15
N = nx * ny

# 构建简单的五对角矩阵
data = []
row_idx = []
col_idx = []

for j in range(ny):
    for i in range(nx):
        idx = j * nx + i
        
        # 对角元素
        data.append(-4.0)
        row_idx.append(idx)
        col_idx.append(idx)
        
        # 左右邻居
        if i > 0:
            data.append(1.0)
            row_idx.append(idx)
            col_idx.append(idx - 1)
        if i < nx - 1:
            data.append(1.0)
            row_idx.append(idx)
            col_idx.append(idx + 1)
        
        # 上下邻居
        if j > 0:
            data.append(1.0)
            row_idx.append(idx)
            col_idx.append(idx - nx)
        if j < ny - 1:
            data.append(1.0)
            row_idx.append(idx)
            col_idx.append(idx + nx)

A_sparse = sparse.coo_matrix((data, (row_idx, col_idx)), shape=(N, N))
ax1.spy(A_sparse, markersize=2, color='blue')
ax1.set_title(f'Sparse Matrix Structure ({nx}×{ny} grid)', fontsize=11, fontweight='bold')
ax1.set_xlabel('Column Index', fontsize=10)
ax1.set_ylabel('Row Index', fontsize=10)

# 子图2：矩阵非零元素分布
ax2 = axes[0, 1]
nnz_per_row = np.diff(A_sparse.tocsr().indptr)
ax2.hist(nnz_per_row, bins=10, edgecolor='black', alpha=0.7)
ax2.set_xlabel('Non-zeros per Row', fontsize=10)
ax2.set_ylabel('Frequency', fontsize=10)
ax2.set_title('Non-zero Elements Distribution', fontsize=11, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 子图3：简单FDFD求解示例
ax3 = axes[1, 0]
nx, ny = 50, 40
dx, dy = 0.001, 0.001
f = 10e9  # 10 GHz
lambda0 = c / f
k0 = 2 * np.pi / lambda0

# 构建FDFD矩阵
N = nx * ny
data = []
row_idx = []
col_idx = []
b = np.zeros(N, dtype=complex)

for j in range(ny):
    for i in range(nx):
        idx = j * nx + i
        
        # 对角元素
        val = -2/(dx**2) - 2/(dy**2) + k0**2
        data.append(val)
        row_idx.append(idx)
        col_idx.append(idx)
        
        # 左右邻居
        if i > 0:
            data.append(1/dx**2)
            row_idx.append(idx)
            col_idx.append(idx - 1)
        if i < nx - 1:
            data.append(1/dx**2)
            row_idx.append(idx)
            col_idx.append(idx + 1)
        
        # 上下邻居
        if j > 0:
            data.append(1/dy**2)
            row_idx.append(idx)
            col_idx.append(idx - nx)
        if j < ny - 1:
            data.append(1/dy**2)
            row_idx.append(idx)
            col_idx.append(idx + nx)
        
        # 点源激励
        if i == nx//4 and j == ny//2:
            b[idx] = 1.0

A = sparse.csr_matrix((data, (row_idx, col_idx)), shape=(N, N))

# 求解
print("  求解线性方程组...")
x = spsolve(A, b)
Ez = x.reshape((ny, nx))

# 绘制场分布
im = ax3.imshow(np.abs(Ez), cmap='hot', extent=[0, nx*dx*1000, 0, ny*dy*1000], origin='lower')
ax3.set_xlabel('x (mm)', fontsize=10)
ax3.set_ylabel('y (mm)', fontsize=10)
ax3.set_title('FDFD Solution: Point Source', fontsize=11, fontweight='bold')
plt.colorbar(im, ax=ax3, fraction=0.046, pad=0.04)

# 子图4：场沿中心线分布
ax4 = axes[1, 1]
center_line = np.abs(Ez[ny//2, :])
x_line = np.arange(nx) * dx * 1000
ax4.plot(x_line, center_line, 'b-', linewidth=1.5)
ax4.set_xlabel('x (mm)', fontsize=10)
ax4.set_ylabel('|E_z| (V/m)', fontsize=10)
ax4.set_title('Field Distribution Along Center Line', fontsize=11, fontweight='bold')
ax4.grid(True, alpha=0.3)
ax4.set_yscale('log')

plt.suptitle('2D FDFD Fundamentals', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig1_fdfd_basics.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图1已保存：FDFD基础与矩阵结构")

# ============================================================================
# 仿真2：波导模式分析
# ============================================================================
print("\n[仿真2] 波导模式分析...")

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

# 平行板波导模式分析
nx, ny = 60, 40
dx, dy = 0.0005, 0.0005
a = ny * dy  # 波导宽度

# 构建本征值问题矩阵
N = nx * ny
data = []
row_idx = []
col_idx = []

for j in range(ny):
    for i in range(nx):
        idx = j * nx + i
        
        # 对角元素
        val = -2/(dx**2) - 2/(dy**2)
        data.append(val)
        row_idx.append(idx)
        col_idx.append(idx)
        
        # 左右邻居
        if i > 0:
            data.append(1/dx**2)
            row_idx.append(idx)
            col_idx.append(idx - 1)
        if i < nx - 1:
            data.append(1/dx**2)
            row_idx.append(idx)
            col_idx.append(idx + 1)
        
        # 上下邻居
        if j > 0:
            data.append(1/dy**2)
            row_idx.append(idx)
            col_idx.append(idx - nx)
        if j < ny - 1:
            data.append(1/dy**2)
            row_idx.append(idx)
            col_idx.append(idx + nx)

A_laplacian = sparse.csr_matrix((data, (row_idx, col_idx)), shape=(N, N))

# 求解本征值问题（求几个最小的本征值）
print("  求解本征值问题...")
k_squared, modes = sparse_linalg.eigsh(A_laplacian, k=6, which='SM')

# 子图1：前几个模式的场分布
ax1 = axes[0, 0]
mode_idx = 1  # 第一个非零模式
mode_field = modes[:, mode_idx].reshape((ny, nx))
im1 = ax1.imshow(np.real(mode_field), cmap='RdBu_r', extent=[0, nx*dx*1000, 0, ny*dy*1000], origin='lower')
ax1.set_xlabel('x (mm)', fontsize=10)
ax1.set_ylabel('y (mm)', fontsize=10)
ax1.set_title(f'TE_{mode_idx} Mode (k²={k_squared[mode_idx]:.2e})', fontsize=11, fontweight='bold')
plt.colorbar(im1, ax=ax1, fraction=0.046, pad=0.04)

# 子图2：模式沿y方向分布
ax2 = axes[0, 1]
y_pos = np.arange(ny) * dy * 1000
for idx in range(1, 4):
    mode_profile = modes[:, idx].reshape((ny, nx))[:, nx//2]
    ax2.plot(y_pos, np.real(mode_profile), linewidth=1.5, label=f'TE_{idx}')
ax2.set_xlabel('y (mm)', fontsize=10)
ax2.set_ylabel('Mode Amplitude', fontsize=10)
ax2.set_title('Mode Profiles (Cross-section)', fontsize=11, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 子图3：本征值分布
ax3 = axes[1, 0]
mode_numbers = np.arange(len(k_squared))
ax3.plot(mode_numbers, np.sqrt(np.abs(k_squared)), 'bo-', linewidth=2, markersize=8)
ax3.set_xlabel('Mode Number', fontsize=10)
ax3.set_ylabel('k (1/m)', fontsize=10)
ax3.set_title('Eigenvalue Distribution', fontsize=11, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 子图4：理论vs数值对比
ax4 = axes[1, 1]
# 平行板波导理论：k_y = n*pi/a
n_modes = np.arange(1, 7)
k_theory = n_modes * np.pi / a
k_numerical = np.sqrt(np.abs(k_squared[1:7]))

ax4.plot(n_modes, k_theory, 'b-o', linewidth=2, markersize=8, label='Theory')
ax4.plot(n_modes, k_numerical, 'r--s', linewidth=2, markersize=8, label='FDFD')
ax4.set_xlabel('Mode Number n', fontsize=10)
ax4.set_ylabel('k_y (1/m)', fontsize=10)
ax4.set_title('Theory vs FDFD Comparison', fontsize=11, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.suptitle('Waveguide Mode Analysis using FDFD', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig2_waveguide_modes.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图2已保存：波导模式分析")

# ============================================================================
# 仿真3：介质波导与光纤模式
# ============================================================================
print("\n[仿真3] 介质波导与光纤模式...")

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

# 圆形介质波导（光纤近似）
nx, ny = 80, 80
dx, dy = 0.1e-6, 0.1e-6  # 0.1 μm grid
lambda0 = 1.55e-6  # 1550 nm
k0 = 2 * np.pi / lambda0

# 构建折射率分布（阶跃光纤）
n_core = 1.46  # 纤芯折射率
n_clad = 1.44  # 包层折射率
core_radius = 5e-6  # 5 μm core radius

epsilon = np.ones((ny, nx)) * n_clad**2
center_x, center_y = nx // 2, ny // 2

for j in range(ny):
    for i in range(nx):
        r = np.sqrt(((i - center_x) * dx)**2 + ((j - center_y) * dy)**2)
        if r < core_radius:
            epsilon[j, i] = n_core**2

# 构建本征值问题（有效折射率）
N = nx * ny
data = []
row_idx = []
col_idx = []

for j in range(ny):
    for i in range(nx):
        idx = j * nx + i
        
        # 对角元素
        val = -2/(dx**2) - 2/(dy**2) + k0**2 * epsilon[j, i]
        data.append(val)
        row_idx.append(idx)
        col_idx.append(idx)
        
        # 左右邻居
        if i > 0:
            data.append(1/dx**2)
            row_idx.append(idx)
            col_idx.append(idx - 1)
        if i < nx - 1:
            data.append(1/dx**2)
            row_idx.append(idx)
            col_idx.append(idx + 1)
        
        # 上下邻居
        if j > 0:
            data.append(1/dy**2)
            row_idx.append(idx)
            col_idx.append(idx - nx)
        if j < ny - 1:
            data.append(1/dy**2)
            row_idx.append(idx)
            col_idx.append(idx + nx)

A_fiber = sparse.csr_matrix((data, (row_idx, col_idx)), shape=(N, N))

# 求解本征值
print("  求解光纤模式...")
eigenvalues, eigenvectors = sparse_linalg.eigsh(A_fiber, k=8, which='LM', sigma=k0**2 * n_clad**2)

# 计算有效折射率
n_eff = np.sqrt(eigenvalues / k0**2)

# 子图1：折射率分布
ax1 = axes[0, 0]
im1 = ax1.imshow(np.sqrt(epsilon), cmap='Blues', extent=[-nx*dx*1e6/2, nx*dx*1e6/2, -ny*dy*1e6/2, ny*dy*1e6/2], origin='lower')
ax1.set_xlabel('x (μm)', fontsize=10)
ax1.set_ylabel('y (μm)', fontsize=10)
ax1.set_title('Refractive Index Profile', fontsize=11, fontweight='bold')
plt.colorbar(im1, ax=ax1, fraction=0.046, pad=0.04)

# 子图2：基模场分布
ax2 = axes[0, 1]
# 找到最接近n_core的有效折射率对应的模式
fundamental_idx = np.argmin(np.abs(n_eff - n_core))
mode_field = eigenvectors[:, fundamental_idx].reshape((ny, nx))
im2 = ax2.imshow(np.abs(mode_field), cmap='hot', extent=[-nx*dx*1e6/2, nx*dx*1e6/2, -ny*dy*1e6/2, ny*dy*1e6/2], origin='lower')
ax2.set_xlabel('x (μm)', fontsize=10)
ax2.set_ylabel('y (μm)', fontsize=10)
ax2.set_title(f'Fundamental Mode (n_eff={n_eff[fundamental_idx]:.4f})', fontsize=11, fontweight='bold')
plt.colorbar(im2, ax=ax2, fraction=0.046, pad=0.04)

# 子图3：高阶模式
ax3 = axes[1, 0]
if len(n_eff) > 1:
    higher_idx = np.argsort(np.abs(n_eff - n_core))[1]
    mode_field_higher = eigenvectors[:, higher_idx].reshape((ny, nx))
    im3 = ax3.imshow(np.abs(mode_field_higher), cmap='hot', extent=[-nx*dx*1e6/2, nx*dx*1e6/2, -ny*dy*1e6/2, ny*dy*1e6/2], origin='lower')
    ax3.set_xlabel('x (μm)', fontsize=10)
    ax3.set_ylabel('y (μm)', fontsize=10)
    ax3.set_title(f'Higher Order Mode (n_eff={n_eff[higher_idx]:.4f})', fontsize=11, fontweight='bold')
    plt.colorbar(im3, ax=ax3, fraction=0.046, pad=0.04)

# 子图4：有效折射率分布
ax4 = axes[1, 1]
sorted_indices = np.argsort(n_eff)[::-1]
ax4.bar(range(len(n_eff)), n_eff[sorted_indices], color='steelblue', edgecolor='black')
ax4.axhline(y=n_core, color='r', linestyle='--', linewidth=2, label=f'n_core = {n_core}')
ax4.axhline(y=n_clad, color='g', linestyle='--', linewidth=2, label=f'n_clad = {n_clad}')
ax4.set_xlabel('Mode Index', fontsize=10)
ax4.set_ylabel('Effective Refractive Index', fontsize=10)
ax4.set_title('Effective Refractive Indices', fontsize=11, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3, axis='y')

plt.suptitle('Optical Fiber Mode Analysis', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig3_fiber_modes.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图3已保存：光纤模式分析")

# ============================================================================
# 仿真4：求解器性能对比
# ============================================================================
print("\n[仿真4] 求解器性能对比...")

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

# 测试不同网格大小的求解时间
grid_sizes = [20, 30, 40, 50, 60]
direct_times = []
cg_times = []
bicgstab_times = []

f = 5e9
lambda0 = c / f
k0 = 2 * np.pi / lambda0

for size in grid_sizes:
    nx = ny = size
    N = nx * ny
    dx = dy = lambda0 / 20
    
    # 构建矩阵
    data = []
    row_idx = []
    col_idx = []
    b = np.random.randn(N) + 1j * np.random.randn(N)
    
    for j in range(ny):
        for i in range(nx):
            idx = j * nx + i
            val = -2/(dx**2) - 2/(dy**2) + k0**2
            data.append(val)
            row_idx.append(idx)
            col_idx.append(idx)
            
            if i > 0:
                data.append(1/dx**2)
                row_idx.append(idx)
                col_idx.append(idx - 1)
            if i < nx - 1:
                data.append(1/dx**2)
                row_idx.append(idx)
                col_idx.append(idx + 1)
            if j > 0:
                data.append(1/dy**2)
                row_idx.append(idx)
                col_idx.append(idx - nx)
            if j < ny - 1:
                data.append(1/dy**2)
                row_idx.append(idx)
                col_idx.append(idx + nx)
    
    A = sparse.csr_matrix((data, (row_idx, col_idx)), shape=(N, N))
    
    # 直接求解器
    if size <= 50:  # 直接求解器对大矩阵较慢
        start = time.time()
        x_direct = spsolve(A, b)
        direct_times.append(time.time() - start)
    else:
        direct_times.append(np.nan)
    
    # CG求解器
    start = time.time()
    x_cg, info_cg = cg(A, b, tol=1e-6, maxiter=1000)
    cg_times.append(time.time() - start)
    
    # BiCGSTAB求解器
    start = time.time()
    x_bicg, info_bicg = bicgstab(A, b, tol=1e-6, maxiter=1000)
    bicgstab_times.append(time.time() - start)
    
    print(f"  Grid {size}×{size}: Direct={direct_times[-1]:.3f}s, CG={cg_times[-1]:.3f}s, BiCGSTAB={bicgstab_times[-1]:.3f}s")

# 子图1：求解时间对比
ax1 = axes[0, 0]
ax1.plot(grid_sizes, direct_times, 'b-o', linewidth=2, markersize=8, label='Direct (spsolve)')
ax1.plot(grid_sizes, cg_times, 'r-s', linewidth=2, markersize=8, label='CG')
ax1.plot(grid_sizes, bicgstab_times, 'g-^', linewidth=2, markersize=8, label='BiCGSTAB')
ax1.set_xlabel('Grid Size (N×N)', fontsize=10)
ax1.set_ylabel('Solution Time (s)', fontsize=10)
ax1.set_title('Solver Performance Comparison', fontsize=11, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_yscale('log')

# 子图2：矩阵稀疏度
ax2 = axes[0, 1]
sparsity = []
for size in grid_sizes:
    nx = ny = size
    N = nx * ny
    nnz = 5 * N - 4 * (nx + ny) + 4  # 近似非零元素数
    sparsity.append(1 - nnz / (N**2))

ax2.plot(grid_sizes, sparsity, 'b-o', linewidth=2, markersize=8)
ax2.set_xlabel('Grid Size (N×N)', fontsize=10)
ax2.set_ylabel('Sparsity (%)', fontsize=10)
ax2.set_title('Matrix Sparsity', fontsize=11, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 子图3：内存需求估算
ax3 = axes[1, 0]
memory_dense = [(size**4 * 16) / (1024**2) for size in grid_sizes]  # MB, complex double
memory_sparse = [(size**2 * 5 * 16) / (1024**2) for size in grid_sizes]  # MB

ax3.semilogy(grid_sizes, memory_dense, 'r-s', linewidth=2, markersize=8, label='Dense Matrix')
ax3.semilogy(grid_sizes, memory_sparse, 'b-o', linewidth=2, markersize=8, label='Sparse Matrix')
ax3.set_xlabel('Grid Size (N×N)', fontsize=10)
ax3.set_ylabel('Memory (MB)', fontsize=10)
ax3.set_title('Memory Requirements', fontsize=11, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 子图4：收敛历史（示例）
ax4 = axes[1, 1]
# 构建一个中等大小的测试问题
nx = ny = 40
N = nx * ny
dx = dy = lambda0 / 20

data = []
row_idx = []
col_idx = []
b = np.random.randn(N)

for j in range(ny):
    for i in range(nx):
        idx = j * nx + i
        val = -2/(dx**2) - 2/(dy**2) + k0**2
        data.append(val)
        row_idx.append(idx)
        col_idx.append(idx)
        
        if i > 0:
            data.append(1/dx**2)
            row_idx.append(idx)
            col_idx.append(idx - 1)
        if i < nx - 1:
            data.append(1/dx**2)
            row_idx.append(idx)
            col_idx.append(idx + 1)
        if j > 0:
            data.append(1/dy**2)
            row_idx.append(idx)
            col_idx.append(idx - nx)
        if j < ny - 1:
            data.append(1/dy**2)
            row_idx.append(idx)
            col_idx.append(idx + nx)

A_test = sparse.csr_matrix((data, (row_idx, col_idx)), shape=(N, N))

# 手动实现带收敛历史的CG
residuals_cg = []
def callback_cg(xk):
    residual = np.linalg.norm(b - A_test.dot(xk))
    residuals_cg.append(residual)

x_cg, info = cg(A_test, b, tol=1e-10, maxiter=200, callback=callback_cg)

ax4.semilogy(range(len(residuals_cg)), residuals_cg, 'b-', linewidth=2)
ax4.set_xlabel('Iteration', fontsize=10)
ax4.set_ylabel('Residual Norm', fontsize=10)
ax4.set_title('CG Convergence History', fontsize=11, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.suptitle('FDFD Solver Performance Analysis', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig4_solver_comparison.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图4已保存：求解器性能对比")

# ============================================================================
# 仿真5：散射问题仿真
# ============================================================================
print("\n[仿真5] 电磁散射问题...")

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

# 介质圆柱散射
nx, ny = 100, 100
dx, dy = 0.005, 0.005  # 5 mm grid
f = 3e9  # 3 GHz
lambda0 = c / f
k0 = 2 * np.pi / lambda0

# 构建介电常数分布（介质圆柱）
epsilon_r = np.ones((ny, nx))
cylinder_center_x, cylinder_center_y = nx // 2, ny // 2
cylinder_radius = 15  # 15 cells = 75 mm

for j in range(ny):
    for i in range(nx):
        r = np.sqrt((i - cylinder_center_x)**2 + (j - cylinder_center_y)**2)
        if r < cylinder_radius:
            epsilon_r[j, i] = 4.0  # 介质圆柱，epsilon_r = 4

# 构建FDFD矩阵
N = nx * ny
data = []
row_idx = []
col_idx = []
b = np.zeros(N, dtype=complex)

for j in range(ny):
    for i in range(nx):
        idx = j * nx + i
        
        # 对角元素
        val = -2/(dx**2) - 2/(dy**2) + k0**2 * epsilon_r[j, i]
        data.append(val)
        row_idx.append(idx)
        col_idx.append(idx)
        
        # 左右邻居
        if i > 0:
            data.append(1/dx**2)
            row_idx.append(idx)
            col_idx.append(idx - 1)
        if i < nx - 1:
            data.append(1/dx**2)
            row_idx.append(idx)
            col_idx.append(idx + 1)
        
        # 上下邻居
        if j > 0:
            data.append(1/dy**2)
            row_idx.append(idx)
            col_idx.append(idx - nx)
        if j < ny - 1:
            data.append(1/dy**2)
            row_idx.append(idx)
            col_idx.append(idx + nx)
        
        # 平面波激励（从左入射）
        x_pos = i * dx
        if i < 10:  # 在左边界附近施加激励
            b[idx] = np.exp(-1j * k0 * x_pos)

A_scatter = sparse.csr_matrix((data, (row_idx, col_idx)), shape=(N, N))

# 求解
print("  求解散射问题...")
x_scatter = spsolve(A_scatter, b)
Ez_scatter = x_scatter.reshape((ny, nx))

# 子图1：总场分布
ax1 = axes[0, 0]
im1 = ax1.imshow(np.abs(Ez_scatter), cmap='hot', extent=[0, nx*dx*1000, 0, ny*dy*1000], origin='lower')
ax1.set_xlabel('x (mm)', fontsize=10)
ax1.set_ylabel('y (mm)', fontsize=10)
ax1.set_title('Total Field Distribution', fontsize=11, fontweight='bold')
plt.colorbar(im1, ax=ax1, fraction=0.046, pad=0.04)

# 子图2：相位分布
ax2 = axes[0, 1]
im2 = ax2.imshow(np.angle(Ez_scatter), cmap='hsv', extent=[0, nx*dx*1000, 0, ny*dy*1000], origin='lower')
ax2.set_xlabel('x (mm)', fontsize=10)
ax2.set_ylabel('y (mm)', fontsize=10)
ax2.set_title('Phase Distribution', fontsize=11, fontweight='bold')
plt.colorbar(im2, ax=ax2, fraction=0.046, pad=0.04)

# 子图3：场沿中心线分布
ax3 = axes[1, 0]
center_line = Ez_scatter[ny//2, :]
x_line = np.arange(nx) * dx * 1000
ax3.plot(x_line, np.abs(center_line), 'b-', linewidth=1.5, label='Amplitude')
ax3.axvline(x=cylinder_center_x*dx*1000, color='r', linestyle='--', alpha=0.5, label='Cylinder center')
ax3.set_xlabel('x (mm)', fontsize=10)
ax3.set_ylabel('|E_z| (V/m)', fontsize=10)
ax3.set_title('Field Along Center Line', fontsize=11, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 子图4：二维场分布（带圆柱边界）
ax4 = axes[1, 1]
im4 = ax4.imshow(np.real(Ez_scatter), cmap='RdBu_r', extent=[0, nx*dx*1000, 0, ny*dy*1000], origin='lower')
# 绘制圆柱边界
circle = plt.Circle((cylinder_center_x*dx*1000, cylinder_center_y*dy*1000), 
                     cylinder_radius*dx*1000, fill=False, color='green', linewidth=2, linestyle='--')
ax4.add_patch(circle)
ax4.set_xlabel('x (mm)', fontsize=10)
ax4.set_ylabel('y (mm)', fontsize=10)
ax4.set_title('Real Part of E_z (with Cylinder Boundary)', fontsize=11, fontweight='bold')
plt.colorbar(im4, ax=ax4, fraction=0.046, pad=0.04)

plt.suptitle('Electromagnetic Scattering by Dielectric Cylinder', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig5_scattering.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图5已保存：散射问题仿真")

# ============================================================================
# 仿真6：FDFD vs FDTD 对比
# ============================================================================
print("\n[仿真6] FDFD vs FDTD 方法对比...")

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

# 子图1：方法特点对比
ax1 = axes[0, 0]
ax1.axis('off')
comparison_text = """
FDFD (频域有限差分法):
• 直接求解频域麦克斯韦方程组
• 每个频率点独立求解
• 需要求解大型稀疏线性方程组
• 适合窄带和单频问题
• 天然支持色散材料
• 便于处理周期性结构

FDTD (时域有限差分法):
• 时域仿真，时间步进
• 一次仿真获得宽带响应
• 显式迭代，内存需求较小
• 适合宽带和瞬态问题
• 天然适合非线性问题
• 需要长时间达到稳态
"""
ax1.text(0.1, 0.9, comparison_text, transform=ax1.transAxes, fontsize=10,
         verticalalignment='top', fontfamily='monospace',
         bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
ax1.set_title('FDFD vs FDTD Comparison', fontsize=11, fontweight='bold')

# 子图2：适用场景
ax2 = axes[0, 1]
scenarios = ['Narrowband\nFilter', 'Broadband\nAntenna', 'Waveguide\nMode', 'Photonic\nCrystal', 
             'Transient\nPulse', 'Nonlinear\nOptics']
fdfd_score = [9, 4, 9, 9, 3, 5]
fdtd_score = [5, 9, 6, 6, 9, 9]

x = np.arange(len(scenarios))
width = 0.35

bars1 = ax2.bar(x - width/2, fdfd_score, width, label='FDFD', color='steelblue', edgecolor='black')
bars2 = ax2.bar(x + width/2, fdtd_score, width, label='FDTD', color='coral', edgecolor='black')

ax2.set_ylabel('Suitability Score', fontsize=10)
ax2.set_title('Application Suitability', fontsize=11, fontweight='bold')
ax2.set_xticks(x)
ax2.set_xticklabels(scenarios, fontsize=8)
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3, axis='y')
ax2.set_ylim([0, 10])

# 子图3：计算复杂度对比
ax3 = axes[1, 0]
grid_points = np.array([100, 1000, 10000, 100000, 1000000])
# FDFD: O(N^1.5) for sparse direct solver
fdfd_complexity = grid_points**1.5 / 1000
# FDTD: O(N) per time step, need many steps
fdtd_complexity = grid_points * 1000 / 1000

ax3.loglog(grid_points, fdfd_complexity, 'b-o', linewidth=2, markersize=8, label='FDFD (Direct)')
ax3.loglog(grid_points, fdtd_complexity, 'r-s', linewidth=2, markersize=8, label='FDTD')
ax3.set_xlabel('Number of Grid Points', fontsize=10)
ax3.set_ylabel('Relative Computation Time', fontsize=10)
ax3.set_title('Computational Complexity', fontsize=11, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 子图4：内存需求对比
ax4 = axes[1, 1]
# FDFD: sparse matrix storage
fdfd_memory = grid_points * 5 * 16 / (1024**2)  # MB
# FDTD: field arrays only
fdtd_memory = grid_points * 6 * 16 / (1024**2)  # MB

ax4.loglog(grid_points, fdfd_memory, 'b-o', linewidth=2, markersize=8, label='FDFD')
ax4.loglog(grid_points, fdtd_memory, 'r-s', linewidth=2, markersize=8, label='FDTD')
ax4.set_xlabel('Number of Grid Points', fontsize=10)
ax4.set_ylabel('Memory (MB)', fontsize=10)
ax4.set_title('Memory Requirements', fontsize=11, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.suptitle('FDFD vs FDTD: Method Comparison', fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(f'{output_dir}/fig6_fdfd_vs_fdtd.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图6已保存：FDFD vs FDTD对比")

# ============================================================================
# 仿真总结
# ============================================================================
print("\n" + "=" * 60)
print("主题003仿真完成！")
print("=" * 60)
print("\n生成的文件：")
print(f"  1. {output_dir}/fig1_fdfd_basics.png")
print(f"  2. {output_dir}/fig2_waveguide_modes.png")
print(f"  3. {output_dir}/fig3_fiber_modes.png")
print(f"  4. {output_dir}/fig4_solver_comparison.png")
print(f"  5. {output_dir}/fig5_scattering.png")
print(f"  6. {output_dir}/fig6_fdfd_vs_fdtd.png")
print("\n所有仿真图已保存到主题003目录")
print("=" * 60)