Engquist-Majda吸收边界

是之前都是对自由空间中的进行分析，在实际中由于计算量以及内存的问题。必须设置边界。这很好想清楚的（虽然想了很久）。

我们推导得到了不考虑边界的三维的FDTD的基本公式，可以很容易看出我们想得到某点的某时刻的值，可以通过前一个时刻周围点的值去推出该点该时刻的值。所以如果我们有初始时刻的场的分布，那就可以一直推下去。如果是自由空间，那我们就要存储无限个Yee元胞，这怎么可能我们不可能有无限的内存以及计算资源。比如：要计算a点，就要知道a的左右两点的前一个时刻的值，同理嘛，2->4->8->16->32….. 所以想模拟全空间的传播（亦或者说是全空间我们想要弄清楚的空间的传播）就要设置边界，形成一个有限空间，所有YEE元胞的值都能通过这个有限空间的值递推，不需要借助该空间以外的值。这个边界，一定要首先满足的是不会有反射波的产生，不然就会破坏空间里电磁波原本的传播了。所以它就叫吸收边界。

Engquist-Majda吸收边界的解析形式（Mur吸收边界）

这个边界有点像是截断边界，边界以外的区域的场值为0，即该边界将波全吸收了。（可能

考虑二维问题，对无源区域，电磁场满足亥姆霍兹方程（Helmholtz equation）。

$\frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} - \frac{1}{v^2} \frac{\partial^2 f}{\partial t^2} = 0 \tag {1-1}$

由分离变量法可解得，其解系

$\left\{ A_n \exp{[ \jmath (\omega t - k^n_x x - k^n_y y) ]} \right\} \tag {1-2}$

我们选择一个正方形区域：

$\begin{align} 0 \le x \le x_0 \\ 0 \le y \le y_0 \\ \end{align}$

取左边界$x = 0$可以得知，应存在入射波和反射波，则在区域内解应为入射波和反射波的合成。

取左边界部分观察：

$f(x) = A_- \exp{[ \jmath (\omega t + k_x x + k_y y) ]} + A_+ \exp{[ \jmath (\omega t - k_x x + k_y y) ]} ; \quad k^2 = k_x^2 + k_y^2 \; k = \omega / v \tag {1-3}$

其中入射波为：

$f_- (x) = A_- \exp{[ \jmath (\omega t + k_x x + k_y y) ]} \tag {1-4}$

反射波为：

$f_+ (x) = A_+ \exp{[ \jmath (\omega t - k_x x + k_y y) ]} \tag {1-5}$

把式（1-3）代入（1-1）保留对x的导数得：

$\begin{align} \frac{\partial^2 f}{\partial x^2} - k_y^2 f + k^2 f & = 0 \\ \left[ \frac{\partial^2 }{\partial x^2} + (k^2 - k_y^2) \right] f & = 0 \\ (\frac{\partial }{\partial x} - \jmath \sqrt{ k^2 - k_y^2 }) (\frac{\partial }{\partial x} + \jmath \sqrt{ k^2 - k_y^2 }) f & = 0 \tag {1-6} \\ \end{align}$

定义算子$L$：

$L = (\frac{\partial }{\partial x} - \jmath \sqrt{ k^2 - k_y^2 }) (\frac{\partial }{\partial x} + \jmath \sqrt{ k^2 - k_y^2 }) = L_- L_+ \tag {1-7}$

可有：

$\begin{align} Lf = 0 \\ L_- L_+ (f_- + f_+) = 0 \\ L_- L_+ f_- + L_- L_+ f_+ = 0 \tag {1-8} \\ \end{align}$

显然带入计算可得：

$L_+ f_+ = 0 , \quad L_- f_- = 0$

考虑模拟全空间（自由空间），边界处不应该有反射波。对式（1-8）做如下处理：

$L_+ (L_- f_-) + L_+ (L_- f_+) = L_+ (L_- f_+) = 0 \tag {1-9}$

如要不存在反射波，则有$f_+ = 0$。对式（1-9）有：

$L_- f_+ = 0 \Rightarrow \left. \frac{\partial f}{\partial x} - \jmath \sqrt{ k^2 - k_y^2 } f \right |_{x = 0} = 0 \tag {1-10}$

对于时谐场有算子替换：

$\frac{1}{v} \frac{\partial }{\partial t} \leftrightarrow \jmath k \tag {1-11}$

对左边界部分的反射波有算子替换：

$\frac{\partial }{\partial y} \leftrightarrow j k_y \tag {1-12}$

对式（1-9）有：

$\begin{align} \left. \frac{\partial f}{\partial x} - \sqrt{ \jmath^2 k^2 - \jmath^2 k_y^2 } f \right |_{x = 0} = 0 \\ \left. \frac{\partial f}{\partial x} - \sqrt{ \frac{1}{v^2} \frac{\partial^2 }{\partial t^2} - \frac{\partial^2 }{\partial y^2} } f \right |_{x = 0} = 0 \\ \end{align} \tag {1-13}$

同理对右边界$x = x_0$有：

$\left. \frac{\partial f}{\partial x} + \sqrt{ \frac{1}{v^2} \frac{\partial^2 }{\partial t^2} - \frac{\partial^2 }{\partial y^2} } f \right |_{x = x_0} = 0 \tag {1-14}$

同理可得上下边界

一阶近似

边界的解析式中含有根号内的求导运算，不会近似。恰好，学过Taylor展开近似。

由Taylor展开：

$\sqrt{ \jmath^2 k^2 - \jmath^2 k_y^2 } = \jmath k \sqrt{ 1 - \frac{k_y^2}{k^2} } \approx \jmath k - \jmath \frac{1}{2} \frac{k_y^2}{k} + \cdots \tag {1-15}$

对于一阶近似，我们只取$\jmath k$，则式（1-13）有：

$\left. \frac{\partial f}{\partial x} - \jmath k f \right |_{x = 0} = \left. \left( \frac{\partial f}{\partial x} - \frac{1}{v} \frac{\partial f}{\partial t} \right) \right |_{x = 0} = 0 \tag {1-16}$

二阶近似

对于二阶近似，我们只取$\jmath k - \jmath \frac{1}{2} \frac{k_y^2}{k}$，则式（1-13）有：

$\begin{align} \left. \frac{\partial f}{\partial x} - \jmath k f + \frac{\jmath }{2} \frac{k_y^2}{k} f \right |_{x = 0} & = 0 \\ \left. \left( k \frac{\partial f}{\partial x} - \frac{k}{v} \frac{\partial f}{\partial t} - \frac{\jmath}{2} \frac{\partial^2 f}{\partial y^2} \right) \right |_{x = 0} & = 0 \\ \left. \left( \frac{1}{v} \frac{\partial^2 f}{\partial t \partial x} - \frac{1}{v^2} \frac{\partial^2 f}{\partial t^2} + \frac{1}{2} \frac{\partial^2 f}{\partial y^2} \right) \right |_{x = 0} & = 0 \\ \end{align} \tag {1-17}$

汇总

我们可以类似二维推算出一维、三维。

一维Mur吸收边界的解析式

很容易得出，一维无需近似，其解析式为：

$\begin{cases} \left. \left( \frac{\partial f}{\partial z} - \frac{1}{v} \frac{\partial f}{\partial t} \right) \right |_{z = 0} = 0 \\ \left. \left( \frac{\partial f}{\partial z} + \frac{1}{v} \frac{\partial f}{\partial t} \right) \right |_{z = z_0} = 0 \\ \end{cases} \tag {1-18}$

二维Mur吸收边界的解析式

有：

$\begin{array}{c} % 总表格 \begin{array}{cc} % 第一行内分成两列 \begin{array}{|c|c} % 第一列"最小值"数组 & 一阶近似 & 二阶近似 \\ \hline x = 0 & \frac{\partial f}{\partial x} - \frac{1}{v} \frac{\partial f}{\partial t} & \frac{1}{v} \frac{\partial^2 f}{\partial t \partial x} - \frac{1}{v^2} \frac{\partial^2 f}{\partial t^2} + \frac{1}{2} \frac{\partial^2 f}{\partial y^2} = 0 & \\ \hline x = x_0 & \frac{\partial f}{\partial x} + \frac{1}{v} \frac{\partial f}{\partial t} & \frac{1}{v} \frac{\partial^2 f}{\partial t \partial x} + \frac{1}{v^2} \frac{\partial^2 f}{\partial t^2} - \frac{1}{2} \frac{\partial^2 f}{\partial y^2} = 0 & \\ \hline y = 0 & \frac{\partial f}{\partial y} - \frac{1}{v} \frac{\partial f}{\partial t} & x替换为y即可 \\ \hline y = y_0 & \frac{\partial f}{\partial y} + \frac{1}{v} \frac{\partial f}{\partial t} & x替换为y即可 \\ \end{array} \end{array} % 第一行表格组结束 \end{array} % 总表格结束$

三维Mur吸收边界的解析式

不难得出，三维Mur吸收边界的解析式的一阶近似式等同于二维的形式，对于其二阶近似式，仅给出截断面为yoz时的形式，其他同理即可：

$\begin{cases} \frac{1}{v} \frac{\partial^2 f}{\partial t \partial x} - \frac{1}{v^2} \frac{\partial^2 f}{\partial t^2} + \frac{1}{2} \left( \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} \right)f = 0 \quad ,x = 0 \\ \frac{1}{v} \frac{\partial^2 f}{\partial t \partial x} + \frac{1}{v^2} \frac{\partial^2 f}{\partial t^2} - \frac{1}{2} \left( \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} \right)f = 0 \quad ,x = x_0 \\ \end{cases} \tag {1-19}$

Mur吸收边界的FDTD形式

这里我只推二维，因为三维运算量太大了，暂时不动笔。而一维其实就是二维的一阶近似罢了。

考虑TM波，$x = 0$

TM波的$H_x$和$H_y$分量，其左边界量为$(i, j + 1/2)$和$(i+1/2, j )$，由上一节FDTD基本公式可知：

$H^{n+1/2}_x(i,j+1/2) = CP(i,j+1/2) H^{n-1/2}_x(i,j+1/2) - CQ(i,j+1/2) \frac{E^n_z(i+1,j) - E^n_z(i,j)}{\Delta y} \tag {1-20}$

所要求的值不涉及截断边界以外的值，所以无需使用吸收边界，同理$H_y$。这里想一下如果左边界$E_x(i,j)$仍使用FDTD基本公式求解，可以发现其要求使用到截断边界以外的值，所以对$E_z$要设置截断边界。

一阶近似式

对$E_z$，取$E^{n+1/2}(i+1/2, j)$为观察点。有：

$\frac{ E_z^{n+1/2}(i+1,j) - E_z^{n+1/2}(i,j) }{\Delta x} = \frac{E_z^{n+1}(i+1/2,j) - E_z^n(i+1/2,j)}{v \Delta t} \tag {1-22}$

这里由于与$j$无关我们简化一下：$E_z(i+1,j) \leftarrow f(i+1)$

采用线性插值

$v \Delta t \left[ f^{n+1}(i+1) + f^{n}(i+1) - f^{n+1}(i) - f^{n}(i) \right] = \Delta x \left[ f^{n+1}(i+1) + f^{n+1}(i) - f^n(i+1) -f^n(i) \right] \tag {1-22}$

可得：

$E_z^{n+1}(i, j) = E_z^{n}(i+1, j) + \frac{v \Delta t - \Delta x}{v \Delta t + \Delta x}\left[ E_z^{n+1}(i+1,j) - E_z^{n}(i,j) \right] \tag {1-23}$

二阶近似式

这里我们如果直接进行差分，运算有点大，三维情况下对二阶近似式就是直接差分，同时求出的结果也很直观的反映的是那些位置对边界产生了影响，二维也是如此，二维是6个点对其产生了影响，如果画出图，就很容易知道，为什么对于角顶点还要单独讨论。

这里我们不进行直接差分，对二阶近似式中的第三项。这里我们讨论对于无耗介质$\sigma_m = 0$

$\begin{aligned} & \frac{\partial^2 E_z}{\partial y^2} \\ = & \frac{\partial }{\partial y} \left( -\mu \frac{\partial H_x}{\partial t} -\sigma_m H_x \right) \\ = & - \mu \frac{\partial^2 H_x}{\partial t \partial y} \\ \end{aligned} \tag {1-24}$

将式（1-24）带回$\frac{1}{v} \frac{\partial^2 f}{\partial t \partial x} - \frac{1}{v^2} \frac{\partial^2 f}{\partial t^2} + \frac{1}{2} \frac{\partial^2 f}{\partial y^2} = 0$有：

$\begin{aligned} & \frac{1}{v} \frac{\partial^2 f}{\partial t \partial x} - \frac{1}{v^2} \frac{\partial^2 f}{\partial t^2} + \frac{1}{2} \frac{\partial^2 f}{\partial y^2} = 0 \\ & \frac{1}{v} \frac{\partial^2 f}{\partial t \partial x} - \frac{1}{v^2} \frac{\partial^2 f}{\partial t^2} - \mu \frac{\partial^2 H_x}{\partial t \partial y} = 0 \\ & \frac{\partial f}{\partial x} - \frac{1}{v} \frac{\partial f}{\partial t} - v\mu \frac{\partial H_x}{\partial y} = 0 \\ \end{aligned} \tag {1-25}$

头两项依旧是一阶近似式

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

插一句，同理右边界：

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

对式（1-25）右边第三项进行差分，取$H^{n+1/2}(i+1/2, j+1/2)$为观察点，再进行线性插值：

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

可得二阶近似：

$\begin{aligned} & E_z^{n+1}(i, j) \\ = & E_z^{n}(i+1, j) + \frac{v \Delta t - \Delta x}{v \Delta t + \Delta x}\left[ E_z^{n+1}(i+1,j) - E_z^{n}(i,j) \right] \\ & + \frac{v^2 \mu \Delta t \Delta x}{\Delta y (v \Delta t + \Delta x)} \\ & * \left[H_x^{n+1/2}(i+1,j+1/2) + H_x^{n+1/2}(i,j+1/2) - H_x^{n+1/2}(i+1,j-1/2) - H_x^{n+1/2}(i,j-1/2) \right] \end{aligned} \tag {1-29}$