这篇笔记在做什么
这一篇整理的是二维各向异性介质散射问题。原始会话里先从更一般的 Maxwell 对应关系出发，随后又收束到“无磁、仅电各向异性”的情形；这里我把两部分都保留下来，方便你同时看到一般框架和具体约束后的推导结果。
文章后半部分还保留了二维与三维在离散误差、污染误差、自由度规模和求解策略上的系统对比。
二维各向异性介质电磁散射问题的有限元方法
从 Maxwell 方程到 Helmholtz 方程 · 完整参数对应 · 二维/三维误差分析
一、从 Maxwell 方程推导二维标量方程
1.1 问题设置
二维问题：所有场量与 $z$ 无关（$\partial/\partial z = 0$）。
各向异性介质参数（对角张量）：
$$ \bar{\bar{\varepsilon}}_r = \begin{pmatrix}\varepsilon_{xx}&0&0\\0&\varepsilon_{yy}&0\\0&0&\varepsilon_{zz}\end{pmatrix}, \quad \bar{\bar{\mu}}_r = \begin{pmatrix}\mu_{xx}&0&0\\0&\mu_{yy}&0\\0&0&\mu_{zz}\end{pmatrix} $$时间约定：$e^{-j\omega t}$
Maxwell 方程：
$$ \nabla\times\mathbf{E} = j\omega\mu_0\bar{\bar{\mu}}_r\mathbf{H} $$$$ \nabla\times\mathbf{H} = -j\omega\varepsilon_0\bar{\bar{\varepsilon}}_r\mathbf{E} $$1.2 TM 极化（$E_z$ 极化）
场分量：$\mathbf{E} = E_z\hat{z}$，$\mathbf{H} = H_x\hat{x} + H_y\hat{y}$
从 Faraday 定律：
$$ \nabla\times(E_z\hat{z}) = \frac{\partial E_z}{\partial y}\hat{x} - \frac{\partial E_z}{\partial x}\hat{y} = j\omega\mu_0(\mu_{xx}H_x\hat{x}+\mu_{yy}H_y\hat{y}) $$$$ H_x = \frac{1}{j\omega\mu_0\mu_{xx}}\frac{\partial E_z}{\partial y}, \quad H_y = \frac{-1}{j\omega\mu_0\mu_{yy}}\frac{\partial E_z}{\partial x} $$从 Ampere 定律：
$$ \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} = -j\omega\varepsilon_0\varepsilon_{zz}E_z $$代入 $H_x, H_y$：
$$ \frac{\partial}{\partial x}\left(\frac{-1}{j\omega\mu_0\mu_{yy}}\frac{\partial E_z}{\partial x}\right) - \frac{\partial}{\partial y}\left(\frac{1}{j\omega\mu_0\mu_{xx}}\frac{\partial E_z}{\partial y}\right) = -j\omega\varepsilon_0\varepsilon_{zz}E_z $$两边乘以 $-j\omega\mu_0$：
$$ \frac{\partial}{\partial x}\left(\frac{1}{\mu_{yy}}\frac{\partial E_z}{\partial x}\right) + \frac{\partial}{\partial y}\left(\frac{1}{\mu_{xx}}\frac{\partial E_z}{\partial y}\right) + k_0^2\varepsilon_{zz}E_z = 0 $$写成标准 Helmholtz 形式：
$$ \boxed{ \begin{aligned} &-\frac{\partial}{\partial x}\left(\frac{1}{\mu_{yy}}\frac{\partial E_z}{\partial x}\right) \\ &\quad - \frac{\partial}{\partial y}\left(\frac{1}{\mu_{xx}}\frac{\partial E_z}{\partial y}\right) \\ &\quad + (-k_0^2\varepsilon_{zz})E_z = 0 \end{aligned} } \tag{TM} $$1.3 TE 极化（$H_z$ 极化）
场分量：$\mathbf{H} = H_z\hat{z}$，$\mathbf{E} = E_x\hat{x} + E_y\hat{y}$
完全类似的推导给出：
$$ \boxed{ \begin{aligned} &-\frac{\partial}{\partial x}\left(\frac{1}{\varepsilon_{yy}}\frac{\partial H_z}{\partial x}\right) \\ &\quad - \frac{\partial}{\partial y}\left(\frac{1}{\varepsilon_{xx}}\frac{\partial H_z}{\partial y}\right) \\ &\quad + (-k_0^2\mu_{zz})H_z = 0 \end{aligned} } \tag{TE} $$1.4 与通用 Helmholtz 方程的完整对应
通用形式：
$$ -\frac{\partial}{\partial x}\left(\alpha_x\frac{\partial\phi}{\partial x}\right) - \frac{\partial}{\partial y}\left(\alpha_y\frac{\partial\phi}{\partial y}\right) + \beta\,\phi = f $$| 参数 | TM极化（$E_z$） | TE极化（$H_z$） |
|---|---|---|
| $\phi$ | $E_z$ | $H_z$ |
| $\alpha_x$ | $1/\mu_{yy}$ | $1/\varepsilon_{yy}$ |
| $\alpha_y$ | $1/\mu_{xx}$ | $1/\varepsilon_{xx}$ |
| $\beta$ | $-k_0^2\varepsilon_{zz}$ | $-k_0^2\mu_{zz}$ |
| $f$ | 0（齐次） | 0（齐次） |
无磁介质（$\mu_r = 1$，即 $\mu_{xx}=\mu_{yy}=\mu_{zz}=1$）：
| 参数 | TM极化 | TE极化 |
|---|---|---|
| $\alpha_x$ | $1$ | $1/\varepsilon_{yy}$ |
| $\alpha_y$ | $1$ | $1/\varepsilon_{xx}$ |
| $\beta$ | $-k_0^2\varepsilon_{zz}$ | $-k_0^2$ |
关键观察：
- TM 极化：各向异性只出现在 $\beta$ 中（$\alpha_x=\alpha_y=1$），刚度矩阵各向同性
- TE 极化：各向异性出现在 $\alpha$ 中（$\alpha_x\neq\alpha_y$），刚度矩阵各向异性
二、散射场公式
2.1 总场分解
$$ E_z^{tot} = E_z^{inc} + E_z^{sca} $$入射平面波（$e^{-j\omega t}$，沿 $\hat{k} = \cos\theta_i\hat{x}+\sin\theta_i\hat{y}$ 传播）：
$$ E_z^{inc} = E_0\,e^{jk_0(x\cos\theta_i+y\sin\theta_i)} $$注意 $e^{-j\omega t}$ 下空间因子为 $e^{+j\mathbf{k}\cdot\mathbf{r}}$。
2.2 散射场方程（TM，无磁）
入射场满足自由空间方程：
$$ -\nabla^2 E_z^{inc} - k_0^2 E_z^{inc} = 0 $$总场在整个域满足：
$$ -\nabla^2 E_z^{tot} - k_0^2\varepsilon_{zz}(\mathbf{r})E_z^{tot} = 0 $$其中 $\varepsilon_{zz}(\mathbf{r}) = \varepsilon_{zz}^{(d)}$（介质内）或 $1$（自由空间）。
代入 $E_z^{tot} = E_z^{inc} + E_z^{sca}$，相减：
$$ -\nabla^2 E_z^{sca} - k_0^2\varepsilon_{zz}(\mathbf{r})\,E_z^{sca} = k_0^2\left(\varepsilon_{zz}(\mathbf{r})-1\right)E_z^{inc} $$与通用形式的对应：
$$ \boxed{ \begin{aligned} \phi &= E_z^{sca}, \quad \alpha_x = \alpha_y = 1, \\ \beta &= -k_0^2\varepsilon_{zz}(\mathbf{r}), \\ f &= k_0^2(\varepsilon_{zz}(\mathbf{r})-1)E_z^{inc} \end{aligned} } $$右端项仅在介质区域内非零（$\varepsilon_{zz}\neq 1$）。
类比三维：三维散射场方程 $\nabla\times(\nabla\times\mathbf{E}^{sca})-k_0^2\bar{\bar{\varepsilon}}_r\mathbf{E}^{sca} = k_0^2(\bar{\bar{\varepsilon}}_r-\bar{\bar{I}})\mathbf{E}^{inc}$，结构完全一致。
2.3 ABC 边界条件
二维外行波在 $e^{-j\omega t}$ 下：
$$ E_z^{sca} \sim \frac{e^{jk_0 r}}{\sqrt{r}}\,F(\theta) \quad (r\to\infty) $$二维 Sommerfeld 辐射条件（$e^{-j\omega t}$）：
$$ \lim_{r\to\infty}\sqrt{r}\left(\frac{\partial E_z^{sca}}{\partial r} - jk_0\,E_z^{sca}\right) = 0 $$一阶 ABC（在 $\Gamma_2 = \partial\Omega$ 上）：
$$ \frac{\partial E_z^{sca}}{\partial n} - jk_0\,E_z^{sca} = 0 \quad \text{on } \Gamma_2 $$即：
$$ \frac{\partial E_z^{sca}}{\partial n} + (-jk_0)\,E_z^{sca} = 0 $$与 Robin 条件的对应：
$$ \boxed{ \gamma = -jk_0, \quad q = 0 \quad \text{（TM，无磁，一阶 ABC）} } $$对比 $e^{+j\omega t}$ 约定：$\gamma = +jk_0$，符号相反。
2.4 TE 极化散射场公式
散射场方程（无磁）：
$$ -\frac{\partial}{\partial x}\left(\frac{1}{\varepsilon_{yy}(\mathbf{r})}\frac{\partial H_z^{sca}}{\partial x}\right) - \frac{\partial}{\partial y}\left(\frac{1}{\varepsilon_{xx}(\mathbf{r})}\frac{\partial H_z^{sca}}{\partial y}\right) - k_0^2 H_z^{sca} = \text{源项} $$TE 的源项推导更复杂（涉及 $\alpha$ 的不连续性），但形式为：
$$ f^{TE} = \frac{\partial}{\partial x}\left[\left(\frac{1}{\varepsilon_{yy}}-1\right)\frac{\partial H_z^{inc}}{\partial x}\right] + \frac{\partial}{\partial y}\left[\left(\frac{1}{\varepsilon_{xx}}-1\right)\frac{\partial H_z^{inc}}{\partial y}\right] $$TE 参数对应：
$$ \boxed{ \begin{aligned} \phi &= H_z^{sca}, \\ \alpha_x &= \frac{1}{\varepsilon_{yy}(\mathbf{r})}, \\ \alpha_y &= \frac{1}{\varepsilon_{xx}(\mathbf{r})}, \\ \beta &= -k_0^2, \quad \gamma = -jk_0 \end{aligned} } $$2.5 完整参数对应汇总表
| 参数 | 通用 Helmholtz | TM散射（无磁） | TE散射（无磁） | 3D矢量FEM |
|---|---|---|---|---|
| 未知量 | $\phi$ | $E_z^{sca}$ | $H_z^{sca}$ | $\mathbf{E}^{sca}$ |
| $\alpha_x$ | $\alpha_x$ | $1$ | $1/\varepsilon_{yy}(\mathbf{r})$ | $1/\mu_r=1$ |
| $\alpha_y$ | $\alpha_y$ | $1$ | $1/\varepsilon_{xx}(\mathbf{r})$ | — |
| $\beta$ | $\beta$ | $-k_0^2\varepsilon_{zz}(\mathbf{r})$ | $-k_0^2$ | $-k_0^2\bar{\bar{\varepsilon}}_r$ |
| $\gamma$ | $\gamma$ | $-jk_0$ | $-jk_0$ | $-jk_0$（ABC） |
| $q$ | $q$ | $0$ | $0$ | $0$ |
| $f$ | $f$ | $k_0^2(\varepsilon_{zz}-1)E_z^{inc}$ | 复杂（见上） | $k_0^2(\bar{\bar{\varepsilon}}_r-\bar{\bar{I}})\mathbf{E}^{inc}$ |
| 各向异性位置 | $\alpha$或$\beta$ | 仅$\beta$ | 仅$\alpha$ | 仅$T^e$（质量） |
| 刚度矩阵 | $[K]$ | 各向同性 | 各向异性 | 各向同性（无磁） |
| 质量矩阵 | $[M]$ | 各向异性（复$\beta$） | 各向同性 | 各向异性 |
补充：收束到无磁、电各向异性情形
这一部分沿着上面的通用框架继续，但把材料假设进一步收束到“无磁、仅电各向异性”，并据此重写后续方程、边界条件与误差分析。
一、从 Maxwell 方程推导二维标量方程
1.1 问题设置
无磁（$\mu_r = 1$），电各向异性介质：
$$ \bar{\bar{\varepsilon}}_r = \begin{pmatrix}\varepsilon_{xx}&0&0\\0&\varepsilon_{yy}&0\\0&0&\varepsilon_{zz}\end{pmatrix}, \quad \mu_r = 1 $$二维：$\partial/\partial z = 0$，时间约定 $e^{-j\omega t}$。
1.2 TM 极化推导
场分量：$\mathbf{E} = E_z\hat{z}$，$\mathbf{H} = H_x\hat{x} + H_y\hat{y}$
Faraday 定律（$\mu_r=1$）：
$$ \nabla\times(E_z\hat{z}) = j\omega\mu_0\mathbf{H} $$$$ \frac{\partial E_z}{\partial y}\hat{x} - \frac{\partial E_z}{\partial x}\hat{y} = j\omega\mu_0(H_x\hat{x} + H_y\hat{y}) $$$$ H_x = \frac{1}{j\omega\mu_0}\frac{\partial E_z}{\partial y}, \quad H_y = \frac{-1}{j\omega\mu_0}\frac{\partial E_z}{\partial x} $$Ampere 定律：
$$ \nabla\times\mathbf{H} = -j\omega\varepsilon_0\bar{\bar{\varepsilon}}_r\mathbf{E} $$$z$ 分量：
$$ \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} = -j\omega\varepsilon_0\varepsilon_{zz}E_z $$代入 $H_x, H_y$：
$$ \frac{\partial}{\partial x}\left(\frac{-1}{j\omega\mu_0}\frac{\partial E_z}{\partial x}\right) - \frac{\partial}{\partial y}\left(\frac{1}{j\omega\mu_0}\frac{\partial E_z}{\partial y}\right) = -j\omega\varepsilon_0\varepsilon_{zz}E_z $$$$ \frac{-1}{j\omega\mu_0}\left(\frac{\partial^2 E_z}{\partial x^2} + \frac{\partial^2 E_z}{\partial y^2}\right) = -j\omega\varepsilon_0\varepsilon_{zz}E_z $$两边乘以 $-j\omega\mu_0$：
$$ \frac{\partial^2 E_z}{\partial x^2} + \frac{\partial^2 E_z}{\partial y^2} + \omega^2\mu_0\varepsilon_0\varepsilon_{zz}E_z = 0 $$$$ \nabla^2 E_z + k_0^2\varepsilon_{zz}E_z = 0 $$写成标准 Helmholtz 形式：
$$ \boxed{ \begin{aligned} &-\frac{\partial}{\partial x}\left(1\cdot\frac{\partial E_z}{\partial x}\right) \\ &\quad - \frac{\partial}{\partial y}\left(1\cdot\frac{\partial E_z}{\partial y}\right) \\ &\quad + (-k_0^2\varepsilon_{zz})E_z = 0 \end{aligned} } \tag{TM} $$TM 结论：$\mu_r=1$ 时，$\alpha_x = \alpha_y = 1$（各向同性梯度！），各向异性仅通过 $\varepsilon_{zz}$ 进入 $\beta$ 项。
1.3 TE 极化推导
场分量：$\mathbf{H} = H_z\hat{z}$，$\mathbf{E} = E_x\hat{x} + E_y\hat{y}$
Ampere 定律：
$$ \nabla\times(H_z\hat{z}) = -j\omega\varepsilon_0\bar{\bar{\varepsilon}}_r\mathbf{E} $$$$ \frac{\partial H_z}{\partial y}\hat{x} - \frac{\partial H_z}{\partial x}\hat{y} = -j\omega\varepsilon_0(\varepsilon_{xx}E_x\hat{x} + \varepsilon_{yy}E_y\hat{y}) $$$$ E_x = \frac{1}{-j\omega\varepsilon_0\varepsilon_{xx}}\frac{\partial H_z}{\partial y}, \quad E_y = \frac{1}{j\omega\varepsilon_0\varepsilon_{yy}}\frac{\partial H_z}{\partial x} $$Faraday 定律（$\mu_r=1$）：
$$ \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y} = j\omega\mu_0 H_z $$代入：
$$ \frac{\partial}{\partial x}\left(\frac{1}{j\omega\varepsilon_0\varepsilon_{yy}}\frac{\partial H_z}{\partial x}\right) - \frac{\partial}{\partial y}\left(\frac{-1}{j\omega\varepsilon_0\varepsilon_{xx}}\frac{\partial H_z}{\partial y}\right) = j\omega\mu_0 H_z $$$$ \frac{1}{j\omega\varepsilon_0}\left[\frac{\partial}{\partial x}\left(\frac{1}{\varepsilon_{yy}}\frac{\partial H_z}{\partial x}\right) + \frac{\partial}{\partial y}\left(\frac{1}{\varepsilon_{xx}}\frac{\partial H_z}{\partial y}\right)\right] = j\omega\mu_0 H_z $$两边乘以 $j\omega\varepsilon_0$：
$$ \frac{\partial}{\partial x}\left(\frac{1}{\varepsilon_{yy}}\frac{\partial H_z}{\partial x}\right) + \frac{\partial}{\partial y}\left(\frac{1}{\varepsilon_{xx}}\frac{\partial H_z}{\partial y}\right) + k_0^2 H_z = 0 $$标准 Helmholtz 形式：
$$ \boxed{ \begin{aligned} &-\frac{\partial}{\partial x}\left(\frac{1}{\varepsilon_{yy}}\frac{\partial H_z}{\partial x}\right) \\ &\quad - \frac{\partial}{\partial y}\left(\frac{1}{\varepsilon_{xx}}\frac{\partial H_z}{\partial y}\right) \\ &\quad + (-k_0^2)H_z = 0 \end{aligned} } \tag{TE} $$TE 结论：各向异性通过 $\varepsilon_{xx}, \varepsilon_{yy}$ 进入 $\alpha_x, \alpha_y$（各向异性梯度！），$\beta$ 是各向同性的。
1.4 TM 与 TE 的对偶性
| TM：$E_z$ | TE：$H_z$ | |
|---|---|---|
| 各向异性影响刚度项？ | 否（$\alpha_x=\alpha_y=1$） | 是（$\alpha_x=1/\varepsilon_{yy}$，$\alpha_y=1/\varepsilon_{xx}$） |
| 各向异性影响质量项？ | 是（$\beta=-k_0^2\varepsilon_{zz}$） | 否（$\beta=-k_0^2$） |
| 各向同性退化 | $\alpha=1,\,\beta=-k_0^2\varepsilon_r$ | $\alpha=1/\varepsilon_r,\,\beta=-k_0^2$ |
二、散射场公式与边界条件
2.1 TM 散射场方程
总场 $E_z^{tot} = E_z^{inc} + E_z^{sca}$。
入射场（$e^{-j\omega t}$，沿 $\hat{k}=\cos\theta_i\hat{x}+\sin\theta_i\hat{y}$ 传播）：
$$ E_z^{inc} = E_0\,e^{jk_0(x\cos\theta_i+y\sin\theta_i)} $$入射场满足：$-\nabla^2 E_z^{inc} - k_0^2 E_z^{inc} = 0$
总场满足：$-\nabla^2 E_z^{tot} - k_0^2\varepsilon_{zz}(\mathbf{r})E_z^{tot} = 0$
相减：
$$ \boxed{ -\nabla^2 E_z^{sca} - k_0^2\varepsilon_{zz}(\mathbf{r})E_z^{sca} = k_0^2[\varepsilon_{zz}(\mathbf{r})-1]E_z^{inc} } $$2.2 TE 散射场方程
入射场：$H_z^{inc} = H_0\,e^{jk_0(x\cos\theta_i+y\sin\theta_i)}$
入射场满足：$-\nabla^2 H_z^{inc} - k_0^2 H_z^{inc} = 0$
总场满足 TE 方程。散射场方程：
$$ -\frac{\partial}{\partial x}\left(\frac{1}{\varepsilon_{yy}}\frac{\partial H_z^{sca}}{\partial x}\right) - \frac{\partial}{\partial y}\left(\frac{1}{\varepsilon_{xx}}\frac{\partial H_z^{sca}}{\partial y}\right) - k_0^2 H_z^{sca} $$$$ = \frac{\partial}{\partial x}\left[\left(\frac{1}{\varepsilon_{yy}}-1\right)\frac{\partial H_z^{inc}}{\partial x}\right] + \frac{\partial}{\partial y}\left[\left(\frac{1}{\varepsilon_{xx}}-1\right)\frac{\partial H_z^{inc}}{\partial y}\right] $$TE 源项涉及入射场的导数，这在弱形式中会通过分部积分化简。
2.3 ABC 边界条件（$e^{-j\omega t}$）
二维外行波 $\sim e^{jk_0 r}/\sqrt{r}$，Sommerfeld 条件给出一阶 ABC：
$$ \frac{\partial\phi^{sca}}{\partial n} - jk_0\,\phi^{sca} = 0 \quad \text{on } \Gamma_2 $$即 Robin 条件中：
$$ (\alpha_x\frac{\partial\phi}{\partial x}\hat{x}+\alpha_y\frac{\partial\phi}{\partial y}\hat{y})\cdot\hat{n} + \gamma\,\phi = q $$TM（$\alpha_x=\alpha_y=1$）：
$$ \frac{\partial E_z^{sca}}{\partial n} + \underbrace{(-jk_0)}_{\gamma}\,E_z^{sca} = \underbrace{0}_{q} $$TE（$\alpha_x=1/\varepsilon_{yy}$，$\alpha_y=1/\varepsilon_{xx}$）：
在ABC边界上（自由空间中），$\varepsilon_{xx}=\varepsilon_{yy}=1$，所以 $\alpha_x=\alpha_y=1$：
$$ \frac{\partial H_z^{sca}}{\partial n} + (-jk_0)\,H_z^{sca} = 0 $$ABC 条件在自由空间中，TM 和 TE 的 $\gamma$ 相同。
2.4 完整参数对应表
$$ \boxed{ \begin{array}{c|c|c|c} \text{参数} & \text{通用 Helmholtz} & \text{TM（无磁）} & \text{TE（无磁）} \\\hline \phi & \phi & E_z^{sca} & H_z^{sca} \\ \alpha_x & \alpha_x & 1 & 1/\varepsilon_{yy}(\mathbf{r}) \\ \alpha_y & \alpha_y & 1 & 1/\varepsilon_{xx}(\mathbf{r}) \\ \beta & \beta & -k_0^2\varepsilon_{zz}(\mathbf{r}) & -k_0^2 \\ \gamma & \gamma & -jk_0 & -jk_0 \\ q & q & 0 & 0 \\ f & f & k_0^2(\varepsilon_{zz}-1)E_z^{inc} & \text{TE源项（含导数）} \end{array} } $$三、弱形式与矩阵元素
3.1 TM 弱形式（完整展开）
$$ \begin{aligned} &\underbrace{\int_\Omega\left(\frac{\partial w}{\partial x}\frac{\partial E_z^{sca}}{\partial x}+\frac{\partial w}{\partial y}\frac{\partial E_z^{sca}}{\partial y}\right)d\Omega}_{[K]\text{：各向同性刚度}} \\ &\quad \underbrace{-\,k_0^2\int_\Omega\varepsilon_{zz}(\mathbf{r})\,w\,E_z^{sca}\,d\Omega}_{[M]\text{：各向异性质量}} \\ &\quad \underbrace{-\,jk_0\oint_{\Gamma_2}w\,E_z^{sca}\,d\Gamma}_{[R]\text{：ABC边界}} \\ &= \underbrace{k_0^2\int_{\Omega_d}(\varepsilon_{zz}-1)\,w\,E_z^{inc}\,d\Omega}_{\{b\}\text{：仅介质区域}} \end{aligned} $$3.2 TE 弱形式
$$ \begin{aligned} &\underbrace{\int_\Omega\left(\frac{1}{\varepsilon_{yy}}\frac{\partial w}{\partial x}\frac{\partial H_z^{sca}}{\partial x}+\frac{1}{\varepsilon_{xx}}\frac{\partial w}{\partial y}\frac{\partial H_z^{sca}}{\partial y}\right)d\Omega}_{[K]\text{：各向异性刚度}} \\ &\quad \underbrace{-\,k_0^2\int_\Omega w\,H_z^{sca}\,d\Omega}_{[M]\text{：各向同性质量}} \\ &\quad \underbrace{-\,jk_0\oint_{\Gamma_2}w\,H_z^{sca}\,d\Gamma}_{[R]\text{：ABC边界}} \\ &= \{b_{TE}\} \end{aligned} $$3.3 单元矩阵的物理意义
TM 极化的单元矩阵：
$$ K^e_{ji} = \frac{1}{4A^e}(b_j b_i + c_j c_i) \quad \text{（各向同性！与 $\varepsilon$ 无关）} $$$$ M^e_{ji} = -k_0^2\varepsilon_{zz}^{(e)}\cdot\frac{(1+\delta_{ji})A^e}{12} \quad \text{（$\varepsilon_{zz}$ 只在这里！）} $$TE 极化的单元矩阵：
$$ K^e_{ji} = \frac{1}{4A^e}\left(\frac{b_j b_i}{\varepsilon_{yy}^{(e)}} + \frac{c_j c_i}{\varepsilon_{xx}^{(e)}}\right) \quad \text{（各向异性在这里！）} $$$$ M^e_{ji} = -k_0^2\cdot\frac{(1+\delta_{ji})A^e}{12} \quad \text{（与 $\varepsilon$ 无关）} $$四、二维/三维方法的系统对比
4.1 方程层面
| 2D TM标量 | 2D TE标量 | 3D矢量 | |
|---|---|---|---|
| 方程 | $-\nabla^2 E_z - k_0^2\varepsilon_{zz}E_z = f$ | $-\nabla\cdot(\bar{\bar{\alpha}}\nabla H_z)-k_0^2 H_z=f$ | $\nabla\times(\nabla\times\mathbf{E})-k_0^2\bar{\bar{\varepsilon}}_r\mathbf{E}=\mathbf{f}$ |
| 未知数 | 1个标量 | 1个标量 | 3个矢量分量 |
| 单元 | 三角形 | 三角形 | 四面体 |
| 基函数 | 节点元（$\lambda_k$） | 节点元（$\lambda_k$） | 棱边元（$\lambda_a\nabla\lambda_b-\lambda_b\nabla\lambda_a$） |
| DOF/单元 | 3 | 3 | 6 |
| 伪解 | 无 | 无 | 节点元有，棱边元无 |
4.2 离散化层面
| 2D 节点元 | 3D 棱边元 | |
|---|---|---|
| 全局基函数限制 | $\mathcal{N}_I\big\|_{\Omega^e}=N_{i}^{(e)}$ | $\boldsymbol{\mathcal{N}}_I\big\|_{\Omega^e}=s_i^{(e)}\mathbf{N}_i^{(e)}$ |
| 方向符号 | 不需要 | $s_i^{(e)}=\pm 1$ |
| 组装公式 | $A_{JI}\mathrel{+}=K^e_{ji}$ | $A_{JI}\mathrel{+}=s_js_iK^e_{ji}$ |
| 边界积分 | 线段上 $2\times 2$ | 三角面上 $3\times 3$ |
| 边界DOF/段 | 2 | 3 |
五、完整误差分析
5.1 误差的来源
有限元解 $\phi_h$（或 $\mathbf{E}_h$）与精确解 $\phi$（或 $\mathbf{E}$）之间的误差来自四个独立的来源：
┌─────────────────────────────────────────────────────┐
│ 总误差 = e₁ + e₂ + e₃ + e₄ │
│ │
│ e₁：离散化误差（有限维逼近） │
│ e₂：计算域截断误差（ABC/PML） │
│ e₃：数值积分误差（高斯积分） │
│ e₄：线性方程组求解误差（有限精度/迭代） │
│ │
└─────────────────────────────────────────────────────┘
下面逐项分析。
5.2 离散化误差（最核心）
5.2.1 二维节点元（标量 Helmholtz）
函数空间：精确解 $\phi \in H^1(\Omega)$，有限元解 $\phi_h \in V_h \subset H^1(\Omega)$
其中 $V_h$ 是分片线性函数空间（$p=1$ 阶）。
Céa 引理：Galerkin 方法是"最佳逼近"（在适当范数下）：
$$ \|\phi - \phi_h\|_{H^1} \leq C\inf_{v_h\in V_h}\|\phi - v_h\|_{H^1} $$其中 $C$ 依赖于双线性形式 $a(\cdot,\cdot)$ 的连续性和强制性常数。
插值误差估计：设 $\Pi_h\phi$ 是 $\phi$ 的分片线性插值，则（标准结果）：
$$ \|\phi - \Pi_h\phi\|_{L^2(\Omega)} \leq C_1 h^2 |\phi|_{H^2(\Omega)} $$$$ \|\phi - \Pi_h\phi\|_{H^1(\Omega)} \leq C_2 h |\phi|_{H^2(\Omega)} $$其中 $h$ 是最大单元尺寸，$|\phi|_{H^k}$ 是 $H^k$ 半范数。
结合 Céa 引理：
$$ \boxed{ \|\phi - \phi_h\|_{H^1(\Omega)} \leq C\,h\,|\phi|_{H^2(\Omega)} \quad \text{（2D节点元，$p=1$阶）} } $$$L^2$ 误差（通过 Aubin-Nitsche 对偶论证提高一阶）：
$$ \boxed{ \|\phi - \phi_h\|_{L^2(\Omega)} \leq C\,h^2\,|\phi|_{H^2(\Omega)} \quad \text{（2D节点元，$p=1$阶）} } $$$p$ 阶一般结果（$p$ 阶多项式基函数）：
$$ \|\phi-\phi_h\|_{H^1} \leq C\,h^p\,|\phi|_{H^{p+1}}, \quad \|\phi-\phi_h\|_{L^2} \leq C\,h^{p+1}\,|\phi|_{H^{p+1}} $$5.2.2 三维棱边元（矢量 Helmholtz）
函数空间：精确解 $\mathbf{E} \in H(\text{curl};\Omega)$，有限元解 $\mathbf{E}_h \in V_h \subset H(\text{curl};\Omega)$
$H(\text{curl})$ 范数：
$$ \|\mathbf{E}\|_{H(\text{curl})} = \left(\|\mathbf{E}\|^2_{L^2} + \|\nabla\times\mathbf{E}\|^2_{L^2}\right)^{1/2} $$一阶 Nédélec 元的插值误差（$p=1$）：
$$ \|\mathbf{E} - \Pi_h\mathbf{E}\|_{L^2} \leq C_1\,h\,|\mathbf{E}|_{H^1} $$$$ \|\nabla\times(\mathbf{E}-\Pi_h\mathbf{E})\|_{L^2} \leq C_2\,h\,|\nabla\times\mathbf{E}|_{H^1} $$注意：棱边元的 $L^2$ 插值误差是 $O(h)$，而非节点元的 $O(h^2)$。这是因为棱边元的逼近空间不同。
结合 Céa 引理的 $H(\text{curl})$ 版本：
$$ \boxed{ \|\mathbf{E}-\mathbf{E}_h\|_{H(\text{curl})} \leq C\,h\,\left(|\mathbf{E}|_{H^1} + |\nabla\times\mathbf{E}|_{H^1}\right) \quad \text{（3D Nédélec $p=1$）} } $$$p$ 阶一般结果：
$$ \|\mathbf{E}-\mathbf{E}_h\|_{H(\text{curl})} \leq C\,h^p\,\left(|\mathbf{E}|_{H^p} + |\nabla\times\mathbf{E}|_{H^p}\right) $$5.2.3 二维/三维误差的系统对比
| 2D 节点元 ($p$阶) | 3D 棱边元 ($p$阶) | |
|---|---|---|
| 函数空间 | $H^1(\Omega)$ | $H(\text{curl};\Omega)$ |
| 范数 | $\|\phi\|_{H^1}=(\|\phi\|^2_{L^2}+\|\nabla\phi\|^2_{L^2})^{1/2}$ | $\|\mathbf{E}\|_{H(\text{curl})}=(\|\mathbf{E}\|^2_{L^2}+\|\nabla\times\mathbf{E}\|^2_{L^2})^{1/2}$ |
| 能量范数误差 | $O(h^p)$ | $O(h^p)$ |
| $L^2$ 误差 | $O(h^{p+1})$（Aubin-Nitsche） | $O(h^p)$（无超收敛！） |
| 正则性要求 | $\phi\in H^{p+1}$ | $\mathbf{E}\in H^p$，$\nabla\times\mathbf{E}\in H^p$ |
| $p=1$ 能量范数 | $O(h)$ | $O(h)$ |
| $p=1$ $L^2$ 范数 | $O(h^2)$ | $O(h)$ |
重要差异：节点元的 $L^2$ 误差比能量范数误差高一阶（超收敛），但棱边元没有这个超收敛。这是 $H(\text{curl})$ 空间的固有限制。
5.3 波数相关的误差——“污染效应”
5.3.1 Helmholtz 方程的特殊困难
对于 Helmholtz 方程（$\beta = -k_0^2\varepsilon$），双线性形式 $a(\cdot,\cdot)$ 不是强制的（coercive）。标准 Céa 引理给出的常数 $C$ 依赖于 $k_0$：
$$ C = C(k_0) \to \infty \quad \text{当 } k_0 \to \infty $$精确的误差估计（Babuška-Sauter, Melenk-Sauter 等）：
$$ \|\phi-\phi_h\|_{H^1} \leq C_1(1+k_0 h)^p\,\underbrace{(k_0 h)^p}_{\text{插值误差}} + C_2\underbrace{k_0(k_0 h)^{2p}}_{\text{污染误差}} $$物理解释：
- 插值误差 $(k_0 h)^p$：局部逼近误差，由每个波长内的单元数决定
- 污染误差 $k_0(k_0 h)^{2p}$：全局相位误差的累积，随传播距离增加而增长
5.3.2 “每波长单元数"准则
定义每波长单元数：
$$ N_\lambda = \frac{\lambda}{h} = \frac{2\pi}{k_0 h} $$则 $k_0 h = 2\pi/N_\lambda$，插值误差变为：
$$ \text{插值误差} \sim \left(\frac{2\pi}{N_\lambda}\right)^p $$工程准则：
| 精度要求 | $p=1$（线性元） | $p=2$（二次元） |
|---|---|---|
| 1% 误差 | $N_\lambda \geq 10$ | $N_\lambda \geq 5$ |
| 0.1% 误差 | $N_\lambda \geq 30$ | $N_\lambda \geq 7$ |
5.3.3 二维/三维的污染效应对比
| 2D | 3D | |
|---|---|---|
| DOF 总数 | $N \sim (k_0 L)^2\cdot N_\lambda^2$ | $N \sim (k_0 L)^3\cdot N_\lambda^3$ |
| 插值误差 | $(k_0 h)^p$ | $(k_0 h)^p$ |
| 污染误差 | $k_0 L\cdot(k_0 h)^{2p}$ | $(k_0 L)^2\cdot(k_0 h)^{2p}$ |
| 污染主导条件 | $k_0 L \gg 1$（电大问题） | $k_0 L \gg 1$（更严重） |
| 控制污染所需 $N_\lambda$ | $N_\lambda \sim (k_0 L)^{1/(2p-1)}$ | $N_\lambda \sim (k_0 L)^{2/(2p-1)}$ |
关键结论：三维问题的污染效应比二维更严重，因为波传播的距离更长、累积的相位误差更大。
数值示例：散射体尺寸 $L = 10\lambda$，$p=1$
| 2D | 3D | |
|---|---|---|
| $k_0 L$ | $20\pi\approx 63$ | $20\pi\approx 63$ |
| 控制污染所需 $N_\lambda$ | $\sim 63^{1/1} = 63$ | $\sim 63^{2/1} = 3969$（不可行！） |
| 若 $p=2$ | $\sim 63^{1/3}\approx 4$ | $\sim 63^{2/3}\approx 16$ |
工程启示：对电大问题，高阶元（$p\geq 2$）几乎是必须的，尤其在三维中。
5.3.4 色散误差的物理图像
有限元法中，数值波的传播速度与真实波不同，导致累积相位误差：
$$ \phi_{\text{数值}}(x) = e^{jk_h x}, \quad k_h = k_0(1+\delta_k) $$其中 $\delta_k$ 是相对色散误差：
二维线性节点元：
$$ \delta_k \approx \frac{1}{24}(k_0 h)^2 + O(k_0 h)^4 $$$$ \text{经过距离 } L \text{ 后的相位误差} = k_0 L\cdot\delta_k \approx \frac{k_0 L}{24}(k_0 h)^2 $$三维线性棱边元：
$$ \delta_k \approx \frac{1}{24}(k_0 h)^2 + O(k_0 h)^4 \quad \text{（系数可能不同，依赖于网格类型）} $$$$ \text{经过距离 } L \text{ 后的相位误差} \approx \frac{k_0 L}{24}(k_0 h)^2 $$例：$k_0 h = 2\pi/10$（每波长10个单元），$k_0 L = 20\pi$（10个波长）：
$$ \text{相位误差} \approx \frac{20\pi}{24}\left(\frac{2\pi}{10}\right)^2 \approx \frac{20\pi}{24}\times 0.395 \approx 1.03 \text{ rad} \approx 59° $$这是不可接受的！需要增加 $N_\lambda$ 或使用高阶元。
5.4 计算域截断误差（ABC/PML）
5.4.1 一阶 ABC 的反射误差
一阶 ABC 是一个近似——它只能精确吸收法向入射的平面波。
反射系数（平面波以角度 $\theta$ 入射到 ABC 边界）：
$$ R_{ABC}(\theta) = \frac{\cos\theta - 1}{\cos\theta + 1} $$| 入射角 $\theta$ | $\|R_{ABC}\|$ |
|---|---|
| $0°$（法向） | $0$（完美吸收） |
| $30°$ | $0.072$ |
| $45°$ | $0.172$ |
| $60°$ | $0.333$ |
| $80°$ | $0.673$ |
| $90°$（掠射） | $1.0$（完全反射） |
总的 ABC 误差：
$$ \boxed{ \|e_{ABC}\| \leq C\cdot\frac{1}{k_0 d}\cdot\left(\frac{a}{d}\right) } $$其中 $d$ 是 ABC 边界到散射体的距离，$a$ 是散射体尺寸。
5.4.2 二维/三维 ABC 误差对比
| 2D | 3D | |
|---|---|---|
| ABC 类型 | 曲线上 | 曲面上 |
| 一阶 ABC | $\partial\phi/\partial n - jk_0\phi = 0$ | $\hat{n}\times(\nabla\times\mathbf{E})=jk_0\hat{n}\times(\hat{n}\times\mathbf{E})$ |
| 反射系数 | $R(\theta) = \frac{\cos\theta-1}{\cos\theta+1}$ | 同（但立体角积分不同） |
| 远场衰减 | $\sim 1/\sqrt{r}$（2D柱面波） | $\sim 1/r$（3D球面波） |
| ABC 精度 | 较差（衰减慢） | 较好（衰减快） |
| 推荐距离 | $d \geq 0.5\lambda$ | $d \geq 0.25\lambda$ |
2D 的特殊困难：二维中散射波衰减为 $1/\sqrt{r}$，比三维的 $1/r$ 慢，因此需要更远的 ABC 边界或更高阶的 ABC。
5.4.3 PML 的误差
PML（完美匹配层）消除了角度依赖的反射问题，但引入了离散化误差：
PML 理论反射系数（连续情况）：
$$ R_{PML} = e^{-2jk_0\int_0^{d_{PML}}\sigma(\rho)\,d\rho/(j\omega\varepsilon_0)} \quad \text{（理论上为零，对所有角度）} $$离散化后的 PML 反射：
$$ |R_{PML,h}| \approx e^{-2\sigma_{max}d_{PML}/(ck_0)} + C(k_0 h)^{2p} $$第一项来自 PML 吸收不足（$\sigma_{max}$ 有限），第二项来自离散化误差。
| 2D | 3D | |
|---|---|---|
| PML 层数 | 通常 5-10 层单元 | 通常 5-10 层单元 |
| 额外 DOF | $\sim N_\lambda\times$(周长/h) | $\sim N_\lambda\times$(表面积/$h^2$) |
| 总 DOF 增加 | ~20-50% | ~30-80% |
5.5 数值积分误差
5.5.1 线性元的精确积分
关键观察：对于线性基函数，刚度矩阵的被积函数是常数，质量矩阵的被积函数是二次多项式。
2D 三角形：
| 矩阵 | 被积函数阶数 | 所需积分精度 | 1点高斯？ | 3点高斯？ |
|---|---|---|---|---|
| $K^e$（刚度） | 0（常数） | 精确 | ✓ 精确 | ✓ 精确 |
| $M^e$（质量） | 2（$\lambda_i\lambda_j$） | 2阶 | ✗ 不够 | ✓ 精确 |
| $f^e$（源项，常 $f$） | 1（$\lambda_j$） | 1阶 | ✓ 精确 | ✓ 精确 |
| $f^e$（源项，$f=e^{jk_0\cdot r}$） | 超越函数 | 需高阶 | ✗ | 近似 |
3D 四面体：
| 矩阵 | 被积函数阶数 | 所需积分精度 | 1点高斯？ | 4点高斯？ |
|---|---|---|---|---|
| $S^e$（旋度刚度） | 0（常数） | 精确 | ✓ 精确 | ✓ 精确 |
| $T^e$（质量） | 2 | 2阶 | ✗ 不够 | ✓ 精确 |
| $f^e$（源项，$e^{jk_0\cdot r}$） | 超越函数 | 需高阶 | ✗ | 近似 |
5.5.2 质量矩阵集中化（mass lumping）的误差
有时用对角化的质量矩阵代替一致质量矩阵：
一致质量矩阵（精确）：
$$ [M^e]_{\text{consistent}} = \frac{\beta A^e}{12}\begin{pmatrix}2&1&1\\1&2&1\\1&1&2\end{pmatrix} $$集中质量矩阵（近似）：
$$ [M^e]_{\text{lumped}} = \frac{\beta A^e}{3}\begin{pmatrix}1&0&0\\0&1&0\\0&0&1\end{pmatrix} $$集中化引入的误差：
$$ \|[M]_{\text{consistent}} - [M]_{\text{lumped}}\| = O(h^2) $$这不影响整体收敛阶数（对线性元），但可能增大常数。
5.5.3 源项积分的误差（波数相关）
对于含 $e^{jk_0(x\cos\theta+y\sin\theta)}$ 的源项，需要特别注意：
单元中心近似（最低阶）：
$$ \int_{\Omega^e}\lambda_j\,e^{jk_0\mathbf{k}\cdot\mathbf{r}}\,d\Omega \approx e^{jk_0\mathbf{k}\cdot\mathbf{r}_c}\int_{\Omega^e}\lambda_j\,d\Omega = e^{jk_0\mathbf{k}\cdot\mathbf{r}_c}\frac{A^e}{3} $$误差：$O(k_0 h)^2$
$n$ 点高斯积分：
$$ \text{积分误差} \sim \frac{(k_0 h)^{2n}}{(2n)!} $$推荐：
| $k_0 h$ | 推荐积分点数 |
|---|---|
| $< 0.3$（$N_\lambda > 20$） | 1点（中心）即可 |
| $0.3-0.6$（$N_\lambda \sim 10$） | 3-4点 |
| $0.6-1.0$（$N_\lambda \sim 6$） | 7点以上 |
5.6 线性方程组求解误差
5.6.1 条件数分析
Helmholtz 方程全局矩阵的条件数：
$$ \kappa([A]) = \frac{|\lambda_{max}|}{|\lambda_{min}|} $$| 2D 节点元 | 3D 棱边元 | |
|---|---|---|
| Laplace（$\beta=0$） | $\kappa\sim h^{-2}$ | $\kappa\sim h^{-2}$ |
| Helmholtz（$k_0\neq 0$） | $\kappa\sim k_0^2 h^{-2}$ 或更差 | $\kappa\sim k_0^2 h^{-2}$ 或更差 |
| 近共振 | $\kappa\to\infty$ | $\kappa\to\infty$ |
条件数对求解误差的影响：
$$ \frac{\|\delta\phi\|}{\|\phi\|} \leq \kappa([A])\cdot\frac{\|\delta b\|}{\|b\|} $$其中 $\delta b$ 是舍入误差。
5.6.2 直接法 vs 迭代法
| 方法 | 求解误差 | 适用规模 | 2D | 3D |
|---|---|---|---|---|
| LU分解（双精度） | $\sim 10^{-15}\kappa$ | $N<10^5$ | ✓ | 勉强 |
| GMRES（预条件） | 依赖收敛准则 | $N<10^7$ | ✓ | ✓ |
| 多重网格 | 依赖光滑性 | $N<10^8$ | ✓ | ✓ |
5.6.3 迭代法的收敛性
对于 Helmholtz 方程，标准迭代法（如无预条件的 GMRES）收敛困难：
$$ \text{迭代次数} \sim k_0^2\cdot N^{1/d} \quad \text{（无预条件）} $$| 2D ($d=2$) | 3D ($d=3$) | |
|---|---|---|
| 每次迭代复杂度 | $O(N)$ | $O(N)$ |
| 迭代次数 | $\sim k_0^2\sqrt{N}$ | $\sim k_0^2 N^{1/3}$ |
| 总复杂度 | $O(k_0^2 N^{3/2})$ | $O(k_0^2 N^{4/3})$ |
5.7 总误差的综合分析
5.7.1 各误差源的量级
设 $N_\lambda$ 为每波长单元数，$k_0 L$ 为电尺寸，$d$ 为 ABC 距离，$p$ 为基函数阶数：
$$ \boxed{ \begin{aligned} e_{total} &\leq \underbrace{C_1\left(\frac{2\pi}{N_\lambda}\right)^p}_{e_1:\text{插值误差}} + \underbrace{C_2\,k_0 L\left(\frac{2\pi}{N_\lambda}\right)^{2p}}_{e_1':\text{污染误差}} \\ &+ \underbrace{C_3\,\frac{a}{k_0 d^2}}_{e_2:\text{ABC误差}} + \underbrace{C_4\left(\frac{2\pi}{N_\lambda}\right)^{2n_g}}_{e_3:\text{积分误差}} + \underbrace{C_5\,\kappa\cdot\varepsilon_{machine}}_{e_4:\text{求解误差}} \end{aligned} } $$5.7.2 误差平衡原则
为使总误差最小，各项应量级相当（“平衡设计”）：
$$ e_1 \approx e_2 \approx e_3 \approx e_4 $$推荐参数选择（$p=1$ 线性元）：
| 参数 | 推荐值 | 依据 |
|---|---|---|
| $N_\lambda$ | $10\sim 20$ | $e_1 < 5\%$ |
| ABC 距离 $d$ | $\geq\lambda/2$ | $e_2 < 5\%$ |
| 高斯积分点数 | $3\sim 4$（2D），$4\sim 5$（3D） | $e_3 < e_1$ |
| 求解精度 | $10^{-10}$（直接法） | $e_4 \ll e_1$ |
5.7.3 二维/三维的完整误差对比表
| 误差源 | 2D（$p=1$） | 3D（$p=1$） | 备注 |
|---|---|---|---|
| 插值误差 | $O(h)$ 能量范数 | $O(h)$ $H(\text{curl})$范数 | 阶数相同 |
| $O(h^2)$ $L^2$ 范数 | $O(h)$ $L^2$ 范数 | 2D有超收敛，3D没有！ | |
| 污染误差 | $O(k_0 L\cdot(k_0 h)^2)$ | $O((k_0 L)^2\cdot(k_0 h)^2)$ | 3D更严重 |
| ABC 误差 | $O(1/(k_0 d))$ | $O(1/(k_0 d))$ | 类似，但2D衰减慢 |
| 积分误差 | 刚度精确，质量需3点 | 刚度精确，质量需4点 | |
| 条件数 | $O(k_0^2/h^2)$ | $O(k_0^2/h^2)$ | 类似 |
| DOF数（$N_\lambda$固定） | $O((k_0 L)^2)$ | $O((k_0 L)^3)$ | 维度灾难！ |
| 矩阵非零元（每行） | $\sim 7$ | $\sim 20$ | 3D更稠密 |
| 直接法复杂度 | $O(N^{3/2})$（嵌套剖分） | $O(N^2)$（嵌套剖分） | |
| 内存 | $O(N\log N)$ | $O(N^{4/3})$ |
5.7.4 计算量的具体估算
例：散射体 $10\lambda\times 10\lambda$（2D）或 $10\lambda\times 10\lambda\times 10\lambda$（3D），$N_\lambda=10$：
| 2D | 3D | |
|---|---|---|
| 计算域（含ABC） | $\sim 12\lambda\times 12\lambda$ | $\sim 12\lambda\times 12\lambda\times 12\lambda$ |
| 单元数 | $\sim 2\times 120^2 \approx 29,000$ | $\sim 6\times 120^3 \approx 10,000,000$ |
| DOF 数 | $\sim 14,000$（节点数） | $\sim 50,000,000$（棱边数） |
| 矩阵非零元 | $\sim 100,000$ | $\sim 10^9$ |
| 直接法时间 | $\sim$ 1秒 | 不可行 |
| 迭代法时间 | $\sim$ 0.1秒 | $\sim$ 数小时 |
| 内存 | $\sim$ 10 MB | $\sim$ 100 GB |
5.8 收敛性验证方法
5.8.1 Richardson 外推法
用三个不同网格密度 $h_1 > h_2 > h_3$ 计算，估计收敛阶数：
$$ p_{\text{observed}} = \frac{\ln\frac{\|\phi_{h_1}-\phi_{h_2}\|}{\|\phi_{h_2}-\phi_{h_3}\|}}{\ln\frac{h_1/h_2}{h_2/h_3}} $$预期结果：
| 2D（$L^2$） | 2D（$H^1$） | 3D（$L^2$） | 3D（$H(\text{curl})$） | |
|---|---|---|---|---|
| $p=1$ 元 | $p_{obs}\approx 2$ | $p_{obs}\approx 1$ | $p_{obs}\approx 1$ | $p_{obs}\approx 1$ |
| $p=2$ 元 | $p_{obs}\approx 3$ | $p_{obs}\approx 2$ | $p_{obs}\approx 2$ | $p_{obs}\approx 2$ |
5.8.2 收敛性验证的完整代码框架
import numpy as np
def solve_fem(N_lambda, params):
"""用给定的每波长单元数求解"""
# ... 网格生成、组装、求解 ...
return phi_h, mesh
def compute_error(phi_h, phi_exact, mesh, norm_type='L2'):
"""计算误差范数"""
error = 0.0
for e in mesh.elements:
if norm_type == 'L2':
# ∫|φ_h - φ_exact|² dΩ
error += gauss_integrate(
lambda r: abs(phi_h(r) - phi_exact(r))**2, e)
elif norm_type == 'H1':
# ∫|φ_h - φ_exact|² + |∇(φ_h - φ_exact)|² dΩ
error += gauss_integrate(
lambda r: abs(phi_h(r) - phi_exact(r))**2 +
norm(grad_phi_h(r) - grad_phi_exact(r))**2, e)
return np.sqrt(error)
# 收敛性测试
N_lambdas = [5, 10, 20, 40, 80]
errors_L2 = []
errors_H1 = []
hs = []
for N_lam in N_lambdas:
phi_h, mesh = solve_fem(N_lam, params)
h = 1.0 / N_lam # 波长归一化的单元尺寸
hs.append(h)
errors_L2.append(compute_error(phi_h, phi_exact, mesh, 'L2'))
errors_H1.append(compute_error(phi_h, phi_exact, mesh, 'H1'))
# 计算收敛阶
for i in range(1, len(hs)):
rate_L2 = np.log(errors_L2[i-1]/errors_L2[i]) / np.log(hs[i-1]/hs[i])
rate_H1 = np.log(errors_H1[i-1]/errors_H1[i]) / np.log(hs[i-1]/hs[i])
print(f"h={hs[i]:.4f}: L2 rate={rate_L2:.2f}, H1 rate={rate_H1:.2f}")
# 预期输出（2D线性节点元）：
# h=0.1000: L2 rate=1.95, H1 rate=0.98
# h=0.0500: L2 rate=2.01, H1 rate=1.00
# h=0.0250: L2 rate=2.00, H1 rate=1.00
# h=0.0125: L2 rate=2.00, H1 rate=1.00
六、误差分析总结
6.1 完整的误差层次图
总误差
│
┌───────────┼───────────┬──────────┐
│ │ │ │
离散化误差 截断误差 积分误差 求解误差
(最主要) (ABC/PML) (高斯) (机器精度)
│
┌─────┴─────┐
│ │
插值误差 污染误差
O(h^p) O(k₀L·(k₀h)^{2p})
(局部) (全局累积)
6.2 最终推荐
| 参数 | 2D 推荐 | 3D 推荐 |
|---|---|---|
| 基函数阶数 $p$ | 1（小问题），2（大问题） | 1（小问题），2（必须用于电大问题） |
| $N_\lambda$ | 10-20 | 10-20（$p=1$），6-10（$p=2$） |
| ABC/PML | PML 优先，5-10层 | PML 优先，5-10层 |
| 求解器 | 直接法（$N<10^5$） | 迭代法+预条件 |
| 验证 | 收敛性测试 + 光学定理 | 收敛性测试 + 光学定理 + 互易性 |