Van Leer’s analysis of these methods is a tour-de-force especially considering the sort of resources available to him at that time.  He did not have use of tools such as Mathematica upon which I rely.  Nonetheless, the modern tools do allow more flexibility in analysis, and exploration of extensions of the methods without undue effort.  This is not to say that the analysis is particularly easy as Mathematica analysis goes, it does require considerable skill, experience and attention.  My hope for this rearticulation of the results is for clarity with a substantial appreciation for the skill expressed in the original exposition.

Here we briefly review the six schemes presented by Van Leer in that 1977 paper in terms of their basic structure, and numerical properties.  We will utilize Von Neumann analysis for both the fully discrete and semi-discrete schemes.  In Von Neumann analysis the variable on the mesh is replaced by the Fourier transform,

$u_j = \exp(\imath j \theta)$ or $u_j^n = A^n\exp(\imath j \theta)$

with $\theta$ being the angle and $A$ being the amplification factor.

The fully discrete analysis was the vehicle of presention for Van Leer, and the semi-discrete analysis will be additional detail.   The fully discrete schemes all fufill the characteristic condition, they recover the exact solution as the CFL number goes to one, which also demarks its’ stability limit.  The semi-discrete analysis only provides the character of the spatial scheme.

Scheme I. The Piecewise Linear Method.

This method uses a piecewise linear reconstruction

$P(\xi) = P_0 + P_1 \xi$ with $\xi = (x-x_j)/\Delta x$, $P_0 = u_j^n$ with slopes computed from local data $P_1 = \frac{1}{2} (u_{j+1} - u_{j-1})$.

The full update is

$u^{n+1}_{j} = u^{n}_{j} - \nu (u^{n+1/2}_{j+1/2} - u^{n+1/2}_{j-1/2})$

with $u^{n+1/2}_{j+1/2} = P(1/2-1/2\nu)$ and $\nu=1 \Delta t/\Delta x$ is the CFL number.

This method is second-order accurate as verified by expanding the method function in a Taylor series.  This can also be verified by expanding the Fourier transformed mesh function in a Taylor series expansion around $\theta = 0$.

Scheme I amplification surface plot, shows large error at the mesh scale

The amplification factor is $A =1-\frac{1}{8} \left(\nu-2 \nu^2+2\nu^3-\nu^4 \right)\theta^4+O\left(\theta^6\right)$

The phase error is $\phi = 1+\frac{1}{12}(1-2\nu+3\nu^2)\theta^2 + O(\theta^4)$ this is the leading order error (being one order higher due to how the phase error is written).  A key aspect of this scheme has zero phase error at CFL numbers of one-half and one.

As with the amplification, scheme I has large errors at the mesh scale.

The semi-discrete form is $u_{j+1/2} - u_{j-1/2} = \imath \theta + \frac{1}{12}\imath \theta^3 + \frac{1}{8}\theta^4+ O(\theta^5)$, which is useful if the scheme is used with a method-of-lines time integration.

Dispersion error for Scheme I is large past $pi/2$

Scheme II. The Piecewise Linear Method with Evolving Slopes.

This method evolves both the variable and its first derivative and requires the analysis of a system of variables.  The scalar analysis of stability is replaced by the analysis of a matrix determined by the behavior of its eigenvalues.  One eigenvalue will define the behavior of the scheme, while the second eigenvalue has be referred to as spurious.  More recently it has been observed that the “spurious’’ eigenvalue describes the subcell behavior at the mesh scale.

The method uses the same reconstruction as Scheme I, but evolves both $P_0 = u_j^n$ and $P_1 = s_j^n$  in time.  The update is then two equations, one for the cell-centered data,

$u^{n+1}_{j} = u^{n}_{j} - \nu (u^{n+1/2}_{j+1/2} - u^{n+1/2}_{j-1/2})$,

which is identical to Scheme I and a second for the slopes

$s^{n+1}_j = u_j^n+\frac{1}{2} s_j^n-u_{j-1}^n - \frac{1}{2} s_{j-1}^n - \nu(s_j^n - s_{j-1}^n)$,

where the usual $s_j^n$ as the initial slope is replaced by the first expression on the right hand side of the equation driving a coupling of the variable and the slope in the update.

The slope update is unusual and causes significant error as the CFL number goes to zero.  This method does not have an obvious method of lines form due to the nature of this update.

The amplification factor is $A =1-\frac{1}{8} \left(\nu^4-2 \nu^3+2 \nu^2\right)\theta^4+O\left(\theta^6\right)$.

The phase error is $\phi = 1+\frac{1}{12}(1-3\nu+2\nu^2)\theta^2 + O(\theta^4)$ and is the leading order error as with all of the first three methods. Again, this method does not lend itself to a well-defined semi-discrete method.

Scheme III. The Piecewise Linear Method with Evolving First Moment.

This method uses a reconstruction that is similar to piecewise linear, but instead of slope (first derivative), the linear term is proportional to the moment of the solution,

$P(\xi) = P_0 + P_1 \xi$

again $P_0=u_j^n$, but $P_1 = m_1 = \int_{-1/2}^{1/2} P(\xi)\xi d\xi$.  Here, as with Scheme II, we evolve both $P_0$, and $P_1$, but the evolution equation for $P_1$ is much different, and naturally yields a coupled system.  This also provides a direct path to writing the scheme in a semi-discrete form.

The update for $P_0$ is the same as the previous two schemes,

$u^{n+1}_{j} = u^{n}_{j} - \nu (u^{n+1/2}_{j+1/2} - u^{n+1/2}_{j-1/2})$.

The update for $P_1$ is based on the weak form of the PDE and uses the definition of the first moment,

$\int_{-1/2}^{1/2} \xi(\partial_t u + \partial_x u)d\xi =0$, $\partial_t m - \int_{-1/2}^{1/2} \partial_x \xi u d\xi + [u(1/2)+u(-1/2)] = 0$.

Discretely this equation can be evaluated for the scalar wave equation to

$m_j^{n+1} = m_j^n - \frac{1}{2}\nu (u^{n+1/2}_{j-1/2}+u^{n+1/2}_{j+1/2}) + \frac{1}{6}(u_j^n + 4 u_j^{n+1/2} + u_j^{n+1})$, where $u^{n+1/2}_{j} = u^{n}_{j} - \frac{1}{2}\nu (u^{n+1/2}_{j+1/2} - u^{n+1/2}_{j-1/2})$.

The amplification factor is $A =1-\frac{1}{72} \left(\nu-2 \nu^2+2 \nu^3-\nu^4\right)\theta^4+O\left(\theta^6\right)$.

The phase error is $\phi = 1+\frac{1}{540}(2-5\nu+5\nu^3-2\nu^4)\theta^4 + O(\theta^4)$ now the amplification error is leading and the method is confirmed as being third-order. This scheme is obviously more accurate than the previous two methods.

The semi-discrete form follows from the standard updates for both the variable and its moment without any problem.  The update has an accuracy of second order with small coefficients, $u_{j+1/2} - u_{j-1/2} = \imath \theta + \frac{1}{72} \theta^4 + \frac{1}{540}\imath \theta^5+ O(\theta^6)$, which the keen observer will notice being third-order accurate instead of the expected second order.  Van Leer made a significant note of this.

Scheme IV. The Piecewise Parabolic Method.

This method uses a piecewise parabolic reconstruction

$P(\xi) = P_0 + P_1 \xi + P_2(\xi^2 - \frac{1}{12})$,

and uses the variables and its two nearest neighbors to determine the fit.  The fit is defined as recovering the integral average in all three cells.  The coefficients are then simply defined by $P_0=u_j^n$$P_1=\frac{1}{2}(u_{j+1}^n - u_{j-1}^n)$, and $P_2=\frac{1}{2}(u_{j+1}^n - 2u_j^n +u_{j-1}^n)$ . This is a third-order method.

The update is simple based on the evaluation of edge and time-centered values that need to be considered slightly differently for third-order methods,

$\bar{u}_{j+1/2} = \int_{j+1/2}^{j+1/2-\nu} P(\xi) d\xi = P_0 + P_1(\frac{1}{2} - \frac{\nu}{2}) + P_2 (\frac{1}{2} - \frac{\nu}{2} + \frac{\nu^2}{3})$.  The update then is simply the standard conservative variable form,

$u_j^{n+1} = u_j^n -\nu(\bar{u}_{j+1/2} - \bar{u}_{j-1/2})$ .

The amplification factor is $A =1-\frac{1}{24} \left(\nu^2-2 \nu^3+\nu^4 \right)\theta^4+O\left(\theta^6\right)$

The phase error is $\phi = 1-\frac{1}{60}(2-5\nu+5\nu^3-2\nu^4)\theta^4 + O(\theta^6)$ this is the leading order error (being one order higher due to how the phase error is written).  A key aspect of this scheme has zero phase error at CFL numbers of one-half and one.

The semi-discrete form is $u_{j+1/2} - u_{j-1/2} = \imath \theta + \frac{1}{12}\theta^4 - \frac{1}{30}\imath\theta^5 + O(\theta^6)$, which is useful if the scheme is used with a method-of-lines time integration.

Scheme V. The Piecewise Parabolic Edge-Based Method.

This method uses a piecewise parabolic reconstruction

$P(\xi) = P_0 + P_1 \xi + P_2(\xi^2 - \frac{1}{12})$,

and uses the variables and its two edge values determining the fit.  The fit recovers the integral of the cell $j$ and the value at the edges $j\pm1/2$.  The coefficients are then simply defined by $P_0=u_j^n$, $P_1 =u_{j+1/2}^n - u_{j-1/2}^n$, and $P_2 = 3(u_{j+1/2}^n-2u_j^n+u_{j-1/2}^n)$ . This is a third-order method.

This method is distinctly different than Scheme IV in that both the cell-centered, conserved quantities, and the edge-centered (non-conserved) quantities.  The form of the updates are distinctly different and derived accordingly.  The time-integrated edge value follows the polynomial form from Scheme IV,

$\bar{u}_{j+1/2} = \int_{j+1/2}^{j+1/2-\nu} P(\xi) d\xi = P_0 + P_1(\frac{1}{2} - \frac{\nu}{2}) + P_2 (\frac{1}{2} - \frac{\nu}{2} + \frac{\nu^2}{3}) = u_{j+1/2}^n -\frac{\nu}{2} P_2 -(\frac{\nu}{2} - \frac{\nu^2}{3})P_3$,

giving a cell centered update, $u_j^{n+1} = u_j^n -\nu(\bar{u}_{j+1/2} - \bar{u}_{j-1/2})$ .  The edge update is defined by a time-integrated, upstream-biased first derivative evaluated at the cell-edge position,

$u_{j+1/2}^{n+1}= u_{j+1/2}^n - \nu P_2-\nu P_3+\nu^2 P_3$ .

The amplification factor is $A =1-\frac{1}{72} \left(\nu-2 \nu^2+2 \nu^3-\nu^4\right)\theta^4+O\left(\theta^6\right)$.

The phase error is $\phi = 1+\frac{1}{540}(2-5\nu+5\nu^3-2\nu^4)\theta^4 + O(\theta^4)$ now the amplification error is leading and the method is confirmed as being third-order. The scheme is very close in performance and formal error form to Scheme III, but is arguably simpler.

The semi-discrete form follows from the standard updates for both the variable and its moment without any problem.  The update has an accuracy of third order with small coefficients, $u_{j+1/2} - u_{j-1/2} = \imath \theta + \frac{1}{72} \theta^4 + \frac{1}{270}\imath \theta^5+ O(\theta^6)$ .

Looking at the recent PPML (piecewise parabolic method – local) scheme of Popov and Norman produces an extremely interesting result.  The interpolation used by PPML is identical to Scheme V, but there is a key difference in the update of the edge values.  The differential form of evolution is not used and instead a semi-Lagrangian type of backwards characteristic interpolation is utilized,

$u^{n+1}_{j+1/2} = P_2(1/2 -\nu)$ .

This produces exactly the same scheme as the original Scheme V for the scalar equation.  In many respects the semi-Lagrangian update is simpler and potentially more extensible to complex systems. Indeed this may be actualized by the PPML authors who have applied the method to integrating the equations of MHD.

Scheme VI. The Piecewise Parabolic Method with Evolving First and Second Moments.

This method uses a piecewise parabolic reconstruction

$P(\xi) = P_0 + P_1 \xi + P_2(\xi^2 - \frac{1}{12})$,

and uses the variables and its two edge values determining the fit.  The fit recovers the integral of the cell $j$ and the value at the edges $j\pm1/2$.  The coefficients are then simply defined by $P_0=u_j^n$, $P_1 =12 m_{1,j}^n$, and $P_2= 180 m_{2,j}^n$ $- 15 u_j^n$ . $m_1$ is the first moment of the solution, and $m_2$ is the second moment.  The moments are computed by integrating the solution multiplied by the coordinate over the mesh cell.  This is a fifth-order method.

For all the power Mathematica provides, this scheme’s detailed analysis has remained elusive, so the results are only presented for the semi-discrete case where Mathematica remained up to the task.  Indeed I did get results for the fully discrete scheme, but the plot of the amplification error was suspicious and for that reason will not be shown, nor will the detailed truncation error associated with the Taylor series.

The semi-discrete form follows from the standard updates for both the variable and its moment without any problem.  The update has an accuracy of second order with small coefficients, $u_{j+1/2} - u_{j-1/2} = \imath \theta + \frac{1}{720} \theta^6 + \frac{1}{4200}\imath \theta^7+ O(\theta^6)$, which the keen observer will notice being fifth-order accurate instead of the expected third order.  The semi-discrete update for the first moment is $\frac{1}{2}(u_{j+1/2} + u_{j-1/2}) - u_j$, and the second moment is

$1/4(u_{j+1/2} - u_{j-1/2}) - 2 m_{1,j}$.