Accuracy of upwind method for discontinuous solutions

Consider the scalar advection equation

$$\frac{\partial u}{\partial t} + a\frac{\partial u}{\partial x} = 0, x \in {\mathbb{R}}, a > 0$$ (1)

with initial condition $u(x,0)=u_o(x)$. The upwind scheme for this pde is obtained by using backward difference formula for the spatial derivative and is given by

$$v^{n+1}_j = v^n_j - \frac{a\Delta t}{\Delta x}( v^n_j - v^n_{j-1})$$ (2)

We are interested in studying the accuracy of this method for discontinuous solutions by studying the solution of the modified partial differential equation (MPDE). Let us first consider the case of smooth solution and assume that there is a function $v(x,t)$ which agrees exactly with the numerical solution at the grid points,

$$v(x,t+\Delta t) = v(x,t) - \frac{a\Delta t}{\Delta x}[ v(x,t) - v(x-\Delta x, t) ]$$

We want to find the pde satisfied by $v(x,t)$. Expanding the terms in the above equation in Taylors series about $(x,t)$ and simplifying we get

$$v_t + a v_x = \frac{1}{2}(a\Delta xv_{xx} - \Delta tv_{tt}) - \frac{1}{6}( a \Delta x^2 v_{xxx} + \Delta t^2 v_{ttt}) + ....$$ (3)

This is the pde satisfied by the function $v$. If we take $\Delta t/\Delta x$ to be fixed (cfl condition), then the terms in the right hand side are $O(\Delta t), O(\Delta t^2)$, etc. so that for small $\Delta t$ we can truncate this series to get a pde that is quite well satisfied by the $v^n_j$.

If we drop all the terms on the right-hand side, we just recover the original advection equation. Since we have then dropped terms of $O(\Delta t)$, we expect the $v^n_j$ to be first order accurate.

Let us retain the $O(\Delta t)$ terms. Then we get

$$v_t + a v_x = \frac{1}{2}(a\Delta xv_{xx} - \Delta tv_{tt})$$ (4)

This involves second derivatives in both $x$ and $t$, but we can derive a slightly different modified equation with the same accuracy by differentiating equation (4) with respect to $t$ to obtain

$$v_{tt} = - av_{xt} + \frac{1}{2}(a\Delta xv_{xxt} - \Delta tv_{ttt})$$

and with respect to $x$ to obtain

$$v_{tx} = - av_{xx} + \frac{1}{2}(a\Delta xv_{xxx} - \Delta tv_{ttx})$$

Combining these gives

$$v_{tt} = a^2 v_{xx} + O(\Delta t)$$

Inserting this in (4) gives

$$v_t + a v_x = \frac{1}{2}(a\Delta xv_{xx} - a^2 \Delta tv_{xx}) + O(\Delta t^2)$$

Since we have already dropped terms of $O(\Delta t^2)$, we can drop these terms here also to obtain

$$v_t + av_x = \frac{1}{2}a\Delta x( 1 - \lambda)v_{xx}$$ (5)

where $\lambda = a\Delta t/\Delta x$ is the Courant number. This is the MPDE for the upwind scheme and is an advection-diffusion equation with diffusion coefficient or viscosity given by

$$\mu = \frac{1}{2}a\Delta x( 1 - \lambda)$$ (6)

• The grid values $v^n_j$ can be viewed as giving a second-order accurate approximation to the above pde, whereas they only give first-order accurate approximations to the advection equation (1).
• The fact that the modified equation for the upwind method is an advection-diffusion equation tells us a great deal about how the numerical solution behaves. Solutions to the advection-diffusion equation translate at the proper speed $a$ but also diffuse and are smeared out. This is also observed for the upwind scheme.
• The diffusion coefficient in equation (4) vanishes when $a\Delta t = \Delta x$. In this case we already know that the exact solution the advection equation is recovered by the upwind scheme.
• The diffusion coefficient is positive only if $0 < \lambda < 1$. This is precisely the stability limit of the upwind method. If it is violated, then the diffusion coefficient in the modified equation is negative, giving an ill-posed backward heat equation with exponentially growing solutions. Hence we see that some information about stability can also be obtained from MPDE.

Discontinuous solutions
The MPDE was derived on the assumption that the solution is smooth. When the exact solution is discontinuous, the numerical solution is still smooth because it is governed by the advection-diffusion equation.

Consider an initial condition with a jump discontinuity,

$$u_o(x) = 2 H(x)$$

where $H$ is the usual Heaviside function. The exact solution is

$$u(x,t) = 2 H(x - at)$$ (7)

that is the initial profile is advected with speed $a$ in the positive $x$ direction without change of shape. The exact solution of the advection-diffusion equation is given by

$$v(x,t) = \textrm{erfc}\left( \frac{ x - at }{\sqrt{4\mu t}} \right)$$ (8)

In figure the dashed line is the exact solution to the advection equation, the solid line is the solution (8) while the symbols are the upwind numerical solution. We see that the numerical solution gives a very good approximation to the MPDE. Hence the difference between equation (7) and (8) will give us an estimate of the error in the numerical solution. We will evaluate this error in the 1-norm,

$$\begin{eqnarray*} \vert\vert u - v \vert\vert _1 &=& \int^\infty_{-\infty} \vert... ...rt \textrm{d}x + \int_{at}^\infty \vert v(x,t)\vert \textrm{d} x \end{eqnarray*}$$

Now the solution $v$ is such that the part for $x > at$ can be obtained from the part $x < at$ by first reflecting it about $x=at$ and then about $v=1$. So that the two integrals in the previous equation are equal because they are just the two areas indicated in figure 2.

Figure 2

Hence the 1-norm error is

$$\begin{eqnarray*} \vert\vert v - u \vert\vert _1 &=& 2\int_{at}^\infty \vert v(x... ...t^\infty_0 \textrm{erfc}(z) \textrm{d}z \\ &=& C_1 \sqrt{\mu t} \end{eqnarray*}$$

for some constant $C_1$ independent of $\mu$, $\Delta x$ and $t$. Using the expression for $\mu$ we get the following estimate for the error

$$\vert\vert v - u \vert\vert _1 \approx C_2 \sqrt{t \Delta x}$$ (9)

as $\Delta x\to 0$ with $\Delta t/\Delta x$ fixed. This indicates that the 1-norm of the error decays only like $(\Delta x)^{1/2}$ even though the method is formally first order accurate based on the local truncation error, which is valid for smooth solutions only. We also see that the thickness of the transition zone (shock structure) increases with time like $t^{1/2}$.

Source
RJ LeVeque, Finite volume methods for hyperbolic problems, Cambridge Univ Press, 2002.