A generalized numerical scheme for the linear hyperbolic equation $u_t + au_x = 0$ has the following form $$ \frac{u_j^{n+1} - u_j^n}{\Delta t} + a\frac{u_{j+1}^{n} - u_{j-1}^n}{2\Delta x} - \chi \frac{u_{j+1}^{n} - 2u_{j}^n + u_{j-1}^n}{\Delta x^2} = 0 . $$ a) Find the values of $\chi$ so that the above scheme is the central-explicit, forward, backward, Lax-Friedrich and Lax-Wendroff scheme.

b) What is the order of truncation for each scheme?

c) What is the stability condition for each of these schemes?

d) In what range of $\chi$ the scheme is unconditionally unstable?

e) Classify the schemes in (a) as (i) dissipative, and (ii) dispersive.


(a) Notice that if $\chi = 0$, then we obtain the central explicit scheme. If we write the numerical discretization as

$$u_j^{n+1} = u_j^n - \frac{ a \Delta t }{2 \Delta x} (u_{j+1}^n - u_{j-1}^n) + \frac{ \chi \Delta t }{\Delta x^2} (u_{j+1}^n - 2 u_j^n + u_{j-1}^n )$$

And one can see that if $\chi = a \Delta t/2$, then we have the LW scheme. If $\chi = a \Delta x/2$, then we have

$$ u_j^{n+1} = u_j^n + \frac{ a \Delta t }{\Delta x} (u_{j+1}^n - u_j^n) $$

which is forward explicit and if $\chi = - a \Delta x/2$, we obtain backward explicit. Now if $\chi = \frac{ \Delta x^2}{2 \Delta t }$, then we obtain Lax-Friedrichs.

We go to part (c) as (b) is trivial. Perhaps we can combine (c) and (d) in a single problem if we apply the discrete fourier transform:

$$ \hat{u}^{n+1} = \hat{u}^n - r( e^{ij \xi} -e^{-ij \xi } ) \hat{u}^n + \frac{ 2 \chi r }{a \Delta x} (e^{ij \xi} -e^{-ij \xi } -2 ) \hat{u}^n$$

Now, am I on the right track here? Is this the correct approach? As for e), I need some suggestions, how would we approach it?


The proposed values of $\chi$ in a) are correct. One notes that the method is a straightforward finite-difference approximation of the advection-diffusion method $u_t + au_x = \chi u_{xx}$. Introducing the Courant number $r = |a| {\Delta t}/{\Delta x}$, the diffusion constants from a) rewrite as $$ \chi \in |a|\frac{\Delta x}{2} \times \left\lbrace \frac{1}{r},\, 1,\, r,\, 0 \right\rbrace , $$ which correspond to the numerical viscosities of the Lax-Friedrichs, upwind, Lax-Wendroff, and central explicit method, respectively. One notes that if $0<r<1$, then those values of $\chi$ are sorted in decreasing order.

As done in this post, let us assume a perturbation of the form $u_j^n = \xi^n e^{-\text i k j \Delta x}$, where $\xi = \text e^{\text i \omega \Delta t}$. Note that $\omega$ is a complex number, with real part $\omega_R$ and imaginary part $\omega_I$. Injecting this Ansatz in the time-stepping formula, Euler's formulas lead to $$ \xi = 1 + \text i r \sin( k \Delta x) + 2\lambda \left( \cos(k \Delta x) -1 \right) , $$ where $r$ is the Courant number and $$ \lambda = {\chi}\frac{\Delta t}{\Delta x^2} \in \left\lbrace \frac{1}{2},\, \frac{r}{2},\, \frac{r^2}{2},\, 0 \right\rbrace $$ denotes the Fourier number. From this equation, one equates the real and imaginary parts: \begin{aligned} \text{Re}\, \xi &= e^{-\omega_I \Delta t}\cos(\omega_R \Delta t) = 1 + 2\lambda \left( \cos(k \Delta x) -1 \right) ,\\ \text{Im}\, \xi &= e^{-\omega_I \Delta t}\sin(\omega_R \Delta t) = r \sin( k \Delta x) \, . \end{aligned} The stability can be analyzed in terms of the squared amplification factor $|\xi|^2$ which has to be smaller than one. Hence, we must have $$ |\xi|^2 - 1 = 4\lambda^2 - 4 \lambda + r^2 + 4 \lambda (1-2\lambda) Z + (4\lambda^2 - r^2) Z^2 \leq 0 $$ for all $Z = \cos( k \Delta x)$ in $[-1, 1]$. The dispersion is analyzed in terms of the phase velocity ${\omega_R}/{k}$ with ${\omega_R}\Delta t = \arctan\left[{\text{Im}\, \xi}/{\text{Re}\, \xi}\right]$. The attenuation is analyzed in terms of ${\omega_I}\Delta t = -\ln |\xi|$ which is non-negative if $|\xi| \leq 1$ (cf. e.g. (1) p. 182).

(1) E. Godlewski, P.-A. Raviart, Numerical Approximation of Hyperbolic Systems of Conservation Laws, Springer, 1996. doi:10.1007/978-1-4612-0713-9

  • $\begingroup$ You are using the fourier modes to analyze stability, correct? Was there any particular reason to use this? We could use the discrete fourier transform as well. $\endgroup$ – James May 17 at 17:36
  • $\begingroup$ @James I guess that DFT would work as well. $\endgroup$ – Harry49 May 17 at 17:39
  • $\begingroup$ If you dont mind asking, can you explain a little bit more how to analyze the dispersion/dissipation? $\endgroup$ – James May 17 at 17:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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