《Level Set Methods and Dynamic Implicit Surfaces》

Level Set的核心思想是用数量场中的等值面表示几何曲面，求解 Level Set 的实质，还是是用数值方法求解一个微分方程。既然涉及到空间微分方程，绕不开的一个问题就是数值微分。最简单的中心差分由于数值震荡的原因无法使用，迎风格式很简单很有效然而精度不高。Level Set 计算中最常用的 WENO 差分格式，然而这格式实在比较复杂，加之网上找到的资料多为流体方程求解所用，花了不少力气才看懂。

$$ \phi_t + \vec{V} \nabla{\phi} =0 $$

其中

$$ \vec{V} = \{ u,v,w\} $$

在一维情况下展开为：

$$ \frac{{\phi}^{n+1}-{{\phi}^n}}{{\Delta}t} + u^n{{\phi}^n_x} = 0$$

可见，求解方程首先需要计算$${{\phi}^n_x}$$，即第n时间步上$${{\phi}_x}$$的数值。 为简单，下面定义$${{\phi}_{i}} = {{\phi}_{x_i}}$$

CFL条件用于控制空间步长与时间步长，保持数值稳定性。一般认为，满足下式时，数值求解是稳定的，误差不会扩大：

$$ \Delta{t} < \frac{\Delta{x}}{max\{ |u| \}} $$

令$\alpha = \Delta{t} \frac{max\{ |u| \}}{\Delta{x}}$.一般，保守的取值是使得 $ \alpha = 0.5 $，激进的优化中可以取到 $ \alpha = 0.9 $.

定义：

$$ D^+{\phi} = {\phi}^+=\frac{{\phi}_{i+1}-{{\phi}_i}}{{\Delta}x}$$ $$ D^-{\phi} = {\phi}^-=\frac{{\phi}_{i}-{{\phi}_{i-1}}}{{\Delta}x}$$

迎风格式的思想是与物理实际相结合，运动方向是哪边来的，就利用哪边的导数信息（或者主要利用哪边的），故有：

$$ {\phi}_{i}=\left\{ \begin{aligned} {\phi}^-_i &\qquad& u_i > 0\\ {\phi}^+_i &\qquad& u_i < 0 \end{aligned} \right. $$

至此，一切都很简单。

Hamilton Jacobi ENO格式，精度比上述迎风格式高，可以到三阶精度。ENO全称Essentially Nonoscillatiry，精选非震荡，听着有点像广告…… 三阶精度不稀奇，难度在于保证光滑性与数值稳定性。保证数值稳定性的思想仍是迎风的思想，根据$u$的符号侧重使用某一边的数值信息，但是ENO格式提高了对${\phi}^+$或者${\phi}^-$的估计精度；对于光滑性的保证，ENO格式的原则是：在临近节点中选择选择可以使得结果最为光滑的点，作为插值模板，构建公式。

1.零阶，定义在网格点上： $$D^0_i{\phi}={\phi}_i$$ 2.一阶，定义在网格中间： $$D^1_{i+1/2}{\phi}=\frac{D^0_{i+1}{\phi}-D^0_{i}{\phi}}{\Delta{x}}$$ 3.二阶，定义在网格点上： $$D^2_{i}{\phi}=\frac{D^1_{i+1/2}{\phi}-D^1_{i-1/2}{\phi}}{2\Delta{x}}$$ 4.三阶，定义在网格点上： $$D^3_{i+1/2}{\phi}=\frac{D^2_{i+1}{\phi}-D^1_{i}{\phi}}{3\Delta{x}}$$

构建${\phi}(x)$的牛顿插值多项式如下： $${\phi}(x) = Q_0(x) + Q_1(x) + Q_2(x) + Q_3(x)$$ 在$x_i$处对上式两边取导数，则有： $${\phi}'(x_i) = {Q_1}'(x) + {Q_2}'(x) + {Q_3}'(x)$$

首先求一阶导数的估计，由于是 $\phi^-_x$的情况， 从$k=i-1$侧估计一阶导数，如下 ： $$ Q_1(x) = (D^1_{k+1/2})(x-x_i) $$ $$ {Q_1}'(x) = D^1_{k+1/2} $$ 下面看$Q_2$，从牛顿多项式的形式上来讲： $$ Q_2(x) = c(x-x_k)(x-x_{k+1}) $$ 两边求导： $$ {Q_2}'(x) = c[2(i-k)-1]\Delta{x} $$ 其中，$c$为二阶差分，但是二阶差分有两个，用哪一个呢？前面说过了：这取决于哪边的光滑……感觉好奇怪，但是情况就是这么个情况： $$c=\left\{ \begin{aligned} D^2_k\phi &\qquad& |{D^2_k\phi}| < |D^2_{k+1}\phi| \\ D^2_{k+1}\phi &\qquad& |{D^2_k\phi}| > |D^2_{k+1}\phi| \end{aligned} \right.$$ 可以看出来，上面所谓“从$k=i-1$侧”开始，其实是定义了一个类似程序指针的变量，定义了下面一步要计算哪一个节点的值。算完$Q_2$之后，需要再定义一个指针变量$k^*$留着后面用： $$k^*=\left\{ \begin{aligned} k-1 &\qquad& |{D^2_k\phi}| < |D^2_{k+1}\phi| \\ k &\qquad& |{D^2_k\phi}| > |D^2_{k+1}\phi| \end{aligned} \right.$$ 第三步，依葫芦画瓢，设$Q_3$的形式如下： $$ Q_3(x) = c^*(x-x_{k^*})(x-x_{k^*+1})(x-x_{k^*+2}) $$ 两边求导： $$ {Q_3}'(x) = c^*[3(i-k^*)^2-6(i-k^*)+2](\Delta{x})^2 $$ 其中： $$c*=\left\{ \begin{aligned} D^3_{k*+1/2}\phi &\qquad& |{D^3_{k*+1/2}\phi}| < |D^3_{k*+3/2}\phi| \\ D^3_{k*+3/2}\phi &\qquad& |{D^3_{k*+1/2}\phi}| > |D^3_{k*+3/2}\phi| \end{aligned} \right.$$ 至此，$Q_{1,2,3}$全部到手，可以计算$\phi^-_x$了

只有一句：开始的时候定义$k=i$即可。

ENO格式成于“essentially”精选，也缺陷于精选，因为“精选”就意味着丢弃。 例如，对于$\phi^-_{x_i}$，HJ_ENO格式使用了如下共计六个节点的数值：

$$\{{\phi_{i-3}},{\phi_{i-2}},{\phi_{i-2}},{\phi_{i}},{\phi_{i+1}},{\phi_{i+2}}\}$$

这些节点每一个都参与了计算，但是精选之后能留下来的有多少呢？

令： $$v_1=D^-{\phi_{i-2}}\qquad v_2=D^-{\phi_{i-1}}\qquad \\ v_3=D^-{\phi_{i}}\qquad v_4=D^-{\phi_{i+1}}\qquad v_5=D^-{\phi_{i+2}}\qquad$$ ENO格式的运算实质上是一个有两个分支的流程，对其整理可以发现可能的执行结果只有下面三种： $$ \begin{aligned} \phi^1_x &=& \frac{v_1}{3}-\frac{7v_2}{6}+\frac{11v_3}{6}\\ \phi^2_x &=& -\frac{v_2}{6} +\frac{5v_3}{6} + \frac{v_4}{3}\\ \phi^3_x &=& \frac{v_3}{3}+\frac{5v_4}{6}-\frac{v_5}{6} \end{aligned} $$ 两个分支下来，最后得到的只可能是这三种结果，详细的条件选择太烦，不写了，总之浪费了计算能力。 WENO格式解决这个问题的思路是：这么多点还是要的，否则不精确，但是计算结果不要丢掉，想想办法利用一下，不能减轻计算量提高点精度也是好的。于是他们做了一个加权平均，WENO之“W”即 weight 。

WENO的几个作者Chi-Wang Shu，Guang-Shan Jiang， Danping Peng，乍看都是中国人，不过署名学校是 University of California at Los Angeles以及Brown University。

对于证明细节也不细写，总之，在函数的光滑区间内，取$\omega_1 = 0.1$，$\omega_2 = 0.6$，$\omega_3 = 0.3$的时候，加权平均： $$ \phi_{x} = \omega_1\phi^1_x + \omega_2\phi^2_x + \omega_3\phi^3_x $$ 可以得到五阶精度的导数值估计。后面又有神人证明，若在每一个$\omega$上附加一项$O((\Delta{x})^2)$，即修改为$\omega_1 = 0.1 + O((\Delta{x})^2)$，$\omega_2 = 0.6 + O((\Delta{x})^2)$，$\omega_3 = 0.3 + O((\Delta{x})^2)$，精度仍然是五阶。 看起来特别简单。问题是，不光滑区间怎样处理？需要构造一种算法来调整，该算法原则如下：

• 在光滑区间内，保证$\omegas = \omega{s-origin} + O((\Delta{x})^2)$
• 在包含奇点的区间上，尽量消除包含奇点的计算模板的影响，也就相当于是自动选择三种ENO候选格式$\phi^i_x$中不受影响的那一个

显然，核心思想还是“找最光滑的”，那么首先要知道每个模板光滑程度如何。分别定义三个模板的光滑度： $$ \begin{aligned} S_1 &=& \frac{13}{12}(v_1-2v_2+v_3)^2 + \frac{1}{4}(v_1-4v_2+3v_3)\\ S_2 &=& \frac{13}{12}(v_2-2v_3+v_4)^2 + \frac{1}{4}(v_2-v_4)\\ S_3 &=& \frac{13}{12}(v_3-2v_4+v_5)^2 + \frac{1}{4}(3v_3-4v_4+v_5) \end{aligned} $$ 继续定义，中间变量： $$ \begin{aligned} \alpha_1 &=& \frac{0.1}{(S_1 + \epsilon)^2}\\ \alpha_2 &=& \frac{0.6}{(S_1 + \epsilon)^2}\\ \alpha_3 &=& \frac{0.3}{(S_1 + \epsilon)^2} \end{aligned} $$ 其中$\epsilon$的存在是为了防止出现被零除，其定义为： $$ \epsilon = 10^{-6}max \left\{ (v_1)^2,(v_2)^2,(v_3)^2,(v_4)^2,(v_5)^2 \right\} +10^{-99} $$ 现在，可以计算得到新的$\omega_i$： $$ \begin{aligned} \omega_1 &=& \frac{\alpha_1}{\alpha_1 + \alpha_2 + \alpha_3}\\ \omega_2 &=& \frac{\alpha_2}{\alpha_1 + \alpha_2 + \alpha_3}\\ \omega_3 &=& \frac{\alpha_3}{\alpha_1 + \alpha_2 + \alpha_3} \end{aligned} $$ 用这个加权平均即可。WENO格式除初始化（选择 $\phi^+_x$或者 $\phi^-_x$ ）以外，没有出现分支，一则不浪费计算能力，二则便于GPU并行化计算（GPU计算在某些时候需要尽量避免程序分支）。 最后，上面计算的都是$\phi^-_x$，对于$\phi^-_x$，将$v_i$定义的网格点向正方形平移即可，如图：

如图，(a)图为函数图像，(b)图为ENO的“权重取值”，只有0和1，因为ENO格式只取一个，©图为WENO算法下的权重。后面两图中，○+×分别代表123项的权重。可以看出，在多数情况下，WENO的权重很接近0.1、0.6、0.3，但是在接近奇点处会自动调整为更合适的值，保证数值计算的稳定性。