How can I derive a bound on the spurious oscillations in the numerical solution of the 1D advection equation?

Question

Suppose I had the following periodic 1D advection problem:

$\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0$ in $\Omega=[0,1]$
$u(0,t)=u(1,t)$
$u(x,0)=g(x)$
where $g(x)$ has a jump discontinuity at $x^*\in (0,1)$.

It is my understanding that for linear finite difference schemes of higher than first order, spurious oscillations occur near the discontinuity as it is advected over time, resulting in a distortion of the solution from its expected wave shape. According to wikipedia explanation, it seems that these oscillations typically occur when a discontinuous function is approximated with a finite fourier series.

For some reason, I can't seem to grasp how a finite fourier series can be observed in the solution of this PDE. In particular, how can I estimate a bound on the "over-shoot" analytically?

Answer 1

The first order upwind method is monotone; it does not introduce spurious oscillations. But it is only first order accurate, resulting in so much numerical diffusion as to be unusable for many purposes. Godunov's Theorem states that linear spatial discretizations of higher than first order cannot be monotone. To rigorously control oscillations, we use Total Variation Diminishing (TVD) schemes. TVD methods are typically limited to second order accuracy. For higher order, we must either relax our request, leading to Total Variation Bounded (TVB) methods like (Weighted) Essential Non-Oscillatory ((W)ENO), or we must relax the definition of TVD to "maximum-principle preserving" or similar, where the initial extrema are in terms of an initial reconstructed solution, resulting in special limiting schemes.

Answer 2

Linear finite difference discretization of a 1D problem with periodic boundaries leads to a discretization of the form $$U^{n+1} = LU^n$$

where $L$ is a circulant matrix. The eigenvectors of any circulant matrix are discrete Fourier modes $$v_j = \exp(ijh\xi)$$ (here $h$ is the grid spacing and $\xi$ is the wavenumber, which ranges from zero up to the highest wavenumber representable on the grid). These eigenvectors form a basis for all functions that can be represented on the grid. If you express the solution in terms of these discrete Fourier modes, then the numerical method is diagonalized, i.e. each Fourier component is multiplied by a (generally complex) scalar factor at each step. The scalar factor is often referred to as the amplification factor, and what I have just described is known as von Neumann analysis. It is analogous to Fourier analysis of linear PDEs, in which one uses a Fourier basis, to "diagonalize" the linear differential operators.

You can find nice explanations, for instance, in the text of Strikwerda or LeVeque.

Answer 3

Not all spurious oscillations are Gibbs phenomena. They look similar, but there are Gibbs oscillations for all finite Fourier approximations of discontinuous functions (they just get smaller as you add more terms). Whereas, there are non-oscillatory representations of discontinuous functions resulting from the solution of finite difference approximations to PDEs that don't require infinite series.

Bathe (Inf–sup testing of upwind methods, PDF) has a paper on this for finite element methods (convection-diffusion, IIRC) in 1-D that involves computing the constant for the $\inf$-$\sup$ condition and relating that to oscillations. You might gain some insight from that.

Answer 4

As for your last question about the connection between finite Fourier series and finite element approximation: In general, if you try to project a function with a jump onto a finite dimensional space whose basis functions are continuous, you get Gibbs phenomenon. This is true if the basis is a finite Fourier series (where the basis functions are the sines and cosines) or if the basis are the usual finite element hat functions -- it's a property of the projection plus the unsuitability of the basis functions.

Answer 5

One approach is via the equivalent equation, that is, the differential equation to which your discrete method gives the closest aproximation. This is never the differential equation that you intended to solve. Then you look at the asymptotic solution of the equivalent equation, for a step function as initial data. Look at Bouche, D., Bonnaud, G. and Ramos, D., 2003. Comparison of numerical schemes for solving the advection equation. Applied mathematics letters, 16(2), pp.147-154.