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

$\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?

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.

• My apologies... for some reason, I got the impression that this was true for the first order scheme as well. I edited the question to reflect this comment.
– Paul
May 28, 2012 at 6:37

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.

• I'm familiar with von neumann analysis. But can I really use this analysis to derive a bound on the spurious oscillations?
– Paul
May 28, 2012 at 6:22
• I was responding mainly to your statement I can't seem to grasp how a finite fourier series can be observed in the solution of this PDE. But yes, you can get such bounds from this analysis. For instance, you could look at the worst case secenario in which all modes constructively interfere. However, this is likely to be a very pessimistic bound. In practice, I haven't seen anyone derive bounds other than TVD or TVB (which are quite strong and don't hold for linear schemes). May 28, 2012 at 6:25
• You could probably get a more interesting bound by looking at the dispersion relation for the highest wavenumber modes. But I've never seen it done. May 28, 2012 at 6:28

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.

• This is a useful paper, but note that inf-sup stability does not provide strong control of oscillations. No amount of inf-sup stability can provide a TVD method, for example. And in light of Godunov's Theorem, it makes no sense to go looking for linear spatial discretizations if we intend to have non-oscillatory solutions of greater than first order. Note that the Peclet number appears in all the methods in this paper, and the methods degrade to first order accuracy as $\mathrm{Pe}\to \infty$, while also not being TVD. May 27, 2012 at 19:00
• These are all true statements. It only really applies to convection-diffusion problems. May 27, 2012 at 21:59

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.

• I'm happy to be proved wrong, since I'm clearly out of practice, but I'm not buying your comment about projections onto hat functions without further qualification. My quick computation using my old 1-D MATLAB code from my first year FEM class shows that the projection of the step function onto $H_0^1$ using hat functions is non-oscillatory. Do you have an example that can show what I'm missing? May 27, 2012 at 22:27
• Never mind. Old code is old. I can reproduce oscillations. Previous comment retracted. May 27, 2012 at 22:45
• I'm glad I could help :-) May 29, 2012 at 17:33

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.