双曲守恒律的ENO格式和WENO格式

问题描述

考虑下述Burgers方程的初值问题

分别使用有限体积法4阶ENO格式，有限体积法3阶和5阶WENO方法求解。时间方向用三阶TVD
Runge-Kutta方法。研究格式在解光滑情况下和有间断的情况下的数值精度并作图。

关键算法和格式描述

精确解

为了计算数值精度，必须要有精确解。对于精确解的产生，我们采用两种方式。

当时间比较短，解还比较光滑时，特征线交错情况比较简单时，我们使用特征线方法求解精确解。
对于上述黎曼问题的精确解，我们用特征线方法来求解精确值。过某个初值点(x0,0)$(x_0,0)$的特征线是一根直线，它的斜率为u0(x)$u_0(x)$，那么给定一点(x1,t1)$(x_1,t_1)$，我们知道它满足，

求得x0$x_0$后，我们可以得到在(x1,t1)$(x_1,t_1)$点处的值为u0(x0)$u_0(x_0)$。
需要注意的是，matlab的solve函数主要是用于多项式求根，对于非线性方程来说，也可以求，但是往往只给出一个最小解。在matlab里面solve命令主要是用来求解代数方程（即多项式）的解,但是也不是说其它方程一个也不能解，不过求解非代数方程的能力相当有限，通常只能给出很特殊的实数解。

第二种方法，精确解采用5阶WENO格式的小步长逼近作为相对近似精确解，这也只是相对而言的精确解。将精确解保存为mat 格式调用（不同步长取对应点）。

5阶WENO格式

空间上采用五阶WENO格式，时间上采用三阶Runge Kutta格式，具体操作描述如下：

空间离散

将空间定义域等分，取每个等分区间上中点的初值作为该区间的初值（不是计算等分节点的值）。那么，N等分就有N个值，对这N个值编号，不妨记为ui$u_i$，使用WENO格式计算下一时间层的值。

Lax-Friedrichs通量分裂

使用Lax-Friedrichs，进行通量分裂，如下：

其中，α=max(|f′(ui)|)$\alpha =\text{max}(|f^{\prime }(u_i)|)$，这个值随着ui$u_i$的变化而变化。定义左值u+j+1/2=12(f(uj)+αuj)$u_{j+1/2}^{+}=\frac{1}{2}(f(u_j)+\alpha u_j)$，定义右值u−j+1/2=12(f(uj+1)−αuj+1)$u_{j+1/2}^{-}=\frac{1}{2}(f(u_{j+1})-\alpha u_{j+1})$，下面要做的就是对这些左右值进行重构更新，再合成通量，求得空间算子。

原函数方法重构左值右值

下面以右值u−j+1/2$u_{j+1/2}^{-}$的重构更新为例（左值的重构与右值呈镜像对称（系数等），不再详述），来说明重构过程。为描述方便，用vi$v_i$来表示右值集合u−j+1/2$u_{j+1/2}^{-}$。

• 选择不同的模板，计算不同模板下的重构值，采用如图所示的模板，计算得的对应模板的三个重构值如下：

  

  v(0)iv(1)iv(2)i===(2vi−2−7vi−1+11vi)/6(−vi−1+5vi+2vi+1)/6(2vi+5vi+1−vi+2)/6$$\begin{array}{rlr}{v}_i^{(0)}& =& (2v_{i-2}-7v_{i-1}+11v_i)/6\\ v_i^{(1)}& =& (-v_{i-1}+5v_i+2v_{i+1})/6\\ v_i^{(2)}& =& (2v_i+5v_{i+1}-v_{i+2})/6\end{array}$$

  

• 计算三个模板得到的重构值对应的权重，逆向罗列公式如下（k=3$k=3$）：
  

  wrαr==αr∑r=0k−1αrdr(10−6+βr)2$$\begin{array}{rlr}{w}_r& =& \frac{{\alpha }_r}{\sum _{r=0}^{k-1}{\alpha }_r}\\ {\alpha }_r& =& \frac{d_r}{({10}^{-6}+{\beta }_r{)}^2}\end{array}$$

  对于五阶WENO格式，有：d0=0.3,d1=0.6,d2=0.1$d_0=0.3,d_1=0.6,d_2=0.1$，且βr${\beta }_r$可表达为：
  

  β0β1β2===13(vi−2vi−1+vi−2)212+(3vi−4vi−1+vi−2)24(vi−1−vi+1)24+13(vi−1−2vi+vi+1)21213(vi−2vi+1+vi+2)212+(3vi−4vi+1+vi+2)24$$\begin{array}{rlr}{\beta }_0& =& \frac{13{\left(v_i-2v_{i-1}+v_{i-2}\right)}^2}{12}+\frac{{\left(3v_i-4v_{i-1}+v_{i-2}\right)}^2}{4}\\ {\beta }_1& =& \frac{{\left(v_{i-1}-v_{i+1}\right)}^2}{4}+\frac{13{\left(v_{i-1}-2v_i+v_{i+1}\right)}^2}{12}\\ {\beta }_2& =& \frac{13{\left(v_i-2v_{i+1}+v_{i+2}\right)}^2}{12}+\frac{{\left(3v_i-4v_{i+1}+v_{i+2}\right)}^2}{4}\end{array}$$

• 更新最后的左值为各个值的加权平均：即u+i+1/2=v^i=∑r=0k−1wrv(r)i$u_{i+1/2}^{+}={\hat{v}}_i=\sum _{r=0}^{k-1}{w}_r{v}_i^{(r)}$。

• 同相同的方法重构左值。原来从左到右计算的，重构左值时从右到左。

空间算子和时间三阶Runge Kutta方法

通过如上方法，计算得到更新后的左值和右值，取左值和右值的和为通量函数近似，即：

那么，可以得到空间离散算子的作用结果，即：

有了空间离散算子，我们就在时间上使用三阶Runge-Kutta方法，求解下一时间层的值。

3阶WENO格式

三阶的WENO格式（k=2）和五阶的WENO格式（k=3）操作方法是类似的，所不同的是，在重构通量左右值时所用的模板和系数是不同的。一般来说，对于2k−1$2k-1$阶的WENO格式，所用的模板（具有对称性），共有k$k$个：

其上重构值的系数（u(r)i+1/2=∑j=0k−1Crjui−r+j$u_{i+1/2}^{(r)}=\sum _{j=0}^{k-1}{C}_{rj}{u}_{i-r+j}$），有人总结成了如下表格（右值）。

且对于k比较小时，dr,βr$d_r,{\beta }_r$有如下结果：

4阶ENO格式

ENO格式和WENO格式的过程大部分是类似的，不同的是，ENO格式左右值的重构，不需要计算全部模板下的重构值后取加权平均，而只要选取其中的某一个模板作为计算左右值重构即可，选取的准则为使差商最小。具体的算法描述如下。

模板的重构系数Crj$C_{rj}$同前表格（k = 4）。

计算结果展示和讨论

三种格式对于光滑解数值精度

五阶WENO格式

对于一个空间k阶格式，时间上是三阶龙格库塔格式，我们得保证，(Δt)3(Δx)k$\frac{(\mathrm{\Delta }t{)}^3}{(\mathrm{\Delta }x{)}^k}$是一个常值，或者让时间因素相比于空间因素，对结果造成的很小。换言之，若固定了空间步长，我们可以这样选取时间步长：

这里，Constant是一个常数，它的选取应该使得时间空间步长满足CFL条件：Δt/Δx∗α=CFL,α=max|f′(u)|$\mathrm{\Delta }t/\mathrm{\Delta }x\ast \alpha =CFL,\alpha =\text{max}|f^{\prime }(u)|$。

选取终止时间为0.1，解依然是连续的。此时，空间上使用五阶WENO格式，时间上采用三阶Runge
Kutta格式，进行了数值实验，最后得到了如下结果，时间关系，程序上的相关的参数设定和细节处理就不再详述了，见程序。得到的结果如下表：

三阶WENO格式

仿照五阶WENO格式，依然选择终止时间为0.1，得到的结果如下表：

收敛阶好像不是太稳定，我也不知道为什么呀。

四阶ENO格式

四阶ENO格式，在解连续的情况下，得到的结果如下表：有些小问题。为了减小时间方向带来的影响，设置Courant系数尽量小。

三种格式对于间断解的数值精度

从数值格式中，并没有看到解随着时间的推移变得间断（lusongno1.cn），所以这个间断是什么意思，没明白。下面我假设在t足够大时，如t=3$t=3$，解开始断开。我们来研究解在有间断的情况下的数值精度。这里的真解采用WENO5格式在空间步长为2000的情况下的解。

五阶WENO格式

三阶WENO格式

四阶ENO格式

不知道是否精确解求解有问题，以上三种格式在解间断时只能达到一阶精度。

程序源码和其他话

源代码见附件，使用的语言为MAMTLAB2014a（盗版）,Windows 10环境。

一些小体会：

• Matlab巧用shiftcirt函数，来处理周期边界问题是特别棒的事情。

• 知道如何进行右值重构后，再进行计算左值重构，我们没必要考虑系数的镜像对称关系。一个简单的做法是，将左值反向排列，然后同样用计算右值的程序计算一遍，得到的结果最后再转个方向即可。

• CFL迎风条件： dt/dx * alpha = CFL

• 对于均匀网格，使用重构系数表Crj$C_{rj}$，重构值直接用已有值使用系数表线性表出，而省去了插值的过程，提高计算速度，代价是泛化能力差。也就是说k变化时，你得重新写一遍系数。

• 初值点的选取，选的不是等分节点处，而是选取单元中点处。使用单元中点值表示这个单元的值。

• 对这些格式，在有些步长下能得到相应的收敛阶，而在有的步长下不能。应该考虑时间步长和空间步长谁是主导因素，主导因素是否会发生转移。

参考文献
Shu C W. Numerical methods for hyperbolic conservation
laws[J]. Lecture notes (AM257), 2006.
Shu C W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[M]//Advanced numerical approximation of nonlinear hyperbolic equations. Springer Berlin Heidelberg, 1998:
325-432.
Jiang G S, Shu C W. Efficient implementation of weighted ENO
schemes[J]. Journal of computational physics, 1996, 126(1): 202-228.

---