额，双曲守恒律

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

问题描述

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

$$\label{q} \left \{ \begin{aligned} & u_t + (u^2/2)_x = 0 & -1\leq x \leq 1,t>0\\ & u(x,0) = u_0(x) = \text{sin}(\pi x ) & \\ &periodic ~~B.C. ~~for ~~x \end{aligned} \right.$$

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

关键算法和格式描述

精确解

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

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

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

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

5阶WENO格式

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

空间离散

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

Lax-Friedrichs通量分裂

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

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

原函数方法重构左值右值

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

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

  

  $$\begin{aligned} 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{aligned}$$

  

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

  $$\begin{aligned} w_r &=& \frac{\alpha_r}{\sum\limits_{r=0}^{k-1}{\alpha_r}} \\ \alpha_r &=& \frac{d_r}{(10^{-6} + \beta_r)^2} \end{aligned}$$

  对于五阶WENO格式，有：$d_0 = 0.3,d_1 = 0.6,d_2 = 0.1$，且$\beta_r$可表达为：
  

  $$\begin{aligned} \beta_0 &=& \frac{13\, {\left(v_{i} - 2\, v_{i-1} + v_{i-2}\right)}^2}{12} + \frac{{\left(3\, v_{i} - 4\, v_{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} - 2\, v_{i} + v_{i+1}\right)}^2}{12} \\ \beta_2 &=&\frac{13\, {\left(v_{i} - 2\, v_{i+1} + v_{i+2}\right)}^2}{12} + \frac{{\left(3\, v_{i} - 4\, v_{i+1} + v_{i+2}\right)}^2}{4} \end{aligned}$$

• 更新最后的左值为各个值的加权平均：即$u_{i+1/2}^{+} = \hat v_{i} = \sum\limits_{r=0}^{k-1}{w_r v_i^{(r)}}$。

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

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

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

$$f_{j+1/2} = u_{j+1/2}^{+} + u_{j+1/2}^{-}.$$

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

$$L(u_i) = - df/dx = - (f_{i+1/2} - f_{i-1/2})/\triangle x$$

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

3阶WENO格式

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

$$S_r(i) = \{x_{i-r},\dots,x_{i-r+k-1}\},r = 0,\dots,k-1.$$
其上重构值的系数（$u_{i+1/2}^{(r)} = \sum\limits_{j=0}^{k-1}{C_{rj}u_{i-r+j}}$），有人总结成了如下表格（右值）。

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

4阶ENO格式

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

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

计算结果展示和讨论

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

五阶WENO格式

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

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

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

三阶WENO格式

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

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

四阶ENO格式

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

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

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

五阶WENO格式

三阶WENO格式

四阶ENO格式

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

程序源码和其他话

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

一些小体会：

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

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

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

• 对于均匀网格，使用重构系数表$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.

---