请分别尝试 Lax-Friedrichs, Lax-Wendroff,Roe和ENO格式求解一维激波管问题，并比较他们的数值结果。可以参考文献1
问题Lax shock tube定义如下：
( ρ ρ u E ) t + ( ρ u ρ u 2 + p u ( E + p ) ) x = 0 , E = p γ − 1 + 1 2 ρ u 2 ( γ = 1.4 ) \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,
\quad E=\frac{p}{\gamma -1} + \frac{1}{2}\rho u^2 (\gamma = 1.4) ⎝ ⎛ ρ ρ u E ⎠ ⎞ t + ⎝ ⎛ ρ u ρ u 2 + p u ( E + p ) ⎠ ⎞ x = 0 , E = γ − 1 p + 2 1 ρ u 2 ( γ = 1 . 4 )
with the initial condition:
( ρ L u L p L ) = ( 0.445 0.698 3.528 ) x < 0.5 , ( ρ R u R p R ) = ( 0.5 0.0 0.571 ) x > 0.5. \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, \quad
\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. ⎝ ⎛ ρ L u L p L ⎠ ⎞ = ⎝ ⎛ 0 . 4 4 5 0 . 6 9 8 3 . 5 2 8 ⎠ ⎞ x < 0 . 5 , ⎝ ⎛ ρ R u R p R ⎠ ⎞ = ⎝ ⎛ 0 . 5 0 . 0 0 . 5 7 1 ⎠ ⎞ x > 0 . 5 .
所谓的sod激波管问题，好像只是初值条件不一样？
一维Riemann问题具有精确解.
初始条件: ρ ( x < 0.5 ) = 0.445 \rho(x<0.5) = 0.445 ρ ( x < 0 . 5 ) = 0 . 4 4 5 u ( x < 0.5 ) = 0.698 u(x<0.5) = 0.698 u ( x < 0 . 5 ) = 0 . 6 9 8 p ( x < 0.5 ) = 3.528 p(x<0.5) = 3.528 p ( x < 0 . 5 ) = 3 . 5 2 8 ρ ( x > 0.5 ) = 0.5 \rho(x>0.5) = 0.5 ρ ( x > 0 . 5 ) = 0 . 5 u ( x > 0.5 ) = 0 u(x>0.5) = 0 u ( x > 0 . 5 ) = 0 p ( x > 0.5 ) = 0.571 p(x>0.5) = 0.571 p ( x > 0 . 5 ) = 0 . 5 7 1
边界条件: 在x = 0 x=0 x = 0 x = 1 x=1 x = 1 u ( x = 0 ) = u ( x = 1 ) = 0 u(x=0)=u(x=1)=0 u ( x = 0 ) = u ( x = 1 ) = 0
精确解也可以采用5阶WENO格式的小步长逼近作为相对近似精确解，这也只是相对而言的精确解。将精确解保存为
精确解怎么求的，我不大清楚，但我觉得跟标量的时候应该是类似的，大概就是通过迎风线的方法通过求非线性函数零点追踪到初始值的点，然后用这个点值作为当前时刻待求点的值。网上找了一段求精确解的代码，跑了跑。得到结果如下图：
常用的守恒格式的写法为： u j n + 1 = u j n − λ ( f ^ j + 1 2 n − f ^ j − 1 2 n )
u_j^{n+1} = u_j^n - \lambda (\hat {f}_{j+\frac{1}{2}}^n-\hat {f}_{j-\frac{1}{2}}^n) u j n + 1 = u j n − λ ( f ^ j + 2 1 n − f ^ j − 2 1 n ) λ = △ t △ x \lambda = \frac{\triangle t}{\triangle x} λ = △ x △ t △ t \triangle t △ t △ x \triangle x △ x f ^ \hat f f ^
Lax-Friedrichs格式的数值通量函数表达式为：f j + 1 2 n = 1 2 [ f ( u j n ) + f ( u j + 1 n ) − α n ( u j + 1 n − u j n ) ] 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)] f j + 2 1 n = 2 1 [ f ( u j n ) + f ( u j + 1 n ) − α n ( u j + 1 n − u j n ) ] α n = max u n ∣ f ′ ( u ) ∣
\alpha_n = \mathop {\max }\limits_{\rm{u^n}} |f'(u)| α n = u n max ∣ f ′ ( u ) ∣
Roe（迎风）格式的数值通量函数为（为了方便省略了第n n n f ^ j + 1 2 = { f ( u j ) f ( u j + 1 ) − f ( u j ) u j + 1 − u j ≥ 0 f ( u j + 1 ) f ( u j + 1 ) − f ( u j ) u j + 1 − u j < 0
\hat f _{j+\frac{1}{2}} =
\left \{
\begin{aligned}
& f(u_j) & \frac{f(u_{j+1})-f(u_{j})}{u_{j+1}-u_j} \geq 0\\
& f(u_{j+1}) & \frac{f(u_{j+1})-f(u_{j})}{u_{j+1}-u_j} < 0
\end{aligned}
\right. f ^ j + 2 1 = ⎩ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎧ f ( u j ) f ( u j + 1 ) u j + 1 − u j f ( u j + 1 ) − f ( u j ) ≥ 0 u j + 1 − u j f ( u j + 1 ) − f ( u j ) < 0
Lax-Wendroff格式的数值通量表达式为：
(LaxW) f ^ j + 1 2 = 1 2 [ f ( u j ) + f ( u j + 1 ) − λ f ′ ( u j + 1 2 ) ( f ( u j + 1 ) − f ( u j ) ) ] \tag{LaxW}
\hat f _{j+\frac{1}{2}} = \frac{1}{2}[f(u_j)+f(u_{j+1})-\lambda f'(u_{j+\frac{1}{2}})(f(u_{j+1})-f(u_{j}))] f ^ j + 2 1 = 2 1 [ f ( u j ) + f ( u j + 1 ) − λ f ′ ( u j + 2 1 ) ( f ( u j + 1 ) − f ( u j ) ) ] ( L a x W ) λ \lambda λ u j + 1 2 u_{j+\frac{1}{2}} u j + 2 1 u j + 1 2 = 1 2 ( u j + u j + 1 ) u_{j+\frac{1}{2}} = \frac{1}{2}(u_j+u_{j+1}) u j + 2 1 = 2 1 ( u j + u j + 1 )
使用有限体积法4阶ENO格式，时间方向用三阶TVDRunge-Kutta方法。ENO格式是WENO格式的一个特例，我觉得可以先介绍WENO格式。
空间上采用五阶WENO格式，时间上采用三阶Runge Kutta格式，具体操作描述如下：
将空间定义域等分，取每个等分区间上中点的初值作为该区间的初值（不是计算等分节点的值）。那么，N等分就有N个值，对这N个值编号，不妨记为u i u_i u i
使用Lax-Friedrichs，进行通量分裂，如下：f ^ j + 1 / 2 = 1 2 ( f ( u j ) + α u j ) + 1 2 ( f ( u j + 1 ) − α u j + 1 ) \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}) f ^ j + 1 / 2 = 2 1 ( f ( u j ) + α u j ) + 2 1 ( f ( u j + 1 ) − α u j + 1 ) α = max ( ∣ f ′ ( u i ) ∣ ) \alpha = \text{max}(|f'(u_i)|) α = max ( ∣ f ′ ( u i ) ∣ ) u i u_i u i u j + 1 / 2 + = 1 2 ( f ( u j ) + α u j ) u^+_{j+1/2} = \frac{1}{2}(f(u_j)+\alpha u_j) u j + 1 / 2 + = 2 1 ( f ( u j ) + α u j ) u j + 1 / 2 − = 1 2 ( f ( u j + 1 ) − α u j + 1 ) u^-_{j+1/2} = \frac{1}{2}(f(u_{j+1})-\alpha u_{j+1}) u j + 1 / 2 − = 2 1 ( f ( u j + 1 ) − α u j + 1 )
下面以右值u j + 1 / 2 − u^-_{j+1/2} u j + 1 / 2 − v i v_{i} v i u j + 1 / 2 − u^-_{j+1/2} u j + 1 / 2 −
选择不同的模板，计算不同模板下的重构值，采用如图所示的模板，计算得的对应模板的三个重构值如下：
v i ( 0 ) = ( 2 v i − 2 − 7 v i − 1 + 11 v i ) / 6 v i ( 1 ) = ( − v i − 1 + 5 v i + 2 v i + 1 ) / 6 v i ( 2 ) = ( 2 v i + 5 v i + 1 − v i + 2 ) / 6 \begin{aligned}
% \nonumber to remove numbering (before each equation)
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} v i ( 0 ) v i ( 1 ) v i ( 2 ) = = = ( 2 v i − 2 − 7 v i − 1 + 1 1 v i ) / 6 ( − v i − 1 + 5 v i + 2 v i + 1 ) / 6 ( 2 v i + 5 v i + 1 − v i + 2 ) / 6
计算三个模板得到的重构值对应的权重，逆向罗列公式如下（k = 3 k=3 k = 3 w r = α r ∑ r = 0 k − 1 α r α r = d r ( 1 0 − 6 + β r ) 2 \begin{aligned}
% \nonumber to remove numbering (before each equation)
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} w r α r = r = 0 ∑ k − 1 α r α r = ( 1 0 − 6 + β r ) 2 d r
对于五阶WENO格式，有：d 0 = 0.3 , d 1 = 0.6 , d 2 = 0.1 d_0 = 0.3,d_1 = 0.6,d_2 = 0.1 d 0 = 0 . 3 , d 1 = 0 . 6 , d 2 = 0 . 1 β r \beta_r β r β 0 = 13  ( v i − 2  v i − 1 + v i − 2 ) 2 12 + ( 3  v i − 4  v i − 1 + v i − 2 ) 2 4 β 1 = ( v i − 1 − v i + 1 ) 2 4 + 13  ( v i − 1 − 2  v i + v i + 1 ) 2 12 β 2 = 13  ( v i − 2  v i + 1 + v i + 2 ) 2 12 + ( 3  v i − 4  v i + 1 + v i + 2 ) 2 4 \begin{aligned}
% \nonumber to remove numbering (before each equation)
\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} β 0 β 1 β 2 = = = 1 2 1 3 ( v i − 2 v i − 1 + v i − 2 ) 2 + 4 ( 3 v i − 4 v i − 1 + v i − 2 ) 2 4 ( v i − 1 − v i + 1 ) 2 + 1 2 1 3 ( v i − 1 − 2 v i + v i + 1 ) 2 1 2 1 3 ( v i − 2 v i + 1 + v i + 2 ) 2 + 4 ( 3 v i − 4 v i + 1 + v i + 2 ) 2
更新最后的左值为各个值的加权平均：即u i + 1 / 2 + = v ^ i = ∑ r = 0 k − 1 w r v i ( r ) u_{i+1/2}^{+} = \hat v_{i} = \sum\limits_{r=0}^{k-1}{w_r v_i^{(r)}} u i + 1 / 2 + = v ^ i = r = 0 ∑ k − 1 w r v i ( r )
同相同的方法重构左值。原来从左到右计算的，重构左值时从右到左。
通过如上方法，计算得到更新后的左值和右值，取左值和右值的和为通量函数近似，即：f j + 1 / 2 = u j + 1 / 2 + + u j + 1 / 2 − . f_{j+1/2} = u_{j+1/2}^{+} + u_{j+1/2}^{-}. f j + 1 / 2 = u j + 1 / 2 + + u j + 1 / 2 − .
那么，可以得到空间离散算子的作用结果，即：L ( u i ) = − d f / d x = − ( f i + 1 / 2 − f i − 1 / 2 ) / △ x L(u_i) = - df/dx = - (f_{i+1/2} - f_{i-1/2})/\triangle x L ( u i ) = − d f / d x = − ( f i + 1 / 2 − f i − 1 / 2 ) / △ x
有了空间离散算子，我们就在时间上使用三阶Runge-Kutta方法，求解下一时间层的值。
u ( 1 ) = u n + Δ t L ( u n ) u ( 2 ) = 3 4 u n + 1 4 u ( 1 ) + 1 4 Δ t L ( u ( 1 ) ) u n + 1 = 1 3 u n + 2 3 u ( 2 ) + 2 3 Δ t L ( u ( 2 ) ) \begin{aligned} u^{(1)} &=u^{n}+\Delta t L\left(u^{n}\right) \\ u^{(2)} &=\frac{3}{4} u^{n}+\frac{1}{4} u^{(1)}+\frac{1}{4} \Delta t L\left(u^{(1)}\right) \\ u^{n+1} &=\frac{1}{3} u^{n}+\frac{2}{3} u^{(2)}+\frac{2}{3} \Delta t L\left(u^{(2)}\right) \end{aligned} u ( 1 ) u ( 2 ) u n + 1 = u n + Δ t L ( u n ) = 4 3 u n + 4 1 u ( 1 ) + 4 1 Δ t L ( u ( 1 ) ) = 3 1 u n + 3 2 u ( 2 ) + 3 2 Δ t L ( u ( 2 ) )
三阶的WENO格式（k=2）和五阶的WENO格式（k=3）操作方法是类似的，所不同的是，在重构通量左右值时所用的模板和系数是不同的。一般来说，对于2 k − 1 2k-1 2 k − 1 k k 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. S r ( i ) = { x i − r , … , x i − r + k − 1 } , r = 0 , … , k − 1 . u i + 1 / 2 ( r ) = ∑ j = 0 k − 1 C r j u i − r + j u_{i+1/2}^{(r)} = \sum\limits_{j=0}^{k-1}{C_{rj}u_{i-r+j}} u i + 1 / 2 ( r ) = j = 0 ∑ k − 1 C r j u i − r + j
且对于k比较小时，d r , β r d_r,\beta _r d r , β r
d 0 = 1 , k = 1 d 0 = 2 3 , d 1 = 1 3 , k = 2 d 0 = 3 10 , d 1 = 3 5 , d 2 = 1 10 , k = 3 \begin{array}{l}{d_{0}=1, \quad k=1} \\ {d_{0}=\frac{2}{3}, \quad d_{1}=\frac{1}{3}, \quad k=2} \\ {d_{0}=\frac{3}{10}, \quad d_{1}=\frac{3}{5}, \quad d_{2}=\frac{1}{10}, \quad k=3}\end{array} d 0 = 1 , k = 1 d 0 = 3 2 , d 1 = 3 1 , k = 2 d 0 = 1 0 3 , d 1 = 5 3 , d 2 = 1 0 1 , k = 3
ENO格式和WENO格式的过程大部分是类似的，不同的是，ENO格式左右值的重构，不需要计算全部模板下的重构值后取加权平均，而只要选取其中的某一个模板作为计算左右值重构即可，选取的准则为使差商最小。具体的算法描述如下。
模板的重构系数C r j C_{rj} C r j
我的理解是，不管是Lax-Friedrichs, Lax-Wendroff,Roe或者是ENO，格式和标量的时候是相同的，只不过是将u u u f f f u \mathbf{u} u f \mathbf{f} f f ′ ( u ) \mathbf{f}'(\mathbf{u}) f ′ ( u )
常用的流通矢量分裂方法有Lax-Friedrichs（LF）分裂，Steger-Warming（SW）分裂和Van-Leer分裂等。
上面介绍的一些格式都是标量形式的，我们现在需要把它推广到矢量通量的形式。
(Eulerequation) ∂ w ∂ t + ∂ f ∂ x = 0 ,  0 ≤ x ≤ 1.0 \tag{Eulerequation}
\frac{\partial \mathbf{w}}{\partial t}+\frac{\partial \mathbf{f}}{\partial x} = 0,\: 0\leq x\leq 1.0 ∂ t ∂ w + ∂ x ∂ f = 0 , 0 ≤ x ≤ 1 . 0 ( E u l e r e q u a t i o n ) w = [ ρ ρ u E ] ,   f = [ ρ u ρ u 2 + p ( E + p ) u ] \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] w = ⎣ ⎡ ρ ρ u E ⎦ ⎤ , f = ⎣ ⎡ ρ u ρ u 2 + p ( E + p ) u ⎦ ⎤ ρ \rho ρ u u u p p p E E E E = p γ − 1 + 1 2 ρ u 2 ( γ = 1.4 ) E=\frac{p}{\gamma -1} + \frac{1}{2}\rho u^2 (\gamma = 1.4) E = γ − 1 p + 2 1 ρ u 2 ( γ = 1 . 4 ) f \mathbf{f} f w \mathbf{w} w f = ( ρ u 2  E 5 + 4  ρ u 2 5  ρ − ρ u 3 − 7  E  ρ  ρ u 5  ρ 2 ) \mathbf{f} = \left(\begin{array}{c} \mathrm{\rho u}\\ \frac{2\, E}{5} + \frac{4\, {\mathrm{\rho u}}^2}{5\, \mathrm{\rho}}\\ -\frac{{\mathrm{\rho u}}^3 - 7\, E\, \mathrm{\rho}\, \mathrm{\rho u}}{5\, {\mathrm{\rho}}^2} \end{array}\right) f = ⎝ ⎜ ⎛ ρ u 5 2 E + 5 ρ 4 ρ u 2 − 5 ρ 2 ρ u 3 − 7 E ρ ρ u ⎠ ⎟ ⎞ [Eulerequation] 的Jacobian系数矩阵为A ( u ) : = f ′ ( u ) = ( 0 1 0 − 4 5 ( ρ u ρ ) 2 8 5 ( ρ u ρ ) 2 5 − 7 5 ( ρ u E ρ 2 ) + 2 5 ( ρ u ρ ) 3 5 7 ( E ρ ) − 3 5 ( ρ u ρ ) 2 7 5 ( ρ u ρ ) ) 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 u E}{\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) A ( u ) : = f ′ ( u ) = ⎝ ⎜ ⎜ ⎛ 0 − 5 4 ( ρ ρ u ) 2 − 5 7 ( ρ 2 ρ u E ) + 5 2 ( ρ ρ u ) 3 1 5 8 ( ρ ρ u ) 7 5 ( ρ E ) − 5 3 ( ρ ρ u ) 2 0 5 2 5 7 ( ρ ρ u ) ⎠ ⎟ ⎟ ⎞ ρ u \rho u ρ u ρ ∗ u \rho*u ρ ∗ u A \mathbf{A} A
对应的特征向量为e 1 , e 2 , e 3 e_1,e_2,e_3 e 1 , e 2 , e 3 ∂ w ∂ t + A ∂ w ∂ x = 0
\frac{\partial \mathbf{w}}{\partial t}+A\frac{\partial \mathbf{w}}{\partial x} = 0 ∂ t ∂ w + A ∂ x ∂ w = 0 A = R Λ R − 1 A = R\Lambda R^{-1} A = R Λ R − 1 R = ( e 1 , e 2 , e 3 ) R = (e_1,e_2,e_3) R = ( e 1 , e 2 , e 3 ) Λ = d i a g ( λ 1 , λ 2 , λ 3 ) \Lambda = \mathrm{diag}(\lambda_1,\lambda_2,\lambda_3) Λ = d i a g ( λ 1 , λ 2 , λ 3 ) v t + Λ v x = 0 \mathbf{v}_t + \Lambda \mathbf{v}_x = 0 v t + Λ v x = 0 v t = R − 1 ∗ w t \mathbf{v}_t = R^{-1}*\mathbf{w}_t v t = R − 1 ∗ w t
对于时间上的离散，为了方便，我们采用显式欧拉格式。
矢量形式的Lax-Friedrichs格式很简单，其实就是把标量形式中的表达换成向量的表达，再把α n \alpha_n α n Δ x / Δ t \Delta x/\Delta t Δ x / Δ t { u j n + 1 = u j n − Δ t Δ x ( f ^ j + 1 2 n − f ^ j − 1 2 n ) f ^ j + 1 2 n = 1 2 [ f ( u j n ) + f ( u j + 1 n ) − Δ x Δ t ( u j + 1 n − u j n ) ] \left\{\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}\right. ⎩ ⎨ ⎧ u j n + 1 = u j n − Δ x Δ t ( f ^ j + 2 1 n − f ^ j − 2 1 n ) f ^ j + 2 1 n = 2 1 [ f ( u j n ) + f ( u j + 1 n ) − Δ t Δ x ( u j + 1 n − u j n ) ] u j n + 1 = 1 2 ( u j − 1 n + u j + 1 n ) − Δ t 2 △ x ( f j + 1 n − f j − 1 m ) \mathbf{u}_{j}^{n+1}=\frac{1}{2}\left(\mathbf{u}_{j-1}^{n}+\mathbf{u}_{j+1}^{n}\right)-\frac{\Delta t}{2 \triangle x}\left(\mathbf{f}_{j+1}^{n}-\mathbf{f}_{j-1}^{m}\right) u j n + 1 = 2 1 ( u j − 1 n + u j + 1 n ) − 2 △ x Δ t ( f j + 1 n − f j − 1 m )
同样地，Lax-Wendroff格式，可以直接从标量改到向量，只不过是把通量的一个导数值变成了一个雅克比矩阵。{ u j n + 1 = u j n − Δ t Δ x ( f ^ j + 1 2 n − f ^ j − 1 2 n ) f ^ j + 1 2 n = 1 2 [ f ( u j n ) + f ( u j + 1 n ) − Δ t Δ x f ′ ( u j + 1 2 n ) ( f ( u j + 1 n ) − f ( u j n ) ) ] u j + 1 2 n = u j n + u j + 1 n 2 \left\{\begin{aligned} \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{aligned}\right. ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎧ u j n + 1 f ^ j + 2 1 n u j + 2 1 n = u j n − Δ x Δ t ( f ^ j + 2 1 n − f ^ j − 2 1 n ) = 2 1 [ f ( u j n ) + f ( u j + 1 n ) − Δ x Δ t f ′ ( u j + 2 1 n ) ( f ( u j + 1 n ) − f ( u j n ) ) ] = 2 u j n + u j + 1 n f ′ ( u ) = ( 0 1 0 − 4 5 ( ρ u ρ ) 2 8 5 ( ρ u ρ ) 2 5 − 7 5 ( ρ u E ρ 2 ) + 2 5 ( ρ u ρ ) 3 5 7 ( E ρ ) − 3 5 ( ρ u ρ ) 2 7 5 ( ρ 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 u E}{\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) f ′ ( u ) = ⎝ ⎜ ⎜ ⎛ 0 − 5 4 ( ρ ρ u ) 2 − 5 7 ( ρ 2 ρ u E ) + 5 2 ( ρ ρ u ) 3 1 5 8 ( ρ ρ u ) 7 5 ( ρ E ) − 5 3 ( ρ ρ u ) 2 0 5 2 5 7 ( ρ ρ u ) ⎠ ⎟ ⎟ ⎞
Roe格式稍稍有些复杂。考虑：∂ w ∂ t + ∂ f ∂ x = 0 , 0 ≤ x ≤ 1.0 \frac{\partial \mathbf{w}}{\partial t}+\frac{\partial \mathbf{f}}{\partial x}=0,0 \leq x \leq 1.0 ∂ t ∂ w + ∂ x ∂ f = 0 , 0 ≤ x ≤ 1 . 0 w = [ ρ ρ u E ] , f = [ ρ u ρ u 2 + p ( E + p ) u ] \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] w = ⎣ ⎡ ρ ρ u E ⎦ ⎤ , f = ⎣ ⎡ ρ u ρ u 2 + p ( E + p ) u ⎦ ⎤ ρ \rho ρ u u u p p p E E E E = ρ ( e + 1 2 u 2 ) E=\rho\left(e+\frac{1}{2} u^{2}\right) E = ρ ( e + 2 1 u 2 ) e e e e = p ( γ − 1 ) ρ ⟹ p = ( γ − 1 ) ρ e = ( γ − 1 ) [ E − 1 2 ρ u 2 ] e=\frac{p}{(\gamma-1) \rho} \Longrightarrow p=(\gamma-1) \rho e=(\gamma-1)\left[E-\frac{1}{2} \rho u^{2}\right] e = ( γ − 1 ) ρ p ⟹ p = ( γ − 1 ) ρ e = ( γ − 1 ) [ E − 2 1 ρ u 2 ] A = [ 0 1 0 1 2 ( γ − 3 ) u 2 ( 3 − γ ) u γ − 1 u ( 1 2 ( γ − 1 ) u 2 − H ) H − ( γ − 1 ) u 2 γ u ] \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] A = ⎣ ⎡ 0 2 1 ( γ − 3 ) u 2 u ( 2 1 ( γ − 1 ) u 2 − H ) 1 ( 3 − γ ) u H − ( γ − 1 ) u 2 0 γ − 1 γ u ⎦ ⎤ H H H H = H + p ρ = c 2 γ − 1 + 1 2 u 2 H=\frac{H+p}{\rho}=\frac{c^{2}}{\gamma-1}+\frac{1}{2} u^{2} H = ρ H + p = γ − 1 c 2 + 2 1 u 2 c = γ p ρ = γ ρ ( γ − 1 ) ( E − 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)} c = ρ γ p = ρ γ ( γ − 1 ) ( E − 2 1 ρ u 2 ) A A A λ 1 = u − c , λ 2 = u , λ 3 = u + c \lambda_{1}=u-c, \quad \lambda_{2}=u, \quad \lambda_{3}=u+c λ 1 = u − c , λ 2 = u , λ 3 = u + c e 1 = [ 1 u − c H − u c ] , e 2 = [ 1 u 1 2 u 2 ] , e 3 = [ 1 u + c H + u c ] \mathbf{e}_{1}=\left[ \begin{array}{c}{1} \\ {u-c} \\ {H-u c}\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+u c}\end{array}\right] e 1 = ⎣ ⎡ 1 u − c H − u c ⎦ ⎤ , e 2 = ⎣ ⎡ 1 u 2 1 u 2 ⎦ ⎤ , e 3 = ⎣ ⎡ 1 u + c H + u c ⎦ ⎤ ρ ‾ = [ ( ρ L + ρ R ) / 2 ] 2 u ‾ = ( ρ L u L + ρ R u R ) / ( ρ L + ρ R ) H ‾ = ( ρ L H L + ρ R H R ) / ( ρ L + ρ R ) \begin{aligned} \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{aligned} ρ u H = [ ( ρ L + ρ R ) / 2 ] 2 = ( ρ L u L + ρ R u R ) / ( ρ L + ρ R ) = ( ρ L H L + ρ R H R ) / ( ρ L + ρ R ) A ~ ( u l , u r ) \tilde{A}\left(\mathbf{u}_{l}, \mathbf{u}_{r}\right) A ~ ( u l , u r ) A ( u ) A(\mathbf{u}) A ( u ) A ~ \tilde{A} A ~ { u j n + 1 = u j n − Δ t Δ x ( f ^ j + 1 2 n − f ^ j − 1 2 n ) f ^ j + 1 2 n = 1 2 [ f ( u j n ) + f ( u j + 1 n ) − ∑ n = 1 3 ∣ λ ~ k ∣ α ~ k e k ] \left\{\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}\right. ⎩ ⎨ ⎧ u j n + 1 = u j n − Δ x Δ t ( f ^ j + 2 1 n − f ^ j − 2 1 n ) f ^ j + 2 1 n = 2 1 [ f ( u j n ) + f ( u j + 1 n ) − ∑ n = 1 3 ∣ ∣ ∣ λ ~ k ∣ ∣ ∣ α ~ k e k ]
这里的λ ~ k \tilde \lambda_k λ ~ k e ~ k \tilde e_k e ~ k A ~ \tilde A A ~ α ~ k \tilde \alpha_k α ~ k u r − u l = ∑ k = 1 3 α ~ k e ~ k \mathbf{u}_{r}-\mathbf{u}_{l}=\sum_{k=1}^{3} \tilde{\alpha}_{k} \tilde{e}_{k} u r − u l = k = 1 ∑ 3 α ~ k e ~ k f ( u ) = 1 2 [ f ( u l ) + f ( u r ) − ∣ A ~ ∣ ( u r − u l ) ] \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] f ( u ) = 2 1 [ f ( u l ) + f ( u r ) − ∣ A ~ ∣ ( u r − u l ) ] A ~ \tilde A A ~ f ^ = 1 2 [ f ( u l ) + f ( u r ) − ∑ k = 1 3 ∣ λ ~ k ∣ α ~ k e ~ k ] \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] f ^ = 2 1 [ f ( u l ) + f ( u r ) − k = 1 ∑ 3 ∣ ∣ ∣ λ ~ k ∣ ∣ ∣ α ~ k e ~ k ]
矢量Eno格式本质上可以直接从标量的形式推广过来。不同的在于，我们将原来的α = max ( ∣ f ′ ( u i ) ∣ ) \alpha=\max \left(\left|f^{\prime}\left(u_{i}\right)\right|\right) α = max ( ∣ f ′ ( u i ) ∣ )
事实上，从标量到矢量，以上的格式除了Roe格式发生了一个质的变化之外，其他格式几乎都是直接从标量格式自然过渡过来加一些些微调，比如把通量导数变成网格比。
对以上提到的各个格式，分别做以下工作：1、计算各个格式在t=0.1时的误差和收敛阶。2、画出t=0.1时的数值解的图，每种格式单独画一个图。精确解可以参考最前面那个图，就不再画到图上去了，显得有些乱。
计算实践表明，由于Lax差分格式的黏性较大，计算精度不高，在计算激波时，激波会被拉宽和抹平，一般都不采用Lax差分格式计算激波。它在时间和空间上都是一阶精度，Lax差分格式是一个非线性守恒型差分格式。对它进行线性化后,可得到它的稳定性条件∣ w ∣ max Δ t / Δ x ≤ 1 |\mathbf{w}|_{\max } \Delta t / \Delta x \leq 1 ∣ w ∣ max Δ t / Δ x ≤ 1
这个格式，表现出来的特征就是耗散得比较厉害，粘性过大，导致解被磨平了，体现不出来间断的特征（精确解其实就是间断的传播）。
Roe格式的收敛阶和图如下所示：
四阶ENO格式，在解连续的情况下，得到的结果如下表：有些小问题。为了减小时间方向带来的影响，设置Courant系数尽量小。
源代码见附件，使用的语言为MATLAB2014a（盗版）,Windows 10环境。
一点小感悟：
一阶收敛的数值格式是很容易实现的，即使在编程过程中，那个小地方写错了，也可能会达到一阶。也就是说，一阶格式有一阶收敛程序不一定是对的，但是没有一阶收敛，那么肯定是错的。这是一件神奇的事情。
试图将这些程序写成通用的解决一阶标量双曲守恒方程的形式，也就是说只是通量函数f f f f f f
Matlab巧用shiftcirt函数，来处理周期边界问题是特别棒的事情。
知道如何进行右值重构后，再进行计算左值重构，我们没必要考虑系数的镜像对称关系。一个简单的做法是，将左值反向排列，然后同样用计算右值的程序计算一遍，得到的结果最后再转个方向即可。
CFL迎风条件： dt/dx * alpha = CFL
对于均匀网格，使用重构系数表C r j C_{rj} C r j
初值点的选取，选的不是等分节点处，而是选取单元中点处。使用单元中点值表示这个单元的值。
对这些格式，在有些步长下能得到相应的收敛阶，而在有的步长下不能。应该考虑时间步长和空间步长谁是主导因素，主导因素是否会发生转移。
一个高阶收敛的格式，刚开始的时候收敛阶可能表现得不大稳定，比如上面提到Lax-Wendroff
1 2 3 4 5 6 7 8
Lax P D . Weak solutions of nonlinear hyperbolic equations and their numerical computation[J]. Comm. Pure Appl. Math. 1954, 7(1):159-193. ↩︎ ↩︎
Shu C W. Numerical methods for hyperbolic conservation laws[J].↩︎
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. ↩︎
Anderson J D, Wendt J Computational fluid dynamics[M]. New York: McGraw-Hill, 1995. ↩︎
Landau L D, Lifshitz E M. Fluid Mechanics: Volume 6 (Pergamon International↩︎
庄礼贤,严协远,流体力学,中国科技大学出版社,1992 ↩︎
LeVeque R J. Finite difference methods for ordinary and partial differential equations:steady-state and time-dependent problems[M]. Siam, 2007. ↩︎
