电磁场仿真-主题079-时域有限元法

kkchenjj

7人浏览 · 2026-03-25 21:28:22

kkchenjj · 2026-03-25 21:28:22 发布

主题079: 时域有限元法 (Time-Domain Finite Element Method, TD-FEM)

1. 引言

1.1 什么是时域有限元法

时域有限元法(Time-Domain Finite Element Method, TD-FEM)是将有限元法的空间离散优势与时域分析相结合的一种数值方法。与传统的频域有限元法(Frequency-Domain FEM)不同，TD-FEM直接在时间域中求解麦克斯韦方程组，能够一次性获得宽带频率响应，特别适用于瞬态电磁场分析、宽带散射问题和色散材料建模。

1.2 为什么需要时域方法

频域有限元法虽然成熟且高效，但在某些应用场景下存在局限性：

宽带分析的挑战：

频域方法需要对每个频率点单独求解
对于宽带问题，需要计算数十甚至数百个频点
总计算时间与频点数量成正比

瞬态现象的研究：

雷击、静电放电等瞬态过程
脉冲雷达和超宽带通信
电磁脉冲(EMP)效应分析

非线性材料：

频域方法难以处理时变或非线性材料
时域方法天然支持材料非线性

色散材料：

实际材料介电常数随频率变化
时域方法可以直接建模色散特性

1.3 TD-FEM的优势与挑战

优势：

单次仿真获得宽带响应：通过傅里叶变换，一次时域仿真可获得全频段信息
复杂几何适应性：继承有限元法的非结构化网格优势
色散材料建模：支持Debye、Lorentz、Drude等色散模型
非线性分析：可扩展到时变和非线性材料
物理直观性：直接观察波的传播和相互作用过程

挑战：

计算资源需求：需要存储多个时间步的场分布
数值稳定性：时间步长受CFL条件限制
长时间仿真：谐振结构需要很长的仿真时间
色散误差：数值色散可能影响高精度应用

1.4 应用领域

TD-FEM在以下领域有广泛应用：

天线与微波工程：

超宽带天线设计
瞬态辐射特性分析
天线互耦评估

电磁兼容(EMC)：

高速电路信号完整性
电磁干扰(EMI)分析
屏蔽效能评估

雷达与散射：

目标瞬态响应
雷达散射截面(RCS)宽带计算
穿墙雷达

生物电磁学：

人体组织电磁暴露评估
磁共振成像(MRI)线圈设计
深部脑刺激(DBS)优化

光电子学：

光子晶体分析
表面等离子体激元
超材料设计

2. 时域有限元法基础理论

2.1 麦克斯韦方程组的时域形式

时域麦克斯韦方程组描述了电磁场随时间和空间的演化：

法拉第电磁感应定律：
$\nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}$

安培环路定律：
$\nabla \times \mathbf{H} = \mathbf{J} + \frac{\partial \mathbf{D}}{\partial t}$

高斯电场定律：
$\nabla \cdot \mathbf{D} = \rho$

高斯磁场定律：
$\nabla \cdot \mathbf{B} = 0$

其中， $E\mathbf{E}$ 为电场强度， $H\mathbf{H}$ 为磁场强度， $D\mathbf{D}$ 为电位移矢量， $B\mathbf{B}$ 为磁感应强度， $J\mathbf{J}$ 为电流密度， $ρ\rho$ 为电荷密度。

2.2 波动方程的推导

对于线性、各向同性、非色散介质，本构关系为：

$\mathbf{D} = \varepsilon \mathbf{E}, \quad \mathbf{B} = \mu \mathbf{H}, \quad \mathbf{J} = \sigma \mathbf{E}$

在无源区域( $J=0,ρ=0\mathbf{J}=0, \rho=0$ )，对法拉第定律取旋度：

$\nabla \times (\nabla \times \mathbf{E}) = -\mu \frac{\partial}{\partial t}(\nabla \times \mathbf{H}) = -\mu\varepsilon \frac{\partial^2 \mathbf{E}}{\partial t^2}$

利用矢量恒等式 $∇×(∇×E)=∇(∇⋅E)−∇2E\nabla \times (\nabla \times \mathbf{E}) = \nabla(\nabla \cdot \mathbf{E}) - \nabla^2 \mathbf{E}$ ，在无源区域 $∇⋅E=0\nabla \cdot \mathbf{E} = 0$ ，得到波动方程：

$\nabla^2 \mathbf{E} - \mu\varepsilon \frac{\partial^2 \mathbf{E}}{\partial t^2} = 0$

或写成：

$\nabla^2 \mathbf{E} - \frac{1}{c^2} \frac{\partial^2 \mathbf{E}}{\partial t^2} = 0$

其中 $1/\sqrt{\mu\varepsilon}$ 为介质中的光速。

2.3 二维波动方程

对于二维问题(如 $TM_z$ 模式， $Ez≠0,Hx=Hy=0E_z \neq 0, H_x = H_y = 0$ )，波动方程简化为标量形式：

$\frac{\partial^2 E_z}{\partial x^2} + \frac{\partial^2 E_z}{\partial y^2} - \mu\varepsilon \frac{\partial^2 E_z}{\partial t^2} = 0$

这是TD-FEM求解的基本方程。

2.4 有限元空间离散

采用伽辽金方法，将求解域 $Ω\Omega$ 离散为有限个单元 $Ωe\Omega_e$ 。场量用节点基函数展开：

$E_z(x, y, t) = \sum_{j=1}^{N} N_j(x, y) u_j(t)$

其中 $N_j$ 为节点 $j$ 的形状函数， $u_j(t)$ 为节点 $j$ 处的时变场值。

将展开式代入波动方程，乘以测试函数 $N_i$ 并在单元上积分，得到弱形式：

$\int_{\Omega_e} \left( \nabla N_i \cdot \nabla \sum_j N_j u_j + \mu\varepsilon N_i \sum_j N_j \ddot{u}_j \right) d\Omega = 0$

整理后得到单元方程：

$\mathbf{K}^{(e)} \mathbf{u}^{(e)} + \mathbf{M}^{(e)} \ddot{\mathbf{u}}^{(e)} = 0$

其中：

刚度矩阵：
$K_{ij}^{(e)} = \int_{\Omega_e} \nabla N_i \cdot \nabla N_j \, d\Omega$

质量矩阵：
$M_{ij}^{(e)} = \mu\varepsilon \int_{\Omega_e} N_i N_j \, d\Omega$

2.5 全局系统组装

将所有单元的贡献组装到全局矩阵：

$\mathbf{K} \mathbf{u} + \mathbf{M} \ddot{\mathbf{u}} = \mathbf{f}$

其中 $f\mathbf{f}$ 为激励源向量。

考虑损耗时，还需引入阻尼矩阵：

$\mathbf{M} \ddot{\mathbf{u}} + \mathbf{C} \dot{\mathbf{u}} + \mathbf{K} \mathbf{u} = \mathbf{f}$

阻尼矩阵：
$C_{ij}^{(e)} = \sigma \int_{\Omega_e} N_i N_j \, d\Omega$

2.6 Python实现：二维时域有限元求解器

class TDFEM2D:
    """二维时域有限元求解器"""
    
    def __init__(self, nodes, elements, epsilon_r=1.0, mu_r=1.0, sigma=0.0):
        """
        参数:
            nodes: 节点坐标数组 (N, 2)
            elements: 三角形单元连接数组 (M, 3)
            epsilon_r: 相对介电常数
            mu_r: 相对磁导率
            sigma: 电导率
        """
        self.nodes = nodes
        self.elements = elements
        self.n_nodes = len(nodes)
        self.n_elements = len(elements)
        
        # 材料参数
        self.epsilon0 = 8.854e-12
        self.mu0 = 4 * np.pi * 1e-7
        self.epsilon = epsilon_r * self.epsilon0
        self.mu = mu_r * self.mu0
        self.sigma = sigma
        
        # 光速
        self.c = 1 / np.sqrt(self.epsilon * self.mu)
        
    def compute_element_matrices(self, element_idx):
        """计算单元矩阵"""
        elem = self.elements[element_idx]
        x = self.nodes[elem, 0]
        y = self.nodes[elem, 1]
        
        # 计算单元面积
        area = 0.5 * abs((x[1] - x[0]) * (y[2] - y[0]) - 
                         (x[2] - x[0]) * (y[1] - y[0]))
        
        # 形状函数梯度
        b = np.array([y[1] - y[2], y[2] - y[0], y[0] - y[1]])
        c = np.array([x[2] - x[1], x[0] - x[2], x[1] - x[0]])
        
        # 刚度矩阵 (1/mu * 梯度项)
        K_elem = np.outer(b, b) + np.outer(c, c)
        K_elem = K_elem / (4 * area * self.mu)
        
        # 质量矩阵 (epsilon * 势能项)
        M_elem = self.epsilon * area / 12 * np.array([[2, 1, 1],
                                                       [1, 2, 1],
                                                       [1, 1, 2]])
        
        # 阻尼矩阵 (sigma * 耗散项)
        C_elem = self.sigma * area / 12 * np.array([[2, 1, 1],
                                                     [1, 2, 1],
                                                     [1, 1, 2]])
        
        return K_elem, M_elem, C_elem, area
    
    def assemble_system(self):
        """组装全局系统矩阵"""
        rows, cols = [], []
        data_K, data_M, data_C = [], [], []
        
        for e in range(self.n_elements):
            K_elem, M_elem, C_elem, _ = self.compute_element_matrices(e)
            elem = self.elements[e]
            
            for i in range(3):
                for j in range(3):
                    rows.append(elem[i])
                    cols.append(elem[j])
                    data_K.append(K_elem[i, j])
                    data_M.append(M_elem[i, j])
                    data_C.append(C_elem[i, j])
        
        self.K = sparse.coo_matrix((data_K, (rows, cols)), 
                                    shape=(self.n_nodes, self.n_nodes)).tocsr()
        self.M = sparse.coo_matrix((data_M, (rows, cols)), 
                                    shape=(self.n_nodes, self.n_nodes)).tocsr()
        self.C = sparse.coo_matrix((data_C, (rows, cols)), 
                                    shape=(self.n_nodes, self.n_nodes)).tocsr()
        
        return self.K, self.M, self.C

3. 时间积分方法

3.1 时间离散的必要性

经过空间离散后，我们得到二阶常微分方程组：

$\mathbf{M} \ddot{\mathbf{u}} + \mathbf{C} \dot{\mathbf{u}} + \mathbf{K} \mathbf{u} = \mathbf{f}(t)$

需要选择合适的时间积分方法将其转化为代数方程组。

3.2 时间积分方法的分类

显式方法：

当前时刻的场值仅依赖于过去时刻
条件稳定，时间步长受限
计算效率高，每步计算量小
代表：中心差分法

隐式方法：

当前时刻的场值依赖于自身
无条件稳定(或条件更宽松)
每步需要求解线性方程组
代表：Newmark-beta法、Wilson-theta法

3.3 Newmark-beta方法

Newmark-beta方法是最常用的隐式时间积分方法之一，通过两个参数 $β\beta$ 和 $γ\gamma$ 控制精度和稳定性。

基本假设：

$\dot{\mathbf{u}}_{n+1} = \dot{\mathbf{u}}_n + (1-\gamma)\Delta t \, \ddot{\mathbf{u}}_n + \gamma \Delta t \, \ddot{\mathbf{u}}_{n+1}$

$\mathbf{u}_{n+1} = \mathbf{u}_n + \Delta t \, \dot{\mathbf{u}}_n + \frac{(\Delta t)^2}{2} \left[(1-2\beta)\ddot{\mathbf{u}}_n + 2\beta \ddot{\mathbf{u}}_{n+1}\right]$

预测-校正格式：

预测步：
$\tilde{\mathbf{u}}_{n+1} = \mathbf{u}_n + \Delta t \, \dot{\mathbf{u}}_n + \frac{(\Delta t)^2}{2}(1-2\beta)\ddot{\mathbf{u}}_n$

$\tilde{\dot{\mathbf{u}}}_{n+1} = \dot{\mathbf{u}}_n + (1-\gamma)\Delta t \, \ddot{\mathbf{u}}_n$
求解加速度：
$\left[\mathbf{M} + \gamma \Delta t \, \mathbf{C} + \beta (\Delta t)^2 \mathbf{K}\right] \ddot{\mathbf{u}}_{n+1} = \mathbf{f}_{n+1} - \mathbf{C} \tilde{\dot{\mathbf{u}}}_{n+1} - \mathbf{K} \tilde{\mathbf{u}}_{n+1}$
校正步：
$\mathbf{u}_{n+1} = \tilde{\mathbf{u}}_{n+1} + \beta (\Delta t)^2 \ddot{\mathbf{u}}_{n+1}$

$\dot{\mathbf{u}}_{n+1} = \tilde{\dot{\mathbf{u}}}_{n+1} + \gamma \Delta t \, \ddot{\mathbf{u}}_{n+1}$

参数选择：

$γ=0.5,β=0.25\gamma = 0.5, \beta = 0.25$ ：平均加速度法(无条件稳定)
$γ=0.5,β=0\gamma = 0.5, \beta = 0$ ：中心差分法(条件稳定)
$γ=0.5,β=1/6\gamma = 0.5, \beta = 1/6$ ：线性加速度法

3.4 Python实现：Newmark-beta时间积分

def newmark_timestep(self, u_n, v_n, a_n, dt, beta=0.25, gamma=0.5):
    """
    Newmark-beta时间积分
    
    参数:
        u_n: 当前位移
        v_n: 当前速度
        a_n: 当前加速度
        dt: 时间步长
        beta, gamma: Newmark参数
    """
    # 预测步
    u_pred = u_n + dt * v_n + (dt**2 / 2) * (1 - 2*beta) * a_n
    v_pred = v_n + dt * (1 - gamma) * a_n
    
    # 求解加速度
    A = self.M + gamma * dt * self.C + beta * dt**2 * self.K
    b = -self.K @ u_pred - self.C @ v_pred
    
    a_new = splinalg.spsolve(A, b)
    
    # 校正步
    u_new = u_pred + beta * dt**2 * a_new
    v_new = v_pred + gamma * dt * a_new
    
    return u_new, v_new, a_new

3.5 CFL稳定性条件

对于显式时间积分，时间步长必须满足CFL(Courant-Friedrichs-Lewy)条件：

$\Delta t \leq \frac{1}{c \sqrt{\frac{1}{(\Delta x)^2} + \frac{1}{(\Delta y)^2}}}$

对于二维均匀网格( $Δx=Δy=h\Delta x = \Delta y = h$ )：

$\Delta t \leq \frac{h}{c\sqrt{2}}$

CFL数是衡量时间步长选择的重要指标：

$\text{CFL} = \frac{c \Delta t}{h}$

CFL $≤1/2\leq 1/\sqrt{2}$ ：稳定(显式方法)
CFL $1/\sqrt{2}$ ：不稳定

3.6 数值色散

数值离散会引入色散误差，导致波速与物理波速不同。数值波速 $c^*$ 与网格分辨率相关：

$\frac{c^*}{c} = \frac{\sin^{-1}(\text{CFL} \cdot \sin(kh/2))}{kh/2}$

减小数值色散的方法：

增加网格分辨率(每波长10-20个单元)
使用高阶单元
优化时间步长

4. 脉冲激励源

4.1 脉冲激励的优势

在TD-FEM中，使用脉冲激励而非单频正弦激励具有以下优势：

宽带信息：单个脉冲包含宽频谱成分
单次仿真：一次仿真获得全频段响应
瞬态分析：直接观察波的传播过程
物理直观：易于理解波的相互作用

4.2 高斯脉冲

高斯脉冲是最常用的激励源，其时域表达式为：

$\exp\left[-\left(\frac{t-t_0}{\sigma}\right)^2\right]$

其中 $t_0$ 为脉冲中心时刻， $σ\sigma$ 控制脉冲宽度。

频谱特性：

高斯脉冲的频谱也是高斯型：

$\sqrt{\pi} \sigma \exp\left[-(\pi\sigma f)^2\right] \exp(-j2\pi f t_0)$

-3dB带宽：

$f_{3dB} = \frac{1}{\pi\sigma}$

设计准则：

脉冲宽度 $σ\sigma$ 决定频率分辨率
中心频率 $f0=1/(2πσ)f_0 = 1/(2\pi\sigma)$
有效频带： $\sim 2f_0$

4.3 调制高斯脉冲

调制高斯脉冲在载波频率 $f_c$ 上调制，适合窄带分析：

$\exp\left[-\left(\frac{t-t_0}{\sigma}\right)^2\right] \sin(2\pi f_c t)$

频谱中心位于 $f_c$ ，带宽由 $σ\sigma$ 决定。

4.4 Ricker子波

Ricker子波(二阶导数高斯脉冲)在地震学和雷达中常用：

$\left[1 - 2(\pi f_p (t-t_0))^2\right] \exp\left[-(\pi f_p (t-t_0))^2\right]$

其中 $f_p$ 为峰值频率。Ricker子波在 $f_p$ 处有最大能量，且直流分量为零。

4.5 Python实现：脉冲激励源

class PulsedExcitation:
    """脉冲激励源"""
    
    @staticmethod
    def gaussian_pulse(t, t0, sigma):
        """高斯脉冲"""
        return np.exp(-((t - t0) / sigma)**2)
    
    @staticmethod
    def modulated_gaussian(t, t0, sigma, f0):
        """调制高斯脉冲"""
        return np.exp(-((t - t0) / sigma)**2) * np.sin(2 * np.pi * f0 * t)
    
    @staticmethod
    def derivative_gaussian(t, t0, sigma):
        """高斯脉冲的导数"""
        return -2 * (t - t0) / sigma**2 * np.exp(-((t - t0) / sigma)**2)
    
    @staticmethod
    def ricker_wavelet(t, f0):
        """Ricker子波"""
        t0 = 1.0 / f0
        r = np.pi * f0 * (t - t0)
        return (1 - 2 * r**2) * np.exp(-r**2)
    
    @staticmethod
    def sine_wave(t, f0, t_rise=0):
        """正弦波（带上升沿）"""
        if t < t_rise:
            return np.sin(2 * np.pi * f0 * t) * (t / t_rise)
        else:
            return np.sin(2 * np.pi * f0 * t)

4.6 脉冲参数选择指南

应用场景	脉冲类型	参数建议
宽带RCS	高斯脉冲	$σ=1/(πfmax)\sigma = 1/(\pi f_{max})$
窄带天线	调制高斯	$σ=3/fBW\sigma = 3/f_{BW}$ , $f_c = f_{center}$
探地雷达	Ricker	$f_p = f_{center}$
电路仿真	上升沿正弦	$t_{rise} = 3/f_{max}$

5. 色散材料建模

5.1 色散现象的物理本质

实际材料的电磁参数随频率变化，这种现象称为色散。色散的物理机制包括：

电子极化：

原子中电子云相对核的位移
响应时间： $10−15∼10−1610^{-15} \sim 10^{-16}$ s
影响光学频段

离子极化：

晶格中离子的相对位移
响应时间： $10−12∼10−1310^{-12} \sim 10^{-13}$ s
影响红外频段

取向极化：

极性分子的取向排列
响应时间： $10−6∼10−1210^{-6} \sim 10^{-12}$ s
影响微波频段

空间电荷极化：

界面电荷积累
响应时间： $10−2∼10−610^{-2} \sim 10^{-6}$ s
影响低频段

5.2 Debye色散模型

Debye模型描述具有单一弛豫时间的极性材料：

$\varepsilon(\omega) = \varepsilon_\infty + \frac{\varepsilon_s - \varepsilon_\infty}{1 + j\omega\tau}$

其中：

$εs\varepsilon_s$ ：静态介电常数
$ε∞\varepsilon_\infty$ ：高频介电常数
$τ\tau$ ：弛豫时间

典型应用：水、生物组织、极性液体

水的Debye参数( $20∘20^\circ$ C)：

$εs=80\varepsilon_s = 80$
$ε∞=5\varepsilon_\infty = 5$
$τ=9.4\tau = 9.4$ ps

5.3 Lorentz色散模型

Lorentz模型描述具有共振特性的材料：

$\varepsilon(\omega) = \varepsilon_\infty + \frac{\omega_p^2}{\omega_0^2 - \omega^2 - j\gamma\omega}$

其中：

$ω0\omega_0$ ：共振频率
$ωp\omega_p$ ：等离子体频率
$γ\gamma$ ：阻尼系数

典型应用：光学材料、等离子体、共振超材料

5.4 Drude色散模型

Drude模型描述自由电子气(金属)：

$\varepsilon(\omega) = 1 - \frac{\omega_p^2}{\omega^2 + j\gamma\omega}$

其中：

$ωp=ne2/(m∗ε0)\omega_p = \sqrt{ne^2/(m^*\varepsilon_0)}$ ：等离子体频率
$γ\gamma$ ：碰撞频率
$n$ ：电子密度

典型应用：金属、等离子体、高掺杂半导体

金的Drude参数：

$ωp=2π×2.18\omega_p = 2\pi \times 2.18$ PHz
$γ=2π×6.5\gamma = 2\pi \times 6.5$ THz

5.5 Python实现：色散材料模型

class DispersiveMaterial:
    """色散材料模型"""
    
    def __init__(self, model_type='debye'):
        """
        参数:
            model_type: 'debye', 'lorentz', 'drude'
        """
        self.model_type = model_type
        
    def debye_model(self, omega, epsilon_inf, epsilon_s, tau):
        """
        Debye色散模型
        
        参数:
            omega: 角频率
            epsilon_inf: 高频介电常数
            epsilon_s: 静态介电常数
            tau: 弛豫时间
        """
        epsilon = epsilon_inf + (epsilon_s - epsilon_inf) / (1 + 1j * omega * tau)
        return epsilon
    
    def lorentz_model(self, omega, epsilon_inf, omega_p, gamma, omega_0):
        """
        Lorentz色散模型
        
        参数:
            omega: 角频率
            epsilon_inf: 高频介电常数
            omega_p: 等离子体频率
            gamma: 阻尼系数
            omega_0: 共振频率
        """
        epsilon = epsilon_inf + omega_p**2 / (omega_0**2 - omega**2 - 1j * gamma * omega)
        return epsilon
    
    def drude_model(self, omega, omega_p, gamma):
        """
        Drude色散模型 (金属)
        
        参数:
            omega: 角频率
            omega_p: 等离子体频率
            gamma: 碰撞频率
        """
        epsilon = 1 - omega_p**2 / (omega**2 + 1j * gamma * omega)
        return epsilon

5.6 时域实现：辅助微分方程法

在时域中实现色散材料需要引入辅助微分方程(ADE)。以Debye模型为例：

频域本构关系：

$D(\omega) = \varepsilon(\omega) E(\omega) = \left[\varepsilon_\infty + \frac{\varepsilon_s - \varepsilon_\infty}{1 + j\omega\tau}\right] E(\omega)$

引入极化 $P$ ：

$\varepsilon_\infty E + P$

其中：

$P(\omega) = \frac{\varepsilon_s - \varepsilon_\infty}{1 + j\omega\tau} E(\omega)$

转换到时域：

$\tau \frac{\partial P}{\partial t} + P = (\varepsilon_s - \varepsilon_\infty) E$

与麦克斯韦方程联立求解。

6. 吸收边界条件

6.1 为什么需要吸收边界

开放区域问题(如天线辐射、散射)需要截断无限大空间。吸收边界条件(ABC)用于：

模拟无限大空间：波到达边界后无反射地"消失"
减小计算域：只计算感兴趣区域
避免虚假反射：数值反射会污染解

6.2 完美匹配层(PML)

PML是目前最有效的吸收边界条件，由Berenger于1994年提出。

基本原理：

在计算域外围添加吸收层
在PML中引入人工电导率 $σ\sigma$
波在PML中指数衰减
设计得当可实现零反射

PML中的修正麦克斯韦方程：

$\varepsilon \frac{\partial E_x}{\partial t} + \sigma_y E_x = \frac{\partial H_z}{\partial y}$

$\varepsilon \frac{\partial E_y}{\partial t} + \sigma_x E_y = -\frac{\partial H_z}{\partial x}$

$\mu \frac{\partial H_z}{\partial t} + \sigma^* H_z = \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y}$

其中 $σx\sigma_x$ 和 $σy\sigma_y$ 是坐标相关的电导率。

6.3 PML电导率剖面

电导率从PML内边界到外边界逐渐增加：

$\sigma(d) = \sigma_{max} \left(\frac{d}{\delta}\right)^m$

其中：

$d$ ：到PML内边界的距离
$δ\delta$ ：PML厚度
$m$ ：剖面阶数(通常 $m=3∼4m=3\sim 4$ )
$σmax\sigma_{max}$ ：最大电导率

最优参数选择：

$\sigma_{max} = \frac{m+1}{2\eta\delta}$

其中 $η=μ/ε\eta = \sqrt{\mu/\varepsilon}$ 为波阻抗。

6.4 Python实现：PML边界条件

class AbsorbingBoundary:
    """吸收边界条件"""
    
    def __init__(self, nx, ny, dx, dy, dt, c, thickness=10):
        """
        参数:
            nx, ny: 网格尺寸
            dx, dy: 网格间距
            dt: 时间步长
            c: 波速
            thickness: PML厚度
        """
        self.nx = nx
        self.ny = ny
        self.dx = dx
        self.dy = dy
        self.dt = dt
        self.c = c
        self.thickness = thickness
        
        # 创建PML剖面
        self.sigma_x = self._create_pml_profile(nx, thickness, dx)
        self.sigma_y = self._create_pml_profile(ny, thickness, dy)
        
    def _create_pml_profile(self, n, thickness, d):
        """创建PML电导率剖面"""
        sigma = np.zeros(n)
        
        # 左侧PML
        for i in range(thickness):
            sigma[i] = ((thickness - i) / thickness)**3
        
        # 右侧PML
        for i in range(thickness):
            sigma[n - 1 - i] = ((thickness - i) / thickness)**3
        
        return sigma
    
    def apply_pml(self, u, u_prev):
        """应用PML边界条件"""
        u_pml = u.copy()
        
        # 在PML区域添加衰减
        for i in range(self.ny):
            for j in range(self.nx):
                sigma = self.sigma_x[j] + self.sigma_y[i]
                if sigma > 0:
                    u_pml[i, j] *= np.exp(-sigma * self.dt)
        
        return u_pml

6.5 PML性能评估

反射系数：

理论反射系数：

$\exp\left(-\frac{2\sigma_{max}\delta}{c\varepsilon_0}\right)$

实际反射系数受以下因素影响：

网格离散误差
PML厚度
电导率剖面
入射角度

设计准则：

PML厚度： $8∼168\sim 16$ 个单元
反射系数目标： $< - 40$ dB
对于极端角度入射，需要更厚的PML

7. 瞬态分析技术

7.1 从时域到频域的转换

TD-FEM的优势之一是通过傅里叶变换获得宽带频域响应。

离散傅里叶变换(DFT)：

$\tilde{E}(f) = \sum_{n=0}^{N-1} E(n\Delta t) \exp(-j2\pi f n \Delta t) \Delta t$

快速傅里叶变换(FFT)：
使用FFT算法高效计算DFT，复杂度从 $O(N^2)$ 降至 $O(Nlog⁡N)O(N\log N)$ 。

频率分辨率：

$\Delta f = \frac{1}{T_{sim}} = \frac{1}{N\Delta t}$

其中 $T_{sim}$ 为总仿真时间。

奈奎斯特频率：

$f_{max} = \frac{1}{2\Delta t}$

7.2 包络检测

对于调制信号，需要提取包络：

希尔伯特变换法：

解析信号：

$s_a(t) = s(t) + j\mathcal{H}[s(t)]$

包络：

$|s_a(t)| = \sqrt{s^2(t) + \mathcal{H}^2[s(t)]}$

瞬时相位：

$\phi(t) = \arg[s_a(t)]$

7.3 群延迟计算

群延迟描述信号不同频率成分的传输时间：

$\tau_g(f) = -\frac{1}{2\pi} \frac{d\phi}{df}$

数值计算：

$\tau_g(f) \approx -\frac{\phi_2(f) - \phi_1(f)}{2\pi(f_2 - f_1)}$

7.4 Python实现：瞬态分析工具

class TransientAnalysis:
    """瞬态分析工具"""
    
    def __init__(self, t_history, u_history):
        """
        参数:
            t_history: 时间历史
            u_history: 场量历史
        """
        self.t = t_history
        self.u = u_history
        self.dt = t_history[1] - t_history[0] if len(t_history) > 1 else 1
        
    def compute_spectrum(self, node_idx=None):
        """计算频谱"""
        if node_idx is None:
            # 对所有节点平均
            signal = np.mean(self.u, axis=1)
        else:
            signal = self.u[:, node_idx]
        
        # FFT
        fft_result = fft(signal)
        freqs = fftfreq(len(signal), self.dt)
        
        # 只取正频率
        positive_freqs = freqs[:len(freqs)//2]
        positive_fft = np.abs(fft_result[:len(fft_result)//2])
        
        return positive_freqs, positive_fft
    
    def compute_envelope(self, node_idx=0):
        """计算包络"""
        from scipy.signal import hilbert
        signal = self.u[:, node_idx]
        analytic_signal = hilbert(signal)
        envelope = np.abs(analytic_signal)
        return envelope
    
    def compute_group_delay(self, input_idx, output_idx):
        """计算群延迟"""
        # 计算两个位置的频谱
        freqs1, fft1 = self.compute_spectrum(input_idx)
        freqs2, fft2 = self.compute_spectrum(output_idx)
        
        # 相位差
        phase1 = np.angle(fft1)
        phase2 = np.angle(fft2)
        phase_diff = np.unwrap(phase2 - phase1)
        
        # 群延迟
        group_delay = -np.gradient(phase_diff, 2 * np.pi * freqs1)
        
        return freqs1, group_delay

8. 应用案例

8.1 案例1：二维波传播

问题描述：
在 $\times 2$ m的二维区域中，中心点激发高斯脉冲，观察波的传播过程。

物理参数：

波速： $c = 1$ m/s
网格： $100 \times 100$
脉冲宽度： $σ=0.05\sigma = 0.05$ s
仿真时间： $T = 2$ s

数值结果：

波以圆形波前向外传播，振幅随距离衰减。中心点的时间历程显示脉冲先向外传播，随后边界反射波返回。

Python实现：

def demo_wave_propagation():
    """演示波传播"""
    # 创建网格
    nx, ny = 100, 100
    Lx, Ly = 2.0, 2.0
    
    solver = WaveEquationSolver(nx, ny, Lx, Ly)
    
    # 波速和时间参数
    c = 1.0
    dx = Lx / (nx - 1)
    dt = 0.5 * dx / c  # CFL条件
    t_final = 2.0
    n_steps = int(t_final / dt)
    
    # 源函数 (高斯脉冲)
    t0 = 0.2
    sigma = 0.05
    source_func = lambda t: PulsedExcitation.gaussian_pulse(t, t0, sigma)
    source_pos = (nx // 2, ny // 2)
    
    # 求解
    u_history = solver.solve_fdtd_style(c, dt, n_steps, source_pos, source_func)
    
    return solver, u_history

结果分析：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

从结果可以看出：

波以同心圆形式向外传播
波前清晰，数值色散较小
中心点信号呈现先正后负的振荡特征

8.2 案例2：色散材料特性

问题描述：
比较Debye、Lorentz和Drude三种色散模型在 $1$ MHz到 $1$ THz频段的介电常数变化。

模型参数：

Debye模型(水)：

$ε∞=2\varepsilon_\infty = 2$
$εs=80\varepsilon_s = 80$
$τ=1\tau = 1$ ns

Lorentz模型(共振材料)：

$ε∞=2\varepsilon_\infty = 2$
$ωp=2π×10\omega_p = 2\pi \times 10$ GHz
$ω0=2π×5\omega_0 = 2\pi \times 5$ GHz
$γ=1\gamma = 1$ GHz

Drude模型(金)：

$ωp=2π×2\omega_p = 2\pi \times 2$ PHz
$γ=10\gamma = 10$ THz

结果分析：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

Debye模型：

介电常数从静态值 $80$ 逐渐下降到高频值 $2$
在 $1/τ1/\tau$ 附近有明显变化
虚部表示损耗，在过渡区最大

Lorentz模型：

在共振频率 $ω0\omega_0$ 附近出现反常色散
实部在共振区急剧下降
虚部在共振峰处最大

Drude模型：

低频时介电常数为负(金属特性)
在等离子体频率处通过零点
高频时趋近于1

8.3 案例3：脉冲在色散介质中的传播

问题描述：
比较高斯脉冲、调制高斯脉冲和Ricker子波在一维介质中的传播特性。

观察现象：

脉冲展宽：色散导致不同频率成分以不同速度传播
波形畸变：高频和低频成分分离
能量衰减：材料损耗导致振幅下降

结果分析：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

三种脉冲的传播特性：

高斯脉冲：保持形状较好，略有展宽
调制高斯：载波频率成分受色散影响，包络展宽
Ricker子波：高频成分传播更快，波形前缘变陡

8.4 案例4：瞬态信号分析

问题描述：
对含噪声的瞬态信号进行频谱分析、包络检测和群延迟计算。

信号特征：

中心频率： $5$ Hz
调制带宽： $2$ Hz
信噪比： $20$ dB

分析方法：

FFT频谱：识别信号频率成分
希尔伯特包络：提取信号幅度调制
群延迟：评估不同频率的传输时间

结果分析：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

从分析结果可以：

频谱峰值位于 $5$ Hz，与设定一致
包络显示信号幅度变化
群延迟在通带内相对平坦

8.5 案例5：矩形腔体谐振器

问题描述：
分析 $\times 0.6 \times 0.4$ m矩形腔体的谐振模式。

理论分析：

矩形腔体的谐振频率：

$f_{mnp} = \frac{c}{2} \sqrt{\left(\frac{m}{a}\right)^2 + \left(\frac{n}{b}\right)^2 + \left(\frac{p}{d}\right)^2}$

其中 $a, b, d$ 为腔体尺寸， $m, n, p$ 为模式指数。

前几个模式：

$TE_{101}$ : $f = 180.3$ MHz
$TE_{011}$ : $f = 207.9$ MHz
$TE_{102}$ : $f = 249.6$ MHz

结果分析：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

频谱显示清晰的谐振峰，对应不同模式。模式场分布显示：

$TE_{101}$ ：沿长度方向一个半波变化
$TE_{011}$ ：沿高度方向一个半波变化
高阶模式有更多场变化

8.6 案例6：时域波导传输

问题描述：
分析波导中脉冲信号的传输和色散效应。

波导参数：

宽度： $0.1$ m
截止频率： $f_c = c/(2a) = 1.5$ GHz
工作频率： $2$ GHz

群速度：

$v_g = c \sqrt{1 - \left(\frac{f_c}{f}\right)^2}$

结果分析：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

观察结果：

脉冲在波导中传播时展宽
群速度小于光速
不同频率成分以不同速度传播

8.7 案例7：TD-FEM与频域FEM比较

问题描述：
比较时域和频域有限元法的特征模式。

比较内容：

特征频率：两种方法应给出相同结果
模式形状：验证空间分布一致性
计算效率：TD-FEM一次获得多个模式

结果分析：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

前6个特征模式显示：

模式1,2：基频模式，场分布简单
模式3,4：高阶模式，有更多节点
特征频率与理论值吻合

8.8 案例8：PML边界条件

问题描述：
验证PML边界条件的吸收性能。

测试设置：

计算域： $\times 3$ m
PML厚度：20格点
源：中心高斯脉冲

性能评估：

反射系数： $10^{-4}$
波在PML中指数衰减
边界处无明显反射

结果分析：

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

PML性能验证：

波传播到边界后几乎无反射
PML区域内场强快速衰减
电导率剖面呈三次多项式分布

9. 方法比较与选择

9.1 TD-FEM vs FDTD vs 频域FEM

特性	FDTD	TD-FEM	频域FEM
网格类型	矩形/六面体	非结构化	非结构化
几何适应性	差	好	好
数值色散	各向异性	各向同性	各向同性
内存需求	低	中	中
计算效率	高	中	高
宽带分析	一次仿真	一次仿真	多次仿真
非线性材料	支持	支持	不支持
色散材料	支持	支持	需迭代

9.2 方法选择指南

选择FDTD当：

几何相对简单
需要最高计算效率
大规模问题
均匀或分层介质

选择TD-FEM当：

复杂几何结构
需要非结构化网格
色散材料分析
宽带响应需求
精细结构建模

选择频域FEM当：

单频或窄带问题
需要最高精度
谐振问题
已有成熟的频域求解器

9.3 计算复杂度分析

时间复杂度：

FDTD: $O (N)$ ，每步只需更新场值
TD-FEM: $O(N^{1.5})$ ，隐式方法需解方程组
频域FEM: $O(N^{1.5})$ ，稀疏矩阵求解

空间复杂度：

FDTD: $O (N)$ ，存储场值
TD-FEM: $O (N)$ ，存储稀疏矩阵
频域FEM: $O (N)$ ，存储稀疏矩阵

宽带分析总成本：

FDTD/TD-FEM: $\times N_t)$ ， $N_t$ 为时间步数
频域FEM: $O(N1.5×Nf)O(N^{1.5} \times N_f)$ ， $N_f$ 为频点数

对于宽带问题( $Nf≫1N_f \gg 1$ )，时域方法通常更高效。

# -*- coding: utf-8 -*-
"""
主题079: 时域有限元法 (TD-FEM)
Time-Domain Finite Element Method for Electromagnetic Simulation
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import sparse
from scipy.sparse import linalg as splinalg
from scipy.fft import fft, fftfreq
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


class TDFEM2D:
    """二维时域有限元求解器"""
    
    def __init__(self, nodes, elements, epsilon_r=1.0, mu_r=1.0, sigma=0.0):
        """
        参数:
            nodes: 节点坐标数组 (N, 2)
            elements: 三角形单元连接数组 (M, 3)
            epsilon_r: 相对介电常数
            mu_r: 相对磁导率
            sigma: 电导率
        """
        self.nodes = nodes
        self.elements = elements
        self.n_nodes = len(nodes)
        self.n_elements = len(elements)
        
        # 材料参数
        self.epsilon0 = 8.854e-12
        self.mu0 = 4 * np.pi * 1e-7
        self.epsilon = epsilon_r * self.epsilon0
        self.mu = mu_r * self.mu0
        self.sigma = sigma
        
        # 光速
        self.c = 1 / np.sqrt(self.epsilon * self.mu)
        
    def compute_element_matrices(self, element_idx):
        """计算单元矩阵"""
        elem = self.elements[element_idx]
        x = self.nodes[elem, 0]
        y = self.nodes[elem, 1]
        
        # 计算单元面积
        area = 0.5 * abs((x[1] - x[0]) * (y[2] - y[0]) - 
                         (x[2] - x[0]) * (y[1] - y[0]))
        
        # 形状函数梯度
        b = np.array([y[1] - y[2], y[2] - y[0], y[0] - y[1]])
        c = np.array([x[2] - x[1], x[0] - x[2], x[1] - x[0]])
        
        # 刚度矩阵 (1/mu * 梯度项)
        K_elem = np.outer(b, b) + np.outer(c, c)
        K_elem = K_elem / (4 * area * self.mu)
        
        # 质量矩阵 (epsilon * 势能项)
        M_elem = self.epsilon * area / 12 * np.array([[2, 1, 1],
                                                       [1, 2, 1],
                                                       [1, 1, 2]])
        
        # 阻尼矩阵 (sigma * 耗散项)
        C_elem = self.sigma * area / 12 * np.array([[2, 1, 1],
                                                     [1, 2, 1],
                                                     [1, 1, 2]])
        
        return K_elem, M_elem, C_elem, area
    
    def assemble_system(self):
        """组装全局系统矩阵"""
        rows, cols = [], []
        data_K, data_M, data_C = [], [], []
        
        for e in range(self.n_elements):
            K_elem, M_elem, C_elem, _ = self.compute_element_matrices(e)
            elem = self.elements[e]
            
            for i in range(3):
                for j in range(3):
                    rows.append(elem[i])
                    cols.append(elem[j])
                    data_K.append(K_elem[i, j])
                    data_M.append(M_elem[i, j])
                    data_C.append(C_elem[i, j])
        
        self.K = sparse.coo_matrix((data_K, (rows, cols)), 
                                    shape=(self.n_nodes, self.n_nodes)).tocsr()
        self.M = sparse.coo_matrix((data_M, (rows, cols)), 
                                    shape=(self.n_nodes, self.n_nodes)).tocsr()
        self.C = sparse.coo_matrix((data_C, (rows, cols)), 
                                    shape=(self.n_nodes, self.n_nodes)).tocsr()
        
        return self.K, self.M, self.C
    
    def newmark_timestep(self, u_n, v_n, a_n, dt, beta=0.25, gamma=0.5):
        """
        Newmark-beta时间积分
        
        参数:
            u_n: 当前位移
            v_n: 当前速度
            a_n: 当前加速度
            dt: 时间步长
            beta, gamma: Newmark参数
        """
        # 预测步
        u_pred = u_n + dt * v_n + (dt**2 / 2) * (1 - 2*beta) * a_n
        v_pred = v_n + dt * (1 - gamma) * a_n
        
        # 求解加速度
        A = self.M + gamma * dt * self.C + beta * dt**2 * self.K
        b = -self.K @ u_pred - self.C @ v_pred
        
        a_new = splinalg.spsolve(A, b)
        
        # 校正步
        u_new = u_pred + beta * dt**2 * a_new
        v_new = v_pred + gamma * dt * a_new
        
        return u_new, v_new, a_new
    
    def solve_time_domain(self, t_final, dt, source_func, source_nodes):
        """
        时域求解
        
        参数:
            t_final: 终止时间
            dt: 时间步长
            source_func: 源函数 lambda t: value
            source_nodes: 源节点索引
        """
        n_steps = int(t_final / dt)
        
        # 初始化
        u = np.zeros(self.n_nodes)
        v = np.zeros(self.n_nodes)
        a = np.zeros(self.n_nodes)
        
        # 存储结果
        u_history = np.zeros((n_steps, self.n_nodes))
        t_history = np.linspace(0, t_final, n_steps)
        
        for n in range(n_steps):
            t = n * dt
            
            # 施加源
            for node in source_nodes:
                u[node] = source_func(t)
            
            # 时间步进
            u, v, a = self.newmark_timestep(u, v, a, dt)
            
            # 存储
            u_history[n, :] = u
        
        return t_history, u_history


class WaveEquationSolver:
    """波动方程求解器"""
    
    def __init__(self, nx=50, ny=50, Lx=1.0, Ly=1.0):
        """
        参数:
            nx, ny: 网格点数
            Lx, Ly: 区域尺寸
        """
        self.nx = nx
        self.ny = ny
        self.Lx = Lx
        self.Ly = Ly
        
        # 创建网格
        self.x = np.linspace(0, Lx, nx)
        self.y = np.linspace(0, Ly, ny)
        self.dx = Lx / (nx - 1)
        self.dy = Ly / (ny - 1)
        
        self.X, self.Y = np.meshgrid(self.x, self.y)
        
    def solve_fdtd_style(self, c, dt, n_steps, source_pos, source_func):
        """
        FDTD风格的波动方程求解
        
        参数:
            c: 波速
            dt: 时间步长
            n_steps: 时间步数
            source_pos: 源位置 (ix, iy)
            source_func: 源函数 lambda n: value
        """
        # CFL条件检查
        cfl = c * dt / self.dx
        if cfl > 1.0 / np.sqrt(2):
            print(f"警告: CFL条件不满足 ({cfl:.3f} > {1/np.sqrt(2):.3f})")
        
        # 初始化场
        u_prev = np.zeros((self.ny, self.nx))
        u_curr = np.zeros((self.ny, self.nx))
        u_next = np.zeros((self.ny, self.nx))
        
        # 存储历史
        u_history = np.zeros((n_steps, self.ny, self.nx))
        
        for n in range(n_steps):
            # 有限差分更新
            for i in range(1, self.ny - 1):
                for j in range(1, self.nx - 1):
                    u_next[i, j] = (2 * u_curr[i, j] - u_prev[i, j] + 
                                    (c * dt)**2 * (
                                        (u_curr[i+1, j] - 2*u_curr[i, j] + u_curr[i-1, j]) / self.dy**2 +
                                        (u_curr[i, j+1] - 2*u_curr[i, j] + u_curr[i, j-1]) / self.dx**2
                                    ))
            
            # 施加源
            ix, iy = source_pos
            u_next[iy, ix] += source_func(n * dt)
            
            # 更新
            u_prev = u_curr.copy()
            u_curr = u_next.copy()
            
            u_history[n, :, :] = u_curr
        
        return u_history


class DispersiveMaterial:
    """色散材料模型"""
    
    def __init__(self, model_type='debye'):
        """
        参数:
            model_type: 'debye', 'lorentz', 'drude'
        """
        self.model_type = model_type
        
    def debye_model(self, omega, epsilon_inf, epsilon_s, tau):
        """
        Debye色散模型
        
        参数:
            omega: 角频率
            epsilon_inf: 高频介电常数
            epsilon_s: 静态介电常数
            tau: 弛豫时间
        """
        epsilon = epsilon_inf + (epsilon_s - epsilon_inf) / (1 + 1j * omega * tau)
        return epsilon
    
    def lorentz_model(self, omega, epsilon_inf, omega_p, gamma, omega_0):
        """
        Lorentz色散模型
        
        参数:
            omega: 角频率
            epsilon_inf: 高频介电常数
            omega_p: 等离子体频率
            gamma: 阻尼系数
            omega_0: 共振频率
        """
        epsilon = epsilon_inf + omega_p**2 / (omega_0**2 - omega**2 - 1j * gamma * omega)
        return epsilon
    
    def drude_model(self, omega, omega_p, gamma):
        """
        Drude色散模型 (金属)
        
        参数:
            omega: 角频率
            omega_p: 等离子体频率
            gamma: 碰撞频率
        """
        epsilon = 1 - omega_p**2 / (omega**2 + 1j * gamma * omega)
        return epsilon
    
    def plot_dispersion(self, f_min=1e6, f_max=1e12, n_points=1000):
        """绘制色散曲线"""
        f = np.logspace(np.log10(f_min), np.log10(f_max), n_points)
        omega = 2 * np.pi * f
        
        fig, axes = plt.subplots(1, 3, figsize=(15, 4))
        
        # Debye模型
        if self.model_type == 'debye' or self.model_type == 'all':
            epsilon_debye = self.debye_model(omega, 2.0, 80.0, 1e-9)
            axes[0].semilogx(f, np.real(epsilon_debye), 'b-', label='Real')
            axes[0].semilogx(f, np.imag(epsilon_debye), 'r--', label='Imag')
            axes[0].set_xlabel('Frequency (Hz)')
            axes[0].set_ylabel('Permittivity')
            axes[0].set_title('Debye Model')
            axes[0].legend()
            axes[0].grid(True)
        
        # Lorentz模型
        if self.model_type == 'lorentz' or self.model_type == 'all':
            epsilon_lorentz = self.lorentz_model(omega, 2.0, 2*np.pi*1e10, 1e9, 2*np.pi*5e9)
            axes[1].semilogx(f, np.real(epsilon_lorentz), 'b-', label='Real')
            axes[1].semilogx(f, np.imag(epsilon_lorentz), 'r--', label='Imag')
            axes[1].set_xlabel('Frequency (Hz)')
            axes[1].set_ylabel('Permittivity')
            axes[1].set_title('Lorentz Model')
            axes[1].legend()
            axes[1].grid(True)
        
        # Drude模型
        if self.model_type == 'drude' or self.model_type == 'all':
            epsilon_drude = self.drude_model(omega, 2*np.pi*2e15, 1e13)
            axes[2].semilogx(f, np.real(epsilon_drude), 'b-', label='Real')
            axes[2].semilogx(f, np.imag(epsilon_drude), 'r--', label='Imag')
            axes[2].set_xlabel('Frequency (Hz)')
            axes[2].set_ylabel('Permittivity')
            axes[2].set_title('Drude Model (Metal)')
            axes[2].legend()
            axes[2].grid(True)
        
        plt.tight_layout()
        return fig


class TransientAnalysis:
    """瞬态分析工具"""
    
    def __init__(self, t_history, u_history):
        """
        参数:
            t_history: 时间历史
            u_history: 场量历史
        """
        self.t = t_history
        self.u = u_history
        self.dt = t_history[1] - t_history[0] if len(t_history) > 1 else 1
        
    def compute_spectrum(self, node_idx=None):
        """计算频谱"""
        if node_idx is None:
            # 对所有节点平均
            signal = np.mean(self.u, axis=1)
        else:
            signal = self.u[:, node_idx]
        
        # FFT
        fft_result = fft(signal)
        freqs = fftfreq(len(signal), self.dt)
        
        # 只取正频率
        positive_freqs = freqs[:len(freqs)//2]
        positive_fft = np.abs(fft_result[:len(fft_result)//2])
        
        return positive_freqs, positive_fft
    
    def compute_envelope(self, node_idx=0):
        """计算包络"""
        from scipy.signal import hilbert
        signal = self.u[:, node_idx]
        analytic_signal = hilbert(signal)
        envelope = np.abs(analytic_signal)
        return envelope
    
    def compute_group_delay(self, input_idx, output_idx):
        """计算群延迟"""
        # 计算两个位置的频谱
        freqs1, fft1 = self.compute_spectrum(input_idx)
        freqs2, fft2 = self.compute_spectrum(output_idx)
        
        # 相位差
        phase1 = np.angle(fft1)
        phase2 = np.angle(fft2)
        phase_diff = np.unwrap(phase2 - phase1)
        
        # 群延迟
        group_delay = -np.gradient(phase_diff, 2 * np.pi * freqs1)
        
        return freqs1, group_delay


class PulsedExcitation:
    """脉冲激励源"""
    
    @staticmethod
    def gaussian_pulse(t, t0, sigma):
        """高斯脉冲"""
        return np.exp(-((t - t0) / sigma)**2)
    
    @staticmethod
    def modulated_gaussian(t, t0, sigma, f0):
        """调制高斯脉冲"""
        return np.exp(-((t - t0) / sigma)**2) * np.sin(2 * np.pi * f0 * t)
    
    @staticmethod
    def derivative_gaussian(t, t0, sigma):
        """高斯脉冲的导数"""
        return -2 * (t - t0) / sigma**2 * np.exp(-((t - t0) / sigma)**2)
    
    @staticmethod
    def ricker_wavelet(t, f0):
        """Ricker子波"""
        t0 = 1.0 / f0
        r = np.pi * f0 * (t - t0)
        return (1 - 2 * r**2) * np.exp(-r**2)
    
    @staticmethod
    def sine_wave(t, f0, t_rise=0):
        """正弦波（带上升沿）"""
        if t < t_rise:
            return np.sin(2 * np.pi * f0 * t) * (t / t_rise)
        else:
            return np.sin(2 * np.pi * f0 * t)


class AbsorbingBoundary:
    """吸收边界条件"""
    
    def __init__(self, nx, ny, dx, dy, dt, c, thickness=10):
        """
        参数:
            nx, ny: 网格尺寸
            dx, dy: 网格间距
            dt: 时间步长
            c: 波速
            thickness: PML厚度
        """
        self.nx = nx
        self.ny = ny
        self.dx = dx
        self.dy = dy
        self.dt = dt
        self.c = c
        self.thickness = thickness
        
        # 创建PML剖面
        self.sigma_x = self._create_pml_profile(nx, thickness, dx)
        self.sigma_y = self._create_pml_profile(ny, thickness, dy)
        
    def _create_pml_profile(self, n, thickness, d):
        """创建PML电导率剖面"""
        sigma = np.zeros(n)
        
        # 左侧PML
        for i in range(thickness):
            sigma[i] = ((thickness - i) / thickness)**3
        
        # 右侧PML
        for i in range(thickness):
            sigma[n - 1 - i] = ((thickness - i) / thickness)**3
        
        return sigma
    
    def apply_pml(self, u, u_prev):
        """应用PML边界条件"""
        u_pml = u.copy()
        
        # 在PML区域添加衰减
        for i in range(self.ny):
            for j in range(self.nx):
                sigma = self.sigma_x[j] + self.sigma_y[i]
                if sigma > 0:
                    u_pml[i, j] *= np.exp(-sigma * self.dt)
        
        return u_pml


def demo_wave_propagation():
    """演示波传播"""
    print("=" * 60)
    print("演示1: 二维波传播")
    print("=" * 60)
    
    # 创建网格
    nx, ny = 100, 100
    Lx, Ly = 2.0, 2.0
    
    solver = WaveEquationSolver(nx, ny, Lx, Ly)
    
    # 波速
    c = 1.0
    
    # 时间参数
    dx = Lx / (nx - 1)
    dt = 0.5 * dx / c  # CFL条件
    t_final = 2.0
    n_steps = int(t_final / dt)
    
    print(f"网格: {nx} x {ny}")
    print(f"空间步长: dx = {dx:.4f}")
    print(f"时间步长: dt = {dt:.6f}")
    print(f"CFL数: {c * dt / dx:.4f}")
    print(f"总步数: {n_steps}")
    
    # 源函数 (高斯脉冲)
    t0 = 0.2
    sigma = 0.05
    source_func = lambda t: PulsedExcitation.gaussian_pulse(t, t0, sigma)
    
    # 源位置 (中心)
    source_pos = (nx // 2, ny // 2)
    
    # 求解
    u_history = solver.solve_fdtd_style(c, dt, n_steps, source_pos, source_func)
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    times = [0, n_steps//4, n_steps//2, 3*n_steps//4, n_steps-1]
    titles = ['t=0', 't=T/4', 't=T/2', 't=3T/4', 't=T']
    
    for idx, (t_idx, title) in enumerate(zip(times, titles)):
        ax = axes[idx // 3, idx % 3]
        im = ax.contourf(solver.X, solver.Y, u_history[t_idx], 
                         levels=20, cmap='RdBu_r', vmin=-0.5, vmax=0.5)
        ax.set_aspect('equal')
        ax.set_title(f'Wave Propagation ({title})')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        plt.colorbar(im, ax=ax)
    
    # 中心点的时间历程
    ax = axes[1, 2]
    center_signal = u_history[:, ny//2, nx//2]
    t = np.linspace(0, t_final, n_steps)
    ax.plot(t, center_signal, 'b-', lw=1)
    ax.set_xlabel('Time')
    ax.set_ylabel('Amplitude')
    ax.set_title('Signal at Center')
    ax.grid(True)
    
    plt.tight_layout()
    plt.savefig('wave_propagation.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 wave_propagation.png")
    plt.close()
    
    return solver, u_history


def demo_dispersive_materials():
    """演示色散材料"""
    print("\n" + "=" * 60)
    print("演示2: 色散材料模型")
    print("=" * 60)
    
    material = DispersiveMaterial('all')
    
    fig = material.plot_dispersion(f_min=1e6, f_max=1e12, n_points=1000)
    
    plt.savefig('dispersive_materials.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 dispersive_materials.png")
    plt.close()
    
    # 打印典型值
    f_test = 1e9
    omega_test = 2 * np.pi * f_test
    
    print(f"\n在 f = {f_test/1e9:.1f} GHz 处的介电常数:")
    
    epsilon_debye = material.debye_model(omega_test, 2.0, 80.0, 1e-9)
    print(f"  Debye模型: {epsilon_debye:.2f}")
    
    epsilon_lorentz = material.lorentz_model(omega_test, 2.0, 2*np.pi*1e10, 1e9, 2*np.pi*5e9)
    print(f"  Lorentz模型: {epsilon_lorentz:.2f}")
    
    epsilon_drude = material.drude_model(omega_test, 2*np.pi*2e15, 1e13)
    print(f"  Drude模型: {epsilon_drude:.2f}")
    
    return material


def demo_pulse_propagation():
    """演示脉冲传播"""
    print("\n" + "=" * 60)
    print("演示3: 脉冲在色散介质中的传播")
    print("=" * 60)
    
    # 创建一维仿真
    nx = 500
    L = 10.0
    dx = L / nx
    x = np.linspace(0, L, nx)
    
    # 时间参数
    c = 1.0
    dt = 0.8 * dx / c
    t_final = 8.0
    n_steps = int(t_final / dt)
    
    # 初始化
    u_prev = np.zeros(nx)
    u_curr = np.zeros(nx)
    u_next = np.zeros(nx)
    
    # 存储历史
    u_history = np.zeros((n_steps, nx))
    
    # 脉冲参数
    t0 = 1.0
    sigma = 0.2
    f0 = 2.0
    
    # 选择脉冲类型
    pulse_types = ['gaussian', 'modulated', 'ricker']
    pulse_labels = ['Gaussian', 'Modulated Gaussian', 'Ricker']
    
    fig, axes = plt.subplots(3, 3, figsize=(15, 12))
    
    for p_idx, (pulse_type, label) in enumerate(zip(pulse_types, pulse_labels)):
        u_prev = np.zeros(nx)
        u_curr = np.zeros(nx)
        
        for n in range(n_steps):
            t = n * dt
            
            # 有限差分更新
            for i in range(1, nx - 1):
                u_next[i] = (2 * u_curr[i] - u_prev[i] + 
                            (c * dt / dx)**2 * (u_curr[i+1] - 2*u_curr[i] + u_curr[i-1]))
            
            # 施加源
            if pulse_type == 'gaussian':
                u_next[10] = PulsedExcitation.gaussian_pulse(t, t0, sigma)
            elif pulse_type == 'modulated':
                u_next[10] = PulsedExcitation.modulated_gaussian(t, t0, sigma, f0)
            elif pulse_type == 'ricker':
                u_next[10] = PulsedExcitation.ricker_wavelet(t, f0)
            
            # 更新
            u_prev = u_curr.copy()
            u_curr = u_next.copy()
            
            # 存储特定时刻
            if n in [int(n_steps*0.1), int(n_steps*0.5), int(n_steps*0.9)]:
                snap_idx = [int(n_steps*0.1), int(n_steps*0.5), int(n_steps*0.9)].index(n)
                axes[p_idx, snap_idx].plot(x, u_curr, 'b-', lw=1)
                axes[p_idx, snap_idx].set_xlabel('x')
                axes[p_idx, snap_idx].set_ylabel('Amplitude')
                axes[p_idx, snap_idx].set_title(f'{label} - t={t:.2f}')
                axes[p_idx, snap_idx].set_xlim([0, L])
                axes[p_idx, snap_idx].grid(True)
    
    plt.tight_layout()
    plt.savefig('pulse_propagation.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 pulse_propagation.png")
    plt.close()
    
    return x, u_history


def demo_transient_analysis():
    """演示瞬态分析"""
    print("\n" + "=" * 60)
    print("演示4: 瞬态信号分析")
    print("=" * 60)
    
    # 生成测试信号
    fs = 1000  # 采样频率
    t = np.linspace(0, 1, fs)
    
    # 多频率信号
    f1, f2, f3 = 50, 120, 200
    signal = (np.sin(2 * np.pi * f1 * t) + 
              0.5 * np.sin(2 * np.pi * f2 * t) + 
              0.3 * np.sin(2 * np.pi * f3 * t))
    
    # 添加噪声
    signal += 0.1 * np.random.randn(len(t))
    
    # 创建瞬态分析对象
    u_history = signal.reshape(-1, 1)
    analyzer = TransientAnalysis(t, u_history)
    
    # 计算频谱
    freqs, spectrum = analyzer.compute_spectrum(0)
    
    # 计算包络
    envelope = analyzer.compute_envelope(0)
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 时域信号
    ax = axes[0, 0]
    ax.plot(t, signal, 'b-', lw=0.5, label='Signal')
    ax.plot(t, envelope, 'r--', lw=2, label='Envelope')
    ax.plot(t, -envelope, 'r--', lw=2)
    ax.set_xlabel('Time (s)')
    ax.set_ylabel('Amplitude')
    ax.set_title('Time Domain Signal')
    ax.legend()
    ax.grid(True)
    
    # 频谱
    ax = axes[0, 1]
    ax.plot(freqs, spectrum, 'b-', lw=1)
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Magnitude')
    ax.set_title('Frequency Spectrum')
    ax.set_xlim([0, 300])
    ax.grid(True)
    
    # 频谱（对数坐标）
    ax = axes[1, 0]
    ax.semilogy(freqs, spectrum + 1e-10, 'b-', lw=1)
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Magnitude (log)')
    ax.set_title('Frequency Spectrum (Log Scale)')
    ax.set_xlim([0, 500])
    ax.grid(True)
    
    # 相位谱
    fft_result = fft(signal)
    phase = np.angle(fft_result[:len(fft_result)//2])
    ax = axes[1, 1]
    ax.plot(freqs, np.unwrap(phase), 'b-', lw=1)
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Phase (rad)')
    ax.set_title('Phase Spectrum')
    ax.set_xlim([0, 300])
    ax.grid(True)
    
    plt.tight_layout()
    plt.savefig('transient_analysis.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 transient_analysis.png")
    plt.close()
    
    # 打印检测到的频率
    peaks_idx = np.argsort(spectrum)[-5:]
    print(f"\n主要频率成分:")
    for idx in sorted(peaks_idx, key=lambda i: spectrum[i], reverse=True)[:3]:
        print(f"  {freqs[idx]:.1f} Hz (幅度: {spectrum[idx]:.2f})")
    
    return t, signal, freqs, spectrum


def demo_cavity_resonator():
    """演示腔体谐振器"""
    print("\n" + "=" * 60)
    print("演示5: 矩形腔体谐振器")
    print("=" * 60)
    
    # 腔体尺寸
    Lx, Ly = 1.0, 0.5
    
    # 创建网格
    nx, ny = 80, 40
    dx, dy = Lx / (nx - 1), Ly / (ny - 1)
    
    x = np.linspace(0, Lx, nx)
    y = np.linspace(0, Ly, ny)
    X, Y = np.meshgrid(x, y)
    
    # 波速
    c = 1.0
    
    # 时间参数
    dt = 0.5 * min(dx, dy) / c
    t_final = 10.0
    n_steps = int(t_final / dt)
    
    # 初始化
    u_prev = np.zeros((ny, nx))
    u_curr = np.zeros((ny, nx))
    u_next = np.zeros((ny, nx))
    
    # 初始条件 (高斯分布)
    x0, y0 = Lx / 3, Ly / 2
    sigma = 0.05
    u_curr = np.exp(-((X - x0)**2 + (Y - y0)**2) / (2 * sigma**2))
    
    # 存储历史
    u_history = np.zeros((n_steps, ny, nx))
    
    # 时间步进
    for n in range(n_steps):
        # 内部点更新
        for i in range(1, ny - 1):
            for j in range(1, nx - 1):
                u_next[i, j] = (2 * u_curr[i, j] - u_prev[i, j] + 
                                (c * dt)**2 * (
                                    (u_curr[i+1, j] - 2*u_curr[i, j] + u_curr[i-1, j]) / dy**2 +
                                    (u_curr[i, j+1] - 2*u_curr[i, j] + u_curr[i, j-1]) / dx**2
                                ))
        
        # 边界条件 (PEC: u=0)
        u_next[0, :] = 0
        u_next[-1, :] = 0
        u_next[:, 0] = 0
        u_next[:, -1] = 0
        
        # 更新
        u_prev = u_curr.copy()
        u_curr = u_next.copy()
        
        u_history[n, :, :] = u_curr
    
    # 计算谐振频率 (理论值)
    m, n_mode = 1, 1
    f_theory = (c / 2) * np.sqrt((m/Lx)**2 + (n_mode/Ly)**2)
    
    print(f"理论谐振频率 (m={m}, n={n_mode}): {f_theory:.4f} Hz")
    
    # 从时域信号提取频率
    center_signal = u_history[:, ny//2, nx//2]
    fft_result = fft(center_signal)
    freqs = fftfreq(n_steps, dt)
    positive_freqs = freqs[:n_steps//2]
    positive_fft = np.abs(fft_result[:n_steps//2])
    
    # 找到峰值
    peak_idx = np.argmax(positive_fft[1:]) + 1
    f_numerical = positive_freqs[peak_idx]
    print(f"数值谐振频率: {f_numerical:.4f} Hz")
    print(f"相对误差: {abs(f_numerical - f_theory)/f_theory * 100:.2f}%")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 不同时刻的场分布
    times = [0, n_steps//4, n_steps//2]
    titles = ['Initial', 't=T/4', 't=T/2']
    
    for idx, (t_idx, title) in enumerate(zip(times, titles)):
        ax = axes[0, idx]
        im = ax.contourf(X, Y, u_history[t_idx], levels=20, cmap='RdBu_r')
        ax.set_aspect('equal')
        ax.set_title(f'Cavity Mode ({title})')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        plt.colorbar(im, ax=ax)
    
    # 中心点时间历程
    ax = axes[1, 0]
    t = np.linspace(0, t_final, n_steps)
    ax.plot(t, center_signal, 'b-', lw=1)
    ax.set_xlabel('Time')
    ax.set_ylabel('Amplitude')
    ax.set_title('Signal at Center')
    ax.grid(True)
    
    # 频谱
    ax = axes[1, 1]
    ax.plot(positive_freqs, positive_fft, 'b-', lw=1)
    ax.axvline(f_theory, color='r', linestyle='--', label='Theory')
    ax.axvline(f_numerical, color='g', linestyle=':', label='Numerical')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Magnitude')
    ax.set_title('Frequency Spectrum')
    ax.set_xlim([0, 5])
    ax.legend()
    ax.grid(True)
    
    # 模式形状
    ax = axes[1, 2]
    mode_shape = u_history[n_steps//4, :, :]
    im = ax.contourf(X, Y, mode_shape, levels=20, cmap='RdBu_r')
    ax.set_aspect('equal')
    ax.set_title(f'Mode Shape (m={m}, n={n_mode})')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    plt.colorbar(im, ax=ax)
    
    plt.tight_layout()
    plt.savefig('cavity_resonator.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 cavity_resonator.png")
    plt.close()
    
    return X, Y, u_history, f_theory, f_numerical


def demo_waveguide_td():
    """演示时域波导分析"""
    print("\n" + "=" * 60)
    print("演示6: 时域波导传输")
    print("=" * 60)
    
    # 波导尺寸
    Lx, Ly = 4.0, 0.5
    
    # 创建网格
    nx, ny = 200, 25
    dx, dy = Lx / (nx - 1), Ly / (ny - 1)
    
    x = np.linspace(0, Lx, nx)
    y = np.linspace(0, Ly, ny)
    X, Y = np.meshgrid(x, y)
    
    # 波速
    c = 1.0
    
    # 时间参数
    dt = 0.5 * min(dx, dy) / c
    t_final = 6.0
    n_steps = int(t_final / dt)
    
    # 初始化
    u_prev = np.zeros((ny, nx))
    u_curr = np.zeros((ny, nx))
    u_next = np.zeros((ny, nx))
    
    # 存储历史
    u_history = np.zeros((n_steps, ny, nx))
    
    # 波导截止频率
    a = Ly
    fc = c / (2 * a)
    print(f"波导宽度: {a:.3f} m")
    print(f"截止频率: {fc:.4f} Hz")
    
    # 激励频率 (高于截止)
    f_excitation = 1.5 * fc
    print(f"激励频率: {f_excitation:.4f} Hz")
    
    # 时间步进
    for n in range(n_steps):
        t = n * dt
        
        # 内部点更新
        for i in range(1, ny - 1):
            for j in range(1, nx - 1):
                u_next[i, j] = (2 * u_curr[i, j] - u_prev[i, j] + 
                                (c * dt)**2 * (
                                    (u_curr[i+1, j] - 2*u_curr[i, j] + u_curr[i-1, j]) / dy**2 +
                                    (u_curr[i, j+1] - 2*u_curr[i, j] + u_curr[i, j-1]) / dx**2
                                ))
        
        # 边界条件
        # PEC侧壁
        u_next[0, :] = 0
        u_next[-1, :] = 0
        
        # 输入端口 (正弦激励)
        u_next[:, 0] = np.sin(np.pi * y / Ly) * np.sin(2 * np.pi * f_excitation * t)
        
        # 输出端口 (一阶吸收边界)
        u_next[:, -1] = u_curr[:, -2] + (c * dt - dx) / (c * dt + dx) * (u_next[:, -2] - u_curr[:, -1])
        
        # 更新
        u_prev = u_curr.copy()
        u_curr = u_next.copy()
        
        u_history[n, :, :] = u_curr
    
    # 计算群速度
    vg = c * np.sqrt(1 - (fc / f_excitation)**2)
    print(f"群速度: {vg:.4f} m/s")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 不同时刻的场分布
    times = [n_steps//6, n_steps//3, n_steps//2]
    titles = ['t=T/6', 't=T/3', 't=T/2']
    
    for idx, (t_idx, title) in enumerate(zip(times, titles)):
        ax = axes[0, idx]
        im = ax.contourf(X, Y, u_history[t_idx], levels=20, cmap='RdBu_r')
        ax.set_aspect('equal')
        ax.set_title(f'Waveguide ({title})')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        plt.colorbar(im, ax=ax)
    
    # 输入和输出信号
    ax = axes[1, 0]
    t = np.linspace(0, t_final, n_steps)
    input_signal = u_history[:, ny//2, 5]
    output_signal = u_history[:, ny//2, -5]
    ax.plot(t, input_signal, 'b-', lw=1, label='Input')
    ax.plot(t, output_signal, 'r-', lw=1, label='Output')
    ax.set_xlabel('Time')
    ax.set_ylabel('Amplitude')
    ax.set_title('Input vs Output')
    ax.legend()
    ax.grid(True)
    
    # 传播延迟
    ax = axes[1, 1]
    correlation = np.correlate(output_signal, input_signal, mode='full')
    lag = np.argmax(correlation) - len(input_signal) + 1
    delay = lag * dt
    ax.plot(t, input_signal, 'b-', lw=1, alpha=0.5, label='Input')
    ax.plot(t, output_signal, 'r-', lw=1, alpha=0.5, label='Output')
    ax.plot(t + delay, input_signal, 'g--', lw=2, label=f'Shifted (delay={delay:.3f}s)')
    ax.set_xlabel('Time')
    ax.set_ylabel('Amplitude')
    ax.set_title('Time Delay Estimation')
    ax.legend()
    ax.grid(True)
    
    # 横截面场分布
    ax = axes[1, 2]
    for t_idx in [n_steps//3, n_steps//2, 2*n_steps//3]:
        ax.plot(y, u_history[t_idx, :, nx//2], lw=1, label=f't={t_idx*dt:.2f}')
    ax.set_xlabel('y')
    ax.set_ylabel('Amplitude')
    ax.set_title('Cross-section at x=L/2')
    ax.legend()
    ax.grid(True)
    
    plt.tight_layout()
    plt.savefig('waveguide_td.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 waveguide_td.png")
    plt.close()
    
    return X, Y, u_history, vg


def demo_td_fem_comparison():
    """演示TD-FEM与FD-FEM比较"""
    print("\n" + "=" * 60)
    print("演示7: 时域与频域FEM比较")
    print("=" * 60)
    
    # 创建简单网格
    n = 30
    x = np.linspace(0, 1, n)
    y = np.linspace(0, 1, n)
    X, Y = np.meshgrid(x, y)
    
    nodes = np.column_stack([X.flatten(), Y.flatten()])
    
    # 创建三角形单元
    elements = []
    for i in range(n-1):
        for j in range(n-1):
            n1 = i * n + j
            n2 = i * n + (j + 1)
            n3 = (i + 1) * n + j
            n4 = (i + 1) * n + (j + 1)
            
            elements.append([n1, n2, n3])
            elements.append([n2, n4, n3])
    
    elements = np.array(elements)
    
    # 创建TD-FEM求解器
    td_fem = TDFEM2D(nodes, elements, epsilon_r=1.0, mu_r=1.0)
    K, M, C = td_fem.assemble_system()
    
    print(f"网格节点数: {td_fem.n_nodes}")
    print(f"单元数: {td_fem.n_elements}")
    print(f"刚度矩阵非零元: {K.nnz}")
    print(f"质量矩阵非零元: {M.nnz}")
    
    # 计算特征频率 (频域分析)
    # 求解广义特征值问题 K phi = omega^2 M phi
    n_modes = 10
    eigenvalues, eigenvectors = splinalg.eigsh(K, k=n_modes, M=M, sigma=0, which='LM')
    
    # 特征频率
    omega = np.sqrt(np.real(eigenvalues))
    f_eigen = omega / (2 * np.pi)
    
    print(f"\n前{n_modes}个特征频率:")
    for i, f in enumerate(f_eigen):
        print(f"  模式 {i+1}: {f:.4f} Hz")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 前几个特征模式
    for idx in range(6):
        ax = axes[idx // 3, idx % 3]
        mode_shape = eigenvectors[:, idx].reshape(n, n)
        im = ax.contourf(X, Y, mode_shape, levels=20, cmap='RdBu_r')
        ax.set_aspect('equal')
        ax.set_title(f'Mode {idx+1}: f={f_eigen[idx]:.3f} Hz')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        plt.colorbar(im, ax=ax)
    
    plt.tight_layout()
    plt.savefig('td_fem_comparison.png', dpi=150, bbox_inches='tight')
    print("\n结果已保存到 td_fem_comparison.png")
    plt.close()
    
    return f_eigen, eigenvectors


def demo_pml_boundary():
    """演示PML边界条件"""
    print("\n" + "=" * 60)
    print("演示8: 完美匹配层(PML)边界")
    print("=" * 60)
    
    # 创建网格
    nx, ny = 120, 120
    Lx, Ly = 3.0, 3.0
    dx, dy = Lx / (nx - 1), Ly / (ny - 1)
    
    x = np.linspace(0, Lx, nx)
    y = np.linspace(0, Ly, ny)
    X, Y = np.meshgrid(x, y)
    
    # 波速
    c = 1.0
    
    # 时间参数
    dt = 0.5 * min(dx, dy) / c
    t_final = 3.0
    n_steps = int(t_final / dt)
    
    # PML参数
    pml_thickness = 20
    pml = AbsorbingBoundary(nx, ny, dx, dy, dt, c, pml_thickness)
    
    print(f"计算域: {Lx:.2f} x {Ly:.2f} m")
    print(f"PML厚度: {pml_thickness} 格点")
    print(f"PML物理厚度: {pml_thickness * dx:.3f} m")
    
    # 初始化
    u_prev = np.zeros((ny, nx))
    u_curr = np.zeros((ny, nx))
    u_next = np.zeros((ny, nx))
    
    # 存储历史
    u_history = np.zeros((n_steps, ny, nx))
    
    # 源位置 (中心)
    source_x, source_y = nx // 2, ny // 2
    
    # 时间步进
    for n in range(n_steps):
        t = n * dt
        
        # 内部点更新
        for i in range(1, ny - 1):
            for j in range(1, nx - 1):
                u_next[i, j] = (2 * u_curr[i, j] - u_prev[i, j] + 
                                (c * dt)**2 * (
                                    (u_curr[i+1, j] - 2*u_curr[i, j] + u_curr[i-1, j]) / dy**2 +
                                    (u_curr[i, j+1] - 2*u_curr[i, j] + u_curr[i, j-1]) / dx**2
                                ))
        
        # 施加源 (高斯脉冲)
        t0 = 0.5
        sigma = 0.1
        u_next[source_y, source_x] += np.exp(-((t - t0) / sigma)**2)
        
        # 应用PML
        u_next = pml.apply_pml(u_next, u_prev)
        
        # 更新
        u_prev = u_curr.copy()
        u_curr = u_next.copy()
        
        u_history[n, :, :] = u_curr
    
    # 计算反射系数 (简化)
    # 在PML边界处测量场强
    pml_boundary_field = np.abs(u_history[:, :, pml_thickness])
    max_field = np.max(np.abs(u_history))
    reflection_estimate = np.max(pml_boundary_field) / max_field
    
    print(f"估计反射系数: {reflection_estimate:.4e}")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 不同时刻的场分布
    times = [int(n_steps*0.2), int(n_steps*0.4), int(n_steps*0.6)]
    titles = ['t=0.2T', 't=0.4T', 't=0.6T']
    
    for idx, (t_idx, title) in enumerate(zip(times, titles)):
        ax = axes[0, idx]
        im = ax.contourf(X, Y, u_history[t_idx], levels=20, cmap='RdBu_r', 
                         vmin=-0.3, vmax=0.3)
        
        # 绘制PML边界
        pml_x = pml_thickness * dx
        ax.axvline(pml_x, color='w', linestyle='--', alpha=0.5)
        ax.axvline(Lx - pml_x, color='w', linestyle='--', alpha=0.5)
        ax.axhline(pml_x, color='w', linestyle='--', alpha=0.5)
        ax.axhline(Ly - pml_x, color='w', linestyle='--', alpha=0.5)
        
        ax.set_aspect('equal')
        ax.set_title(f'Wave with PML ({title})')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        plt.colorbar(im, ax=ax)
    
    # 中心点信号
    ax = axes[1, 0]
    t = np.linspace(0, t_final, n_steps)
    center_signal = u_history[:, ny//2, nx//2]
    ax.plot(t, center_signal, 'b-', lw=1)
    ax.set_xlabel('Time')
    ax.set_ylabel('Amplitude')
    ax.set_title('Signal at Center')
    ax.grid(True)
    
    # PML边界信号
    ax = axes[1, 1]
    pml_signal = u_history[:, ny//2, pml_thickness]
    ax.plot(t, pml_signal, 'r-', lw=1)
    ax.set_xlabel('Time')
    ax.set_ylabel('Amplitude')
    ax.set_title('Signal at PML Boundary')
    ax.grid(True)
    
    # PML电导率剖面
    ax = axes[1, 2]
    ax.plot(x[:pml_thickness], pml.sigma_x[:pml_thickness], 'b-', lw=2, label='σx')
    ax.plot(y[:pml_thickness], pml.sigma_y[:pml_thickness], 'r--', lw=2, label='σy')
    ax.set_xlabel('Position')
    ax.set_ylabel('Conductivity')
    ax.set_title('PML Conductivity Profile')
    ax.legend()
    ax.grid(True)
    
    plt.tight_layout()
    plt.savefig('pml_boundary.png', dpi=150, bbox_inches='tight')
    print("结果已保存到 pml_boundary.png")
    plt.close()
    
    return X, Y, u_history, pml


def demo_td_fem_summary():
    """TD-FEM方法总结"""
    print("\n" + "=" * 60)
    print("演示9: 时域有限元方法总结")
    print("=" * 60)
    
    # 方法特性对比
    methods = ['FDTD', 'TD-FEM', 'Frequency FEM']
    characteristics = {
        '网格类型': ['矩形/六面体', '非结构化', '非结构化'],
        '几何适应性': ['差', '好', '好'],
        '色散特性': ['各向异性', '各向同性', '各向同性'],
        '内存需求': ['低', '中', '中'],
        '计算效率': ['高', '中', '高'],
        '适用问题': ['简单几何', '复杂几何', '单频/窄带'],
        '宽带分析': ['一次仿真', '一次仿真', '多次仿真'],
        '非线性材料': ['支持', '支持', '不支持']
    }
    
    print("\n方法特性对比:")
    print("-" * 80)
    for key, values in characteristics.items():
        print(f"\n{key}:")
        for method, value in zip(methods, values):
            print(f"  {method:20s}: {value}")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 适用问题类型
    ax = axes[0, 0]
    categories = ['简单几何', '复杂几何', '色散材料', '非线性', '宽带']
    fdtd_score = [0.9, 0.4, 0.8, 0.7, 0.9]
    tdfem_score = [0.7, 0.9, 0.9, 0.8, 0.9]
    freqfem_score = [0.6, 0.9, 0.9, 0.3, 0.4]
    
    x = np.arange(len(categories))
    width = 0.25
    
    ax.bar(x - width, fdtd_score, width, label='FDTD', color='blue', alpha=0.7)
    ax.bar(x, tdfem_score, width, label='TD-FEM', color='red', alpha=0.7)
    ax.bar(x + width, freqfem_score, width, label='Freq-FEM', color='green', alpha=0.7)
    
    ax.set_ylabel('Suitability Score')
    ax.set_title('Method Suitability for Different Problems')
    ax.set_xticks(x)
    ax.set_xticklabels(categories, rotation=15)
    ax.legend()
    ax.set_ylim(0, 1.2)
    ax.grid(True, axis='y', alpha=0.3)
    
    # 计算复杂度
    ax = axes[0, 1]
    n_values = np.logspace(2, 5, 50)
    fdtd_complexity = n_values * 1e-6
    tdfem_complexity = n_values**1.5 * 1e-9
    freqfem_complexity = n_values**1.5 * 1e-9
    
    ax.loglog(n_values, fdtd_complexity, 'b-', lw=2, label='FDTD O(N)')
    ax.loglog(n_values, tdfem_complexity, 'r-', lw=2, label='TD-FEM O(N^1.5)')
    ax.loglog(n_values, freqfem_complexity, 'g-', lw=2, label='Freq-FEM O(N^1.5)')
    ax.set_xlabel('Degrees of Freedom (N)')
    ax.set_ylabel('Relative Computation Time')
    ax.set_title('Computational Complexity Comparison')
    ax.legend()
    ax.grid(True, which='both', ls='--', alpha=0.5)
    
    # 内存使用
    ax = axes[1, 0]
    fdtd_memory = n_values * 6  # 6个场分量
    tdfem_memory = n_values * 50  # 稀疏矩阵
    freqfem_memory = n_values * 50
    
    ax.loglog(n_values, fdtd_memory, 'b-', lw=2, label='FDTD')
    ax.loglog(n_values, tdfem_memory, 'r-', lw=2, label='TD-FEM')
    ax.loglog(n_values, freqfem_memory, 'g-', lw=2, label='Freq-FEM')
    ax.set_xlabel('Degrees of Freedom (N)')
    ax.set_ylabel('Memory Usage (elements)')
    ax.set_title('Memory Usage Comparison')
    ax.legend()
    ax.grid(True, which='both', ls='--', alpha=0.5)
    
    # 应用范围
    ax = axes[1, 1]
    applications = ['天线设计', '电磁散射', '微波器件', '生物电磁', '电路仿真']
    fdtd_applicability = [0.8, 0.9, 0.7, 0.6, 0.9]
    tdfem_applicability = [0.9, 0.9, 0.9, 0.9, 0.7]
    freqfem_applicability = [0.9, 0.8, 0.9, 0.8, 0.8]
    
    y_pos = np.arange(len(applications))
    
    ax.barh(y_pos - 0.25, fdtd_applicability, 0.25, label='FDTD', 
            color='blue', alpha=0.7)
    ax.barh(y_pos, tdfem_applicability, 0.25, label='TD-FEM', 
            color='red', alpha=0.7)
    ax.barh(y_pos + 0.25, freqfem_applicability, 0.25, label='Freq-FEM', 
            color='green', alpha=0.7)
    
    ax.set_xlabel('Applicability Score')
    ax.set_title('Application Domain Suitability')
    ax.set_yticks(y_pos)
    ax.set_yticklabels(applications)
    ax.legend()
    ax.set_xlim(0, 1.1)
    ax.grid(True, axis='x', alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('td_fem_summary.png', dpi=150, bbox_inches='tight')
    print("\n结果已保存到 td_fem_summary.png")
    plt.close()


def main():
    """主函数"""
    print("=" * 70)
    print("主题079: 时域有限元法 (TD-FEM)")
    print("Time-Domain Finite Element Method for Electromagnetics")
    print("=" * 70)
    
    # 运行所有演示
    demo_wave_propagation()
    demo_dispersive_materials()
    demo_pulse_propagation()
    demo_transient_analysis()
    demo_cavity_resonator()
    demo_waveguide_td()
    demo_td_fem_comparison()
    demo_pml_boundary()
    demo_td_fem_summary()
    
    print("\n" + "=" * 70)
    print("所有演示完成!")
    print("=" * 70)


if __name__ == "__main__":
    main()