The 1D Burger's equation is a simple (if not the simplest) non-linear hyperbolic equation commonly used as a model equation to illustrate various numerical schemes for non-linear hyperbolic differential equations. It is normally presented as: $$ \begin{equation} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = 0 \tag{8.76} \end{equation} $$
To enable us to present schemes for a greater variety of hyperbolic differential equations and to better handle shocks (i.e discontinuities in the solution), we will present our model equation on conservative form: $$ \begin{equation} \frac{\partial u}{\partial t} + \frac{\partial }{\partial x} \left (\frac{u^2}{2} \right) = 0 \tag{8.77} \end{equation} $$ and by introducing a flux function $$ \begin{equation} F(u) = \frac{u^2}{2} \tag{8.78} \end{equation} $$
the conservative formulation of the Burger's equation may be represented by a generic transport equation: $$ \begin{equation} \frac{\partial u}{\partial t} + \frac{\partial F(u)}{\partial x} = 0 \tag{8.79} \end{equation} $$
First step: $$ \begin{align} & u_{j+1/2}^{n+1/2} = \frac{1}{2} \left (u_{j+1}^{n} + u_{j}^{n} \right ) - \frac{\Delta t}{2 \Delta x} \left (F(u_{j+1}^{n}) - F(u_{j}^{n}) \right ) \tag{8.82}\\ & u_{j-1/2}^{n+1/2} = \frac{1}{2} \left (u_{j}^{n} + u_{j-1}^{n} \right ) - \frac{\Delta t}{2 \Delta x} \left (F(u_{j}^{n}) - F(u_{j-1}^{n}) \right ) \tag{8.83} \end{align} $$ $$ \begin{equation} u_{j}^{n+1} = u_{j}^{n} - \frac{\Delta t}{\Delta x} \left ( F(u_{j+1/2}^{n+1/2}) - F(u_{j-1/2}^{n+1/2}) \right ) \tag{8.84} \end{equation} $$
In the previous section ((8.7.1 Lax-Wendroff for non-linear systems of hyperbolic PDEs)) we showed how the two step Lax-Wendroff method could be condensed to a one step method. The same procedure may be applied to a the general transport equation given by (8.79). However for the nonlinear case (8.43) no longer holds. This may be overcome be rewriting (8.43): $$ \begin{equation} \begin{array}{c} \dfrac{\partial u}{\partial t}\Bigg|_j^n = -\dfrac{\partial F}{\partial x}\Bigg|_j^n \qquad \text{ and} \qquad \dfrac{\partial^2u}{\partial t^2}\Bigg|_j^n = -\dfrac{\partial^2F(u)}{\partial t \partial x }\Bigg|_j^n = -\dfrac{\partial \left(\frac{\partial F(u)}{\partial t}\right)}{\partial x }\Bigg|_j^n \\ = -\dfrac{\partial \left(\frac{\partial F(u)}{\partial u} \frac{\partial u}{\partial t}\right)}{\partial x }\Bigg|_j^n = \dfrac{\partial \left(a(u) \frac{\partial F}{\partial x}\right)}{\partial x }\Bigg|_j^n \end{array} \tag{8.85} \end{equation} $$
Now inserting into the Taylor series we get $$ \begin{equation} u_{j}^{n+1} = u_{j}^{n} - \Delta t \frac{\partial F(u)}{\partial x} + \frac{\Delta t^2}{2} \dfrac{\partial \left(a(u) \frac{\partial F}{\partial x}\right)}{\partial x } \tag{8.86} \end{equation} $$
and further we may obtain the general Lax-Wendroff one-step method for a generic transport equation $$ \begin{equation} u_{j}^{n+1} = u_{j}^{n} - \frac{\Delta t}{2 \Delta x}\left(F_{j+1} - F_{j-1}\right) + \frac{\Delta t^2}{2 \Delta x^2} \Big[a_{j+1/2}\left(F_{j+1} - F_{j}\right) - a_{j-1/2}\left(F_{j} - F_{j-1}\right)\Big] \tag{8.87} \end{equation} $$
where \( a(u) \) is the wavespeed, or the Jacobian of F, \( F'(u) \), which is \( u \) for the burger equation. As indicated \( a(u) \) has to be approximated at the indices \( (j+1/2) \) and \( (j-1/2) \). This may simply be done by averaging the neighbouring values: $$ \begin{equation} a_{j+1/2} = \frac{1}{2}\left(u_{j}^{n} + u_{j+1}^{n}\right) \tag{8.88} \end{equation} $$ for the burger equation. Another method that assure conservation is to use the following approximation $$ \begin{equation} a_{j+1/2}= \left\{ \begin{array}{ll} \frac{F_{j+1}^n-F_{j}^n}{u_{j+1}^n-u_{j}^n} \quad if \quad u_{j+1}\neq u_{j} \\ u_{j} \quad \quad \quad otherwise \end{array} \right. \tag{8.89} \end{equation} $$
For the Advection equation we were able to verify our schemes by comparing with exact solutions, using MES. For the burger equation it is not easy to find an analytical solution, so in order to verify our schemes we use the MMS approach. However this requires that our schemes can handle source terms. The new equation to solve is thus the modified burgers equation $$ \begin{equation} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = Q \tag{8.93} \end{equation} $$
In this chapter we will consider source terms that are a function of \( x \) and \( t \) only. The basic approach to adding source terms to our schemes is to simply add a term \( Q_i^n \) to our discrete equations. The schemes mentioned above with possibility to handle source terms are summarized in Table (tab:Conservation_schemes).
| Name of Scheme | Scheme | order |
| Upwind | \( u_{j}^{n+1} = u_{j}^{n} - \frac{\Delta t}{\Delta x} + \Delta t Q_j^n \) | 1 |
| Lax-Friedrich | \( u_j^{n+1} = \frac{\Delta t}{2}(u^n_{j+1}+u^n_{j-1})-\frac{F^n_{j+1}-F^n_{j-1}}{2 \Delta x} + \Delta t Q_j^n \) | 1 |
| \( u_{j+1/2}^{n+1/2} = \frac{1}{2} \left (u_{j+1}^{n} + u_{j}^{n} \right )- \frac{\Delta t}{2 \Delta x} \left (F(u_{j+1}^{n}) - F(u_{j}^{n}) \right ) + \frac{\Delta t}{2} Q_{j+1/2}^n \) | ||
| Lax-Wendroff | \( u_{j-1/2}^{n+1/2} = \frac{1}{2} \left (u_{j}^{n} + u_{j-1}^{n} \right )- \frac{\Delta t}{2 \Delta x} \left (F(u_{j}^{n}) - F(u_{j-1}^{n}) \right )+ \frac{\Delta t}{2} Q_{j-1/2}^n \) | 2 |
| \( u_{j}^{n+1} = u_{j}^{n} - \frac{\Delta t}{\Delta x} \left ( F(u_{j+1/2}^{n+1/2}) - F(u_{j-1/2}^{n+1/2}) \right )+ \Delta t Q_{j}^n \) | ||
| macCormack | \( u_j^p =u_j^n - \frac{\Delta t}{\Delta x} \left (F_{j+1}^{n} -F_{j}^{n} \right ) + \Delta t Q_{j}^{n+1/2} \) | 2 |
| \( u_j^{n+1} =\frac{1}{2} \left (u_{j}^{n} + u_{j}^{p} \right ) - \frac{1}{2} \frac{\Delta t}{\Delta x} \left (F_{j}^{p} - F_{j-1}^{p} \right )+ \frac{\Delta t}{2} Q_{j}^{n+1/2} \) |
Examples on how to implement these schemes are given below, where we have used \( RHS(x,t) \) instead of \( Q \) for the source term:
ftbs or upwind:
def ftbs(u, t):
"""method that solves u(n+1), for the scalar conservation equation with source term:
du/dt + dF/dx = RHS,
where F = 0.5u^2 for the burger equation
with use of the forward in time backward in space (upwind) scheme
Args:
u(array): an array containg the previous solution of u, u(n). (RHS)
t(float): an array
Returns:
u[1:-1](array): the solution of the interior nodes for the next timestep, u(n+1).
"""
u[1:-1] = u[1:-1] - (dt/dx)*(F(u[1:-1])-F(u[:-2])) + dt*RHS(t-0.5*dt, x[1:-1])
return u[1:-1]
Lax-Friedrich's:
def lax_friedrich_Flux(u, t):
"""method that solves u(n+1), for the scalar conservation equation with source term:
du/dt + dF/dx = RHS,
where F = 0.5u^2 for the burger equation
with use of the lax-friedrich scheme
Args:
u(array): an array containg the previous solution of u, u(n). (RHS)
t(float): an array
Returns:
u[1:-1](array): the solution of the interior nodes for the next timestep, u(n+1).
"""
u[1:-1] = (u[:-2] +u[2:])/2.0 - dt*(F(u[2:])-F(u[:-2]))/(2.0*dx) + dt*(RHS(t, x[:-2]) + RHS(t, x[2:]))/2.0
return u[1:-1]
Lax-Wendroff-Two-step:
def Lax_W_Two_Step(u, t):
"""method that solves u(n+1), for the scalar conservation equation with source term:
du/dt + dF/dx = RHS,
where F = 0.5u^2 for the burger equation
with use of the Two-step Lax-Wendroff scheme
Args:
u(array): an array containg the previous solution of u, u(n).
t(float): time at t(n+1)
Returns:
u[1:-1](array): the solution of the interior nodes for the next timestep, u(n+1).
"""
ujm = u[:-2].copy() #u(j-1)
uj = u[1:-1].copy() #u(j)
ujp = u[2:].copy() #u(j+1)
up_m = 0.5*(ujm + uj) - 0.5*(dt/dx)*(F(uj)-F(ujm)) + 0.5*dt*RHS(t-0.5*dt, x[1:-1] - 0.5*dx) #u(n+0.5dt,j-0.5dx)
up_p = 0.5*(uj + ujp) - 0.5*(dt/dx)*(F(ujp)-F(uj)) + 0.5*dt*RHS(t-0.5*dt, x[1:-1] + 0.5*dx)#u(n+0.5dt,j+0.5dx)
u[1:-1] = uj -(dt/dx)*(F(up_p) - F(up_m)) + dt*RHS(t-0.5*dt, x[1:-1])
return u[1:-1]
macCormack:
def macCormack(u, t):
"""method that solves u(n+1), for the scalar conservation equation with source term:
du/dt + dF/dx = RHS,
where F = 0.5u^2 for the burger equation
with use of the MacCormack scheme
Args:
u(array): an array containg the previous solution of u, u(n). (RHS)
t(float): an array
Returns:
u[1:-1](array): the solution of the interior nodes for the next timestep, u(n+1).
"""
up = u.copy()
up[:-1] = u[:-1] - (dt/dx)*(F(u[1:]) - F(u[:-1])) + dt*RHS(t-0.5*dt, x[:-1])
u[1:] = .5*(u[1:] + up[1:] - (dt/dx)*(F(up[1:]) - F(up[:-1])) + dt*RHS(t-0.5*dt, x[1:]))
return u[1:-1]
In this exercise we want to assess the stability of the Lax-Friedrich scheme for the advection equation as formulated in (8.40).
a) Use the PC-criterion to find a sufficient condition for stability for (8.40).
b) Use the von Neumann method to find a sufficient and necessary condition for stability for (8.40) and compare with the PC-condition