I'm working on modeling a shock wave using the Euler equation with an advanced Equation of state and the fifth order WENO scheme. The equation are set up on the form: \begin{equation} \frac{\partial U}{\partial t} + \frac{\partial F(U)}{\partial x} = 0 \end{equation} where \begin{equation} U = \begin{bmatrix}\rho\\ \rho v\\E\end{bmatrix}, F(U) = \begin{bmatrix}\rho v\\ \rho v^2 + p\\v(E+p)\end{bmatrix} \end{equation} I've formulated the WENO scheme, but have yet to get it to work. What I find difficult to implement is the Lax-Friedrich Flux splitting. Specifically determining the $\alpha$ in the equation:

\begin{equation} f^{\pm} = f(U) \pm \alpha U \end{equation}

I do know that $\alpha = \max \left|\frac{\partial f(U)}{\partial U}\right|$, but for what index range of $U$ does the equation strech when i want to determine:

\begin{equation} \hat{f}_{i+1/2} = \hat{f}_{i+1/2}^+ + \hat{f}_{i+1/2}^- \end{equation}

  • $\begingroup$ The $\hat{f}$ denotes the numerical formulation. In regards to the Equation of state, it is a somewhat complicated EOS adapted to the LJS fluid. I utilize a function which computes the temperature and pressure based on the internal energy and density. This does makes it difficult to determine analytical eigenvalues, compared to the ideal case. $\endgroup$
    – Twm1995
    Feb 15, 2022 at 10:50
  • $\begingroup$ I solve the temporal discretization using 5th explicit Runge Kutta with the programmingfunction solve_ivp() in python. What I specifically want, is to compute the $f^{\pm}$ in order to utilize the WENO scheme, but descriptions in litterature of $\alpha$ doesn't indicate the range of $U$. The spatial formulation is correct. $\endgroup$
    – Twm1995
    Feb 15, 2022 at 12:28
  • $\begingroup$ Ah got it, so your semi-discrete formulation reads $$\frac{\mathrm d \boldsymbol{u}} {\mathrm d t } = \frac{1}{\Delta x} \Big( \hat{ \boldsymbol{f}}_{i-1/2} - \hat{ \boldsymbol{f}}_{i+1/2} \Big)$$? $\endgroup$
    – Dan Doe
    Feb 15, 2022 at 12:35
  • $\begingroup$ That is correct, where the $\hat{f}_{i\pm1/2}$ values are computed using the WENO scheme, with the splitted fluxes ($\hat{f}^{\pm}$). $\endgroup$
    – Twm1995
    Feb 15, 2022 at 12:41

1 Answer 1


I think you mean with Lax-Friedrichs the local Lax-Friedrichs or Rusanov scheme [1], [2] where the wave speeds (left & right-going) are given by the maximum eigenvalue $\lambda_\max := \max_i | \lambda_i(\boldsymbol{U}) | \in \sigma\Big( \boldsymbol{F}'(\boldsymbol{U}) \Big)$ where $\sigma$ denotes the spectrum of the Jacobian $\boldsymbol{F}'(\boldsymbol{U})$.

Note that the Rusanov / Local Lax Friedrich scheme is an approximate Riemann solver: For a conservation/balance law with Riemannian initial data, it provides an approximate solution. In particular, the scheme depends only on the trace (face / edge) values at the cell boundary / interface. In fact, it does not care how you got them - i.e., the scheme is agnostic whether you used simply the cell average value, a linear approximation with limiters or in your case the fifth order WENO scheme to reconstruct the face / trace values.

To finally answer your question, you need to take only the trace values $\boldsymbol{U}^-, \boldsymbol{U}^+$ into account at a certain face. Then,

$$\alpha := \max \bigg\{ \max_i \big \vert \lambda_i (\boldsymbol{U}^-) \big \vert \in \sigma\Big( \boldsymbol{F}'(\boldsymbol{U}) \Big) \Big \vert_{\boldsymbol{U}^-} , \max_i \big \vert \lambda_i (\boldsymbol{U}^+) \big \vert \in \sigma\Big( \boldsymbol{F}'(\boldsymbol{U}) \Big) \Big \vert_{\boldsymbol{U}^+} \bigg\}. $$ In other words at each face $i$ you jest need the trace values $\boldsymbol{U}_i^+, \boldsymbol{U}_i^-$. In that case, your numerical flux (I think you lack the factor if $1/2$ in your formulation and confused the signs) is given by \begin{align} \hat {\boldsymbol{f}}_{i+1/2}& = \frac{1}{2} \Big( \hat {\boldsymbol{f}}_{i+1/2}^+ + \hat {\boldsymbol{f}}_{i+1/2}^- \Big) =\frac{1}{2} \Big( F(\boldsymbol{U}_i^+) - \alpha \boldsymbol{U}_i^+ + F(\boldsymbol{U}_i^-) + \alpha \boldsymbol{U}_i^- \Big) \\ &= \frac{1}{2} \Big( F(\boldsymbol{U}_i^+) + F(\boldsymbol{U}_i^-) + \alpha (\boldsymbol{U}_i^- - \boldsymbol{U}_i^+ \Big) \end{align}

  • $\begingroup$ Thank you! A couple of follow-up questions: So by interface of an index $i$ I only take the maximum approximate eigenvalue between the $i \pm 1/2$ values? Do you have any advice for computing the jacobian i.e. approximate values for values connected to the EOS? $\endgroup$
    – Twm1995
    Feb 15, 2022 at 14:30
  • $\begingroup$ More precise, you take the maximum of the absolute values of the eigenvalues of the Jacobian $\frac{\mathrm d \boldsymbol{F}}{\mathrm d \boldsymbol{U}} = F'(\boldsymbol{U})$ which is evaluated at both sides of the interface $i+ 1/2$. Based on your reconstruction / interpolation method you have the value $\boldsymbol{U}_{i+1/2}^- $ on the left of the $i+1/2$ cell boundary and $\boldsymbol{U}_{i+1/2}^+ $ on the right. The Jacobian / eigenvalues (if analytically known) are evaluated at this states. $\endgroup$
    – Dan Doe
    Feb 16, 2022 at 7:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.