一维激波管（Lax shock tube）问题的数值求解

一维激波管（Lax shock tube）问题的数值求解

问题描述

请分别尝试 Lax-Friedrichs, Lax-Wendroff,Roe和ENO格式求解一维激波管问题，并比较他们的数值结果。可以参考文献¹。

问题Lax shock tube定义如下：

$${\left(\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right)}_t+{\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ u(E+p)\end{array}\right)}_x=0,E=\frac{p}{\gamma -1}+\frac{1}{2}\rho u^2(\gamma =1.4)$$

with the initial condition:

$$\left(\begin{array}{c}{\rho }_L\\ u_L\\ p_L\end{array}\right)=\left(\begin{array}{c}0.445\\ 0.698\\ 3.528\end{array}\right)x<0.5,\left(\begin{array}{c}{\rho }_R\\ u_R\\ p_R\end{array}\right)=\left(\begin{array}{c}0.5\\ 0.0\\ 0.571\end{array}\right)x>0.5.$$

所谓的sod激波管问题，好像只是初值条件不一样？

一般的守恒算法格式

精确解

一维Riemann问题具有精确解.
因此可以用精确解来校验本文后面的数值算法的正确性和精确度。
对于本文的具体问题

• 初始条件: $\rho (x<0.5)=0.445$, $u(x<0.5)=0.698$,
  $p(x<0.5)=3.528$; $\rho (x>0.5)=0.5$, $u(x>0.5)=0$,
  $p(x>0.5)=0.571$。

• 边界条件: 在$x=0$和$x=1$处为反射边界条件, $u(x=0)=u(x=1)=0$。

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

精确解怎么求的，我不大清楚，但我觉得跟标量的时候应该是类似的，大概就是通过迎风线的方法通过求非线性函数零点追踪到初始值的点，然后用这个点值作为当前时刻待求点的值。网上找了一段求精确解的代码，跑了跑。得到结果如下图：

常用守恒格式（先介绍标量）

常用的守恒格式的写法为： $$u_j^{n+1}=u_j^n-\lambda ({\hat{f}}_{j+\frac{1}{2}}^n-{\hat{f}}_{j-\frac{1}{2}}^n)$$
这里的$\lambda =\frac{△t}{△x}$，$△t$和$△x$分别为时间步长和空间步长。
下面提到的几种格式，也是基于这样的一个守恒格式写法，只是数值通量函数$\hat{f}$
有所不同。

Lax-Friedrichs格式

Lax-Friedrichs格式的数值通量函数表达式为：
$$f_{j+\frac{1}{2}}^n=\frac{1}{2}[f(u_j^n)+f(u_{j+1}^n)-{\alpha }_n(u_{j+1}^n-u_j^n)]$$
其中， $${\alpha }_n=\underset{{\mathrm{u}}^{\mathrm{n}}}{\max }|f^{\prime }(u)|$$

Roe格式

Roe（迎风）格式的数值通量函数为（为了方便省略了第$n$时间层上标，下如是）：
$${\hat{f}}_{j+\frac{1}{2}}=\{\begin{array}{rlr}& f(u_j)& \frac{f(u_{j+1})-f(u_j)}{u_{j+1}-u_j}\ge 0\\ & f(u_{j+1})& \frac{f(u_{j+1})-f(u_j)}{u_{j+1}-u_j}<0\end{array}$$

Lax-Wendroff格式

Lax-Wendroff格式的数值通量表达式为：

$${\hat{f}}_{j+\frac{1}{2}}=\frac{1}{2}[f(u_j)+f(u_{j+1})-\lambda f^{\prime }(u_{j+\frac{1}{2}})(f(u_{j+1})-f(u_j))]\text{\hspace{1em}(LaxW)}$$
这里的$\lambda$含义同上，$u_{j+\frac{1}{2}}$的值可以通过插值得到，我选的是线性插值，即，$u_{j+\frac{1}{2}}=\frac{1}{2}(u_j+u_{j+1})$。

ENO格式简介（先介绍标量）

使用有限体积法4阶ENO格式，时间方向用三阶TVDRunge-Kutta方法。ENO格式是WENO格式的一个特例，我觉得可以先介绍WENO格式。

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^{\prime }(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{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$）：
  $$\begin{array}{rl}{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格式，有：$d_0=0.3,d_1=0.6,d_2=0.1$，且${\beta }_r$可表达为：
  $$\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}^{+}={\hat{v}}_i=\sum _{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})/△x$$

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

$$\begin{array}{rl}{u}^{(1)}& =u^n+\Delta tL\left(u^n\right)\\ u^{(2)}& =\frac{3}{4}{u}^n+\frac{1}{4}{u}^{(1)}+\frac{1}{4}\Delta tL\left(u^{(1)}\right)\\ u^{n+1}& =\frac{1}{3}{u}^n+\frac{2}{3}{u}^{(2)}+\frac{2}{3}\Delta tL\left(u^{(2)}\right)\end{array}$$

3阶WENO格式

三阶的WENO格式（k=2）和五阶的WENO格式（k=3）操作方法是类似的，所不同的是，在重构通量左右值时所用的模板和系数是不同的。一般来说，对于$2k-1$阶的WENO格式，所用的模板（具有对称性），共有$k$ 个：
S r ( i ) = { x i − r , … , x i − r + k − 1 } , r = 0 , … , k − 1. S_r(i) = \{x_{i-r},\dots,x_{i-r+k-1}\},r = 0,\dots,k-1. Sr​(i)={xi−r​,…,xi−r+k−1​},r=0,…,k−1.
其上重构值的系数（$u_{i+1/2}^{(r)}=\sum _{j=0}^{k-1}{C}_{rj}{u}_{i-r+j}$），有人总结成了如下表格（右值）。

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

$$\begin{array}{l}{d}_0=1,k=1\\ d_0=\frac{2}{3},d_1=\frac{1}{3},k=2\\ d_0=\frac{3}{10},d_1=\frac{3}{5},d_2=\frac{1}{10},k=3\end{array}$$

4阶ENO格式

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

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

从标量格式到系统的推广

我的理解是，不管是Lax-Friedrichs, Lax-Wendroff,Roe或者是ENO，格式和标量的时候是相同的，只不过是将$u$和$f$看成了三维的向量$\mathbf{u}$和$\mathbf{f}$来处理。不同的是，标量的时候，对通量函数求导，得到的是一个函数，而矢量的时候，${\mathbf{f}}^{\prime }(\mathbf{u})$是一个矩阵，在L-F方法中，通过查看书我发现，这个值发生了本质上的改变，下面会分别介绍。

常用的流通矢量分裂方法有Lax-Friedrichs（LF）分裂，Steger-Warming（SW）分裂和Van-Leer分裂等。

上面介绍的一些格式都是标量形式的，我们现在需要把它推广到矢量通量的形式。
上述问题可用压缩无黏气体Euler方程组来描述。对应于lax激波管问题在在直角坐标系下无量纲的一维Euler方程组为:

$$\frac{\partial \mathbf{w}}{\partial t}+\frac{\partial \mathbf{f}}{\partial x}=0,0\le x\le 1.0\text{\hspace{1em}(Eulerequation)}$$
其中 $$\mathbf{w}=\left[\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right],\mathbf{f}=\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ (E+p)u\end{array}\right]$$ 这里$\rho$是流体的密度, $u$是流体的速度, $p$是流体的压力, $E$是流体单位体积总能:
$$E=\frac{p}{\gamma -1}+\frac{1}{2}\rho u^2(\gamma =1.4)$$
$\mathbf{f}$可以写成$\mathbf{w}$为自变量的形式：
$$\mathbf{f}=\left(\begin{array}{c}\rho \mathrm{u}\\ \frac{2E}{5}+\frac{4{\rho \mathrm{u}}^2}{5\rho }\\ -\frac{{\rho \mathrm{u}}^3-7E\rho \rho \mathrm{u}}{5{\rho }^2}\end{array}\right)$$
式[Eulerequation]的Jacobian系数矩阵为
$$A(\mathbf{u}):={\mathbf{f}}^{\prime }(\mathbf{u})=\left(\begin{array}{ccc}0& 1& 0\\ -\frac{4}{5}{\left(\frac{\rho u}{\rho }\right)}^2& \frac{8}{5}\left(\frac{\rho u}{\rho }\right)& \frac{2}{5}\\ -\frac{7}{5}\left(\frac{\rho uE}{{\rho }^2}\right)+\frac{2}{5}{\left(\frac{\rho u}{\rho }\right)}^3& \frac{5}{7}\left(\frac{E}{\rho }\right)-\frac{3}{5}{\left(\frac{\rho u}{\rho }\right)}^2& \frac{7}{5}\left(\frac{\rho u}{\rho }\right)\end{array}\right)$$
这里的$\rho u$表示一个独立的变量而不是表示$\rho \ast u$。矩阵$\mathbf{A}$的特征值敲起来比较复杂，就不写了。

对应的特征向量为$e_1,e_2,e_3$。 对于 $$\frac{\partial \mathbf{w}}{\partial t}+A\frac{\partial \mathbf{w}}{\partial x}=0$$
由$A=R\Lambda R^{-1}$（其中，$R=(e_1,e_2,e_3)$，$\Lambda =\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\lambda }_1,{\lambda }_2,{\lambda }_3)$），可以得到表达：
$${\mathbf{v}}_t+\Lambda {\mathbf{v}}_x=0$$
其中的${\mathbf{v}}_t=R^{-1}\ast {\mathbf{w}}_t$。以上是我自己的一些思考，通过线性化，对系统进行解耦后再进行处理。但是这个方法并不可行。因为求起来太麻烦了。

对于时间上的离散，为了方便，我们采用显式欧拉格式。

Lax-Friedrichs格式

矢量形式的Lax-Friedrichs格式很简单，其实就是把标量形式中的表达换成向量的表达，再把${\alpha }_n$换成$\Delta x/\Delta t$即可，为什么要这样换，我暂时也不是很明白，我在书中看到的表达式是这个。故，向量形式的Lax-Friedrichs格式就可以表示为：
$$\{\begin{array}{l}{\mathbf{u}}_j^{n+1}={\mathbf{u}}_j^n-\frac{\Delta t}{\Delta x}\left({\hat{\mathbf{f}}}_{j+\frac{1}{2}}^n-{\hat{\mathbf{f}}}_{j-\frac{1}{2}}^n\right)\\ {\hat{\mathbf{f}}}_{j+\frac{1}{2}}^n=\frac{1}{2}\left[\mathbf{f}\left({\mathbf{u}}_j^n\right)+\mathbf{f}\left({\mathbf{u}}_{j+1}^n\right)-\frac{\Delta x}{\Delta t}\left({\mathbf{u}}_{j+1}^n-{\mathbf{u}}_j^n\right)\right]\end{array}$$
把这俩式子合并一下，其实也可以写成：
$${\mathbf{u}}_j^{n+1}=\frac{1}{2}\left({\mathbf{u}}_{j-1}^n+{\mathbf{u}}_{j+1}^n\right)-\frac{\Delta t}{2△x}\left({\mathbf{f}}_{j+1}^n-{\mathbf{f}}_{j-1}^m\right)$$

Lax-Wendroff格式

同样地，Lax-Wendroff格式，可以直接从标量改到向量，只不过是把通量的一个导数值变成了一个雅克比矩阵。
$$\{\begin{array}{rl}{\mathbf{u}}_j^{n+1}& ={\mathbf{u}}_j^n-\frac{\Delta t}{\Delta x}\left({\hat{\mathbf{f}}}_{j+\frac{1}{2}}^n-{\hat{\mathbf{f}}}_{j-\frac{1}{2}}^n\right)\\ {\hat{\mathbf{f}}}_{j+\frac{1}{2}}^n& =\frac{1}{2}\left[\mathbf{f}\left({\mathbf{u}}_j^n\right)+\mathbf{f}\left({\mathbf{u}}_{j+1}^n\right)-\frac{\Delta t}{\Delta x}{\mathbf{f}}^{\prime }\left({\mathbf{u}}_{j+\frac{1}{2}}^n\right)\left(\mathbf{f}\left({\mathbf{u}}_{j+1}^n\right)-\mathbf{f}\left({\mathbf{u}}_j^n\right)\right)\right]\\ {\mathbf{u}}_{j+\frac{1}{2}}^n& =\frac{{\mathbf{u}}_j^n+{\mathbf{u}}_{j+1}^n}{2}\end{array}$$
其中，雅克比矩阵计算得到的结果同上，即：
$${\mathbf{f}}^{\prime }(\mathbf{u})=\left(\begin{array}{ccc}0& 1& 0\\ -\frac{4}{5}{\left(\frac{\rho u}{\rho }\right)}^2& \frac{8}{5}\left(\frac{\rho u}{\rho }\right)& \frac{2}{5}\\ -\frac{7}{5}\left(\frac{\rho uE}{{\rho }^2}\right)+\frac{2}{5}{\left(\frac{\rho u}{\rho }\right)}^3& \frac{5}{7}\left(\frac{E}{\rho }\right)-\frac{3}{5}{\left(\frac{\rho u}{\rho }\right)}^2& \frac{7}{5}\left(\frac{\rho u}{\rho }\right)\end{array}\right)$$

Roe格式

Roe格式稍稍有些复杂。考虑：
$$\frac{\partial \mathbf{w}}{\partial t}+\frac{\partial \mathbf{f}}{\partial x}=0,0\le x\le 1.0$$
其中：
$$\mathbf{w}=\left[\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right],\mathbf{f}=\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ (E+p)u\end{array}\right]$$
这里的$\rho$是流体的密度，$u$是流体的速度，$p$是流体的压力，$E$是单位体积总能：
$$E=\rho \left(e+\frac{1}{2}{u}^2\right)$$
其中$e$为比内能（单位质量物质的内能），对于理想气体，有：
$$e=\frac{p}{(\gamma -1)\rho }⟹p=(\gamma -1)\rho e=(\gamma -1)\left[E-\frac{1}{2}\rho u^2\right]$$
对于理想气体，式中雅克比系数矩阵为：
$$\mathbf{A}=\left[\begin{array}{ccc}0& 1& 0\\ \frac{1}{2}(\gamma -3)u^2& (3-\gamma )u& \gamma -1\\ u\left(\frac{1}{2}(\gamma -1)u^2-H\right)& H-(\gamma -1)u^2& \gamma u\end{array}\right]$$
其中$H$为总比焓：
$$H=\frac{H+p}{\rho }=\frac{c^2}{\gamma -1}+\frac{1}{2}{u}^2$$ 声速为：
$$c=\sqrt{\frac{\gamma p}{\rho }}=\sqrt{\frac{\gamma }{\rho }(\gamma -1)\left(E-\frac{1}{2}\rho u^2\right)}$$
矩阵$A$的特征值为：
$${\lambda }_1=u-c,{\lambda }_2=u,{\lambda }_3=u+c$$
对应的特征向量为：
$${\mathbf{e}}_1=\left[\begin{array}{c}1\\ u-c\\ H-uc\end{array}\right],{\mathbf{e}}_2=\left[\begin{array}{c}1\\ u\\ \frac{1}{2}{u}^2\end{array}\right],{\mathbf{e}}_3=\left[\begin{array}{c}1\\ u+c\\ H+uc\end{array}\right]$$
这样我们就得到了一个Roe平均的计算方式：
$$\begin{array}{rl}\overline{\rho }& ={\left[\left(\sqrt{{\rho }_L}+\sqrt{{\rho }_R}\right)/2\right]}^2\\ \overline{u}& =\left(\sqrt{{\rho }_L}{u}_L+\sqrt{{\rho }_R}{u}_R\right)/\left(\sqrt{{\rho }_L}+\sqrt{{\rho }_R}\right)\\ \overline{H}& =\left(\sqrt{{\rho }_L}{H}_L+\sqrt{{\rho }_R}{H}_R\right)/\left(\sqrt{{\rho }_L}+\sqrt{{\rho }_R}\right)\end{array}$$
有了Roe平均，下面的过程就比较明朗了：
寻找$\tilde{A}\left({\mathbf{u}}_l,{\mathbf{u}}_r\right)$来替代雅克比矩阵$A(\mathbf{u})$，是$\tilde{A}$满足一定的性质，这里头上的弯弯，表示的就是Roe平均。
得到数值格式如下：
$$\{\begin{array}{l}{\mathbf{u}}_j^{n+1}={\mathbf{u}}_j^n-\frac{\Delta t}{\Delta x}\left({\hat{\mathbf{f}}}_{j+\frac{1}{2}}^n-{\hat{\mathbf{f}}}_{j-\frac{1}{2}}^n\right)\\ {\hat{\mathbf{f}}}_{j+\frac{1}{2}}^n=\frac{1}{2}\left[\mathbf{f}\left({\mathbf{u}}_j^n\right)+\mathbf{f}\left({\mathbf{u}}_{j+1}^n\right)-{\sum }_{n=1}^3\left|{\tilde{\lambda }}_k\right|{\tilde{\alpha }}_k{\mathbf{e}}_k\right]\end{array}$$

这里的${\tilde{\lambda }}_k$和${\tilde{e}}_k$分别表示近似矩阵$\tilde{A}$的特征值和特征向量，${\tilde{\alpha }}_k$
来自于间断于特征向量上的分解，即：
$${\mathbf{u}}_r-{\mathbf{u}}_l=\sum _{k=1}^3{\tilde{\alpha }}_k{\tilde{e}}_k$$
考虑黎曼问题的流通量表示：
$$\mathrm{f}(\mathbf{u})=\frac{1}{2}\left[\mathrm{f}\left({\mathbf{u}}_l\right)+\mathrm{f}\left({\mathbf{u}}_r\right)-|\tilde{A}|\left({\mathbf{u}}_r-{\mathbf{u}}_l\right)\right]$$
把间断的特征分解以及考虑$\tilde{A}$的特征值和特征向量，即可得：
$$\hat{\mathbf{f}}=\frac{1}{2}\left[\mathrm{f}\left({\mathbf{u}}_l\right)+\mathbf{f}\left({\mathbf{u}}_r\right)-\sum _{k=1}^3\left|{\tilde{\lambda }}_k\right|{\tilde{\alpha }}_k{\tilde{\mathbf{e}}}_k\right]$$

Eno格式

矢量Eno格式本质上可以直接从标量的形式推广过来。不同的在于，我们将原来的$\alpha =\max \left(\left|f^{\prime }\left(u_i\right)\right|\right)$换成了最大特征值。不再详细描述。

事实上，从标量到矢量，以上的格式除了Roe格式发生了一个质的变化之外，其他格式几乎都是直接从标量格式自然过渡过来加一些些微调，比如把通量导数变成网格比。

数值实验与结果

对以上提到的各个格式，分别做以下工作：1、计算各个格式在t=0.1时的误差和收敛阶。2、画出t=0.1时的数值解的图，每种格式单独画一个图。精确解可以参考最前面那个图，就不再画到图上去了，显得有些乱。

Lax-Friedrichs格式

计算实践表明，由于Lax差分格式的黏性较大，计算精度不高，在计算激波时，激波会被拉宽和抹平，一般都不采用Lax差分格式计算激波。它在时间和空间上都是一阶精度，Lax差分格式是一个非线性守恒型差分格式。对它进行线性化后,可得到它的稳定性条件$|\mathbf{w}{|}_{\max }\Delta t/\Delta x\le 1$。

这个格式，表现出来的特征就是耗散得比较厉害，粘性过大，导致解被磨平了，体现不出来间断的特征（精确解其实就是间断的传播）。

Lax-Wendroff格式

从表中可以看出，Lax-Wendroff格式对间断而言可能是一阶收敛的。事实上，刚开始时收敛阶可能表现得不那么稳定，但是随着步长的缩短，这个格式收敛阶最终是稳定在一阶的。这个格式在间断点处震荡得比较厉害。

Roe格式

Roe格式的收敛阶和图如下所示：


Roe格式总是在拐点处抖得比较厉害，拐得越厉害，抖得也就越厉害，在间断处抖动震荡是最为厉害的。

ENO格式

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

以上四种格式在解间断时只能达到一阶精度。

ENO格式的图是目前看起来最顺眼和舒服的一个图，虽然这几个格式针对这个间断问题都是一阶收敛，但是ENO格式对于局部的刻画，特别是间断处的刻画是最好的。

程序源码和编程小体会

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

一点小感悟：

• 一阶收敛的数值格式是很容易实现的，即使在编程过程中，那个小地方写错了，也可能会达到一阶。也就是说，一阶格式有一阶收敛程序不一定是对的，但是没有一阶收敛，那么肯定是错的。这是一件神奇的事情。

• 试图将这些程序写成通用的解决一阶标量双曲守恒方程的形式，也就是说只是通量函数
  $f$
  改变，只要在fun函数脚本中修改通量函数，那么这些格式依然能用。有一个小问题在于，写成通用的，多了一些求导等不要的计算，加大了运算量。在ENO
  格式中，没有写成通用格式。主要考虑到ENO格式中积分的计算，若对于一般的通量函数$f$，正负部截断函数的原函数不一定能求，那么这是可能要用到一些数值积分手段来求近似积分，这对一维黎曼方程是多余的，引入了额外误差。

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

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

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

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

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

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

• 一个高阶收敛的格式，刚开始的时候收敛阶可能表现得不大稳定，比如上面提到Lax-Wendroff
  格式，它的收敛阶一会高，一会儿低，但是随着网格的加密，它往往会稳定到它真正的收敛阶。

参考文献

^12345678

---

1. Lax P D . Weak solutions of nonlinear hyperbolic equations and their numerical computation[J]. Comm. Pure Appl. Math. 1954, 7(1):159-193. ↩︎ ↩︎

2. Shu C W. Numerical methods for hyperbolic conservation laws[J].
   Lecture notes (AM257), 2006. ↩︎

3. 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. ↩︎

4. Jiang G S, Shu C W. Efficient implementation of weighted ENO schemes[J]. Journal of computational physics, 1996, 126(1): 202-228. ↩︎

5. Anderson J D, Wendt J Computational fluid dynamics[M]. New York: McGraw-Hill, 1995. ↩︎

6. Landau L D, Lifshitz E M. Fluid Mechanics: Volume 6 (Pergamon International
   Library of Science, Technology, Engineering & Social Studies)[J].2013. ↩︎

7. 庄礼贤,严协远,流体力学,中国科技大学出版社,1992 ↩︎

8. LeVeque R J. Finite difference methods for ordinary and partial differential equations:steady-state and time-dependent problems[M]. Siam, 2007. ↩︎