A conservative scheme gives correct shock location

I learnt about this proof from Prof. Gerald Warnecke when he was visiting IISc, Bangalore, sometime in 2002.

Consider a single conservation law of the form

$$\frac{\partial u}{\partial t} + \frac{\partial}{\partial x} f(u) = 0$$ (1)

Let

$$f_{j+1/2} = f(u_{j-l},....,u_j, u_{j+1}, ...., u_{j+k})$$ (2)

be a numerical flux function which is consistent in the sense that

$$f(u_{j-l},....,c, c, ...., u_{j+k}) = f(c), \forall c \in \Re$$ (3)

Then the first order in time update can be written as

$$u^{n+1}_j = u^n_j - \lambda( f^n_{j+1/2} - f^n_{j-1/2}), \lambda = \Delta t/\Delta x$$ (4)

with the quantity $u_j$ being the cell average value. Let us consider a step function as an initial condition

$$u_o(x) = \left\{ \begin{array}{cc} u_l, & x \le 0 \\ u_r, & x > 0 \end{array} \right.$$ (5)

Figure 1:

After time $t$ the profile is shown in figure (1) and the exact shock location is $x_e(t)$, which is given by

$$x_e(t) = \frac{ \textrm{shaded area} }{ \textrm{jump in } u ... ...frac{ \int_{-\infty}^\infty [ u(x,t) - u_o(x) ] }{ u_l - u_r }$$ (6)

Note that for finite time $t$, $u(\cdot, t) - u_o$ has compact support so that the integral makes sense. We now apply the same definition to the numerical solution, except that now the integral is replaced by a summation

$$x_{num}(n+1) - x_{num}(n)= \frac{ \Delta x \sum_{j=-\infty}^\infty (u^{n+1}_j - u^n_j ) }{u_l - u_r}$$ (7)

Using (4) we get

$$x_{num}(n+1) - x_{num}(n)= - \frac{ \Delta t \sum_{j=-\infty}^\infty ( f^n_{j+1/2} - f^n_{j-1/2}) }{ u_l - u_r }$$ (8)

We next assume that the scheme is TVD (this is not really necessary). Then for $j$ very large or very small, we have by flux consistency

$$f_{j+1/2} = f_{j-1/2}$$

Then (8) becomes (see figure (2))

$$x_{num}(n+1) - x_{num}(n)= - \frac{ \Delta t \sum_{j=j_1}^{j_2} ( f^n_{j+1/2} - f^n_{j-1/2}) }{ u_l - u_r }$$ (9)

Figure 2:

From figure (2) we see that there is a cancellation of the fluxes and the previous equation reduces to

$$x_{num}(n+1) - x_{num}(n)= - \frac{ \Delta t (-f^n_{j_1 - 1/2} + f^n_{j_2 + 1/2}) }{ u_l - u_r }$$ (10)

But by consistency of the flux function

$$f^n_{j_1 - 1/2} = f(u_l), f^n_{j_2 + 1/2} = f(u_r)$$

Hence we get

$$\begin{eqnarray*} x_{num}(n+1) - x_{num}(n) &=& - \frac{ \Delta t (-f_l + f_r) }... ...( \frac{ f_l - f_r }{u_l - u_r } \right)\\ &=& \Delta t \cdot S \end{eqnarray*}$$

where $S$ is the exact shock speed. We thus see that the numerical shock moves with the correct speed and hence has the correct shock location. The constraint of TVD can be easily removed. Even if there are wiggles in the solution we can still find a $j_1, j_2$ such that solution is constant for large $\vert j\vert$ and the proof will still hold. Also the solution need not be of a step function and can have any shape except that there should be only one discontinuity. The remarkable fact about this result is that ALL conservative schemes will give correct shock location; CONSERVATION is the only condition required.