von Neumann analysis

Here after, flux of a Finite Volume method are defined for a convection equation from left to right.

• upwind1 is a 1st order extrapolation of left cell point

• center2 is a 2-points centered interpolation at the face

• fromm and quick are 2 2nd order 3-points interpolations with different weights

• upwind3 is the only 2rd order interpolation with 3 points (left shifted parabolic interpolation)

• center4 is a 4th order 4-points center interpolation

import numpy as np

def flux(q, method):
    def fluxk(q, k):
        return q + .25*(1+k)*(np.roll(q,-1)-q) + .25*(1-k)*(q-np.roll(q,1))
    fluxmeth = {
        'upwind1' : lambda q: q,
        'center2' : lambda q: .5*(q+np.roll(q,-1)),
        'fromm' : lambda q: fluxk(q, 0.),
        'quick' : lambda q: fluxk(q, .5),
        'upwind3' : lambda q: fluxk(q, 1./3.),
        'center4' : lambda q: .5*(q+np.roll(q,-1)) - .125*(2./3.)*(np.roll(q,-2)-np.roll(q,-1)-q+np.roll(q,1))
    }
    if method not in fluxmeth.keys():
        KeyError(f"'{method}' unknown")
    return fluxmeth[method](q)

Thought this flux can be used in an actual integration, they are here used to evaluate the complex discretized propagator \(P(kh)\) which is compared to the theoretical one \(e^{-Ikh}\).

def error(kh, method):
    imax = 2 # stencil of 5 points
    cshift = np.exp(1j*kh) # complex shift to get +1 index
    q = np.array([cshift**(i-imax) for i in range(2*imax+1)]) # q[imax] = 1
    rhs = -(flux(q, method)-flux(q/cshift, method)) # rhs on imax index
    err = rhs[imax] - (-1j*kh)
    return abs(err.real)+1j*abs(err.imag)

One can then plot the real and imaginary parts of this error that will be named dissipation and dispersion errors. The errors are plotted against \(\lambda/h = 2\pi k h\) which is the number of points per wavelength.

import matplotlib.pyplot as plt
fig, (axr, axi) = plt.subplots(1, 2, figsize=(10, 4))
kh = np.geomspace(5e-2, np.pi, 100)
for ax, ylab in zip([axr, axi], ['dissipation error', 'dispersion error']):
    ax.set_xlabel('$\lambda/h$') ; ax.set_ylabel(ylab)
    ax.grid()
    ax.set_ylim(1e-4, 1e1)
for method, sty in zip(['upwind1', 'center2', 'fromm', 'quick', 'upwind3', 'center4'],
                       ['-g', '--b', '-r', '-c', '-m', '--y']):
    err = np.array([error(ikh, method) for ikh in kh])
    axr.loglog(2*np.pi/kh, err.real, sty, label=method)
    axi.loglog(2*np.pi/kh, err.imag, sty, label=method)
axr.legend() ; axi.legend();

Caution

centered schemes have no dissipation error Due tu the symmetry property, the real part of the error is zero. These schemes provide only dispersion errors. It can achieve a higher accuracy but this may be not suitable for steep signals.

parity of dissipation and dispersion errors

The dissipation error is even and the dispersion is odd. That is the reason why

• fromm and quick schemes have the dissipation error trend as 3rd order upwind3

• upwind1 and center2 have the same 2nd order dispersion error

• upwind3 and center4 have the same 4th order dispersion error