This part analyses the accuracy and the stability of presenting scheme in solving the linear advection equation. For concreteness, we only give a detailed analysis for RDO-4 since the generalized RDO-M can be analyzed following the same routine.
4.1. Formulation of the Amplification Matrix
Consider the linear scalar advection equation
Suppose the spatial domain is
and is discretized by a uniform mesh spacing
. The initial condition is periodic,
where
denotes
and the scaled wavenumber is
. Then the initial model variables for the upwind RDO-4 on
are
The core in Fourier analysis is reformulating the algorithm into the matrix form
using the periodicity of
and
denotes the amplification matrix. We do this step by step as follows.
Stage 1 described in
Section 2 reconstructs the Hermite polynomial coefficients
from the rescaled model variables
in the cell
,
where
is defined by (
5) and
is the composition of the model variables after the coordinate transform,
In Stage 2, the solution values and derivatives on solution points are, respectively,
and
where
and
can be found in (
13) and (
10). According to (
31), the temporal derivatives on the solution points are
Then the temporal derivatives on solution points are transformed back to the temporal derivatives of model variables, which is the inverse problem of (
37):
where the subscript
j in the left-hand side indicates that the temporal derivatives are computed locally form
j-th cell, and
.
In the linear case,
is an invariant about the choice of solution points. Furthermore, the matrix
has a simple structure
Therefore, it can be easily concluded that RDO schemes using different sets of solution points are identical in solving the linear advection equation as the following theorem tells.
Theorem 1. For RDO-M schemes (), the numerical result is independent of the solution points chosen in solving the linear advection equation on uniform meshes.
Proof. This can be proved by showing that
is invariant about the choice of solution points
. In fact, one can easily check that
Moreover, it is noticed that the in-cell temporal derivatives (
39) only relies on
,
and
.
and
are independent with the choice of solution points. Therefore, we can conclude that the semi-discrete scheme is independent of solution points chosen. □
Equation (
39) gives the in-cell numerical approximation of derivatives of model derivatives. Denote the temporal derivatives
. From the periodicity of the initial profile, we know that
Then for the upwind discretization, the periodicity of the initial profile (
42) tells that
Hence, we can obtain the amplification matrix
for the upwind scheme as follows,
with
Similarly, the temporal derivatives of
for the central RDO-4 are
Combining (
39) and (
47) we have
where
4.2. Accuracy Analysis
We can see that both the central and the upwind RDO-4 have fourth spatial local discretization errors. However, this does not suffice to ensure the third order accuracy in numerical experiments. In fact, the order accuracy is determined by the principal eigenvalues
according to [
18], with
being the principle eigenvalue of the amplification matrix
. Other eigenvalues of
are spurious ones. Specifically, if
can be expanded as a Taylor series in
w,
then the scheme is
m-th order accurate. Since it is difficult to calculate an expression manually in practice, we use the symbolic computation tools in Matlab
.
For the upwind RDO-4, we can obtain the principal eigenvalue of the amplification matrix
is
This indicates that the upwind RDO-4 is third order accurate in space. Besides, the term introduces a small amount of dissipation which tends to smear the solution, which is certificated in the numerical experiments.
The principal eigenvalue of the amplification matrix
for central RDO-4 is
and we can see that only second order spatial accuracy is maintained for the central scheme. It is also worthy noting that
introduces slight dispersion, which means that phase speed varies as
w changes. Therefore, during the long-term simulation of a complicated initial profile that contains waves of different frequency, the profile will be distorted finally. This is also verified in numerical experiments in
Section 5.
The accuracy analysis for higher-order RDO schemes is exactly the same as RDO-4. However, a concise Taylor expansion of principal eigenvalue like (
51) and (
52) is very difficult (if not impossible) to obtain even using the mathematical software. Here we try to calculate the truncation order of
of a scheme by the following approximation. First, choose a proper
w, say
. Then the error is
Next, we halve the wave number
w and obtain the error corresponding to
Noticing that
for the
m-th order accurate scheme, one has
Hence the order of accuracy can be approximated by a simple logarithmic operation
The modulus operation is required here since the error result can be complex valued. Errors and estimated accuracy orders for the upwind and the central schemes can be seen in
Table 1 and
Table 2. The errors for the upwind schemes agree with the local discretization error well. The upwind scheme using
M model variables per cell achieves
-th accuracy order strictly. However, the accuracy orders for the central schemes varies quite irregularly. Specifically, when
and 9, the accuracy orders for central schemes remain consistent with that of upwind schemes. When
and 7, central schemes achieve higher accuracy orders than upwind schemes. Nevertheless, there also some choices of
M like
and 8 such that the accuracy order of the central scheme is lower than the corresponding upwind scheme. It is remarkable that when
, the central scheme enjoys two orders higher accuracy than that of the local discretization error.
Table 1 and
Table 2 and also contain abundant information about the dissipation and dispersion of the schemes. It should be noted that a positive real part of
will cause the scheme eventually blow up if no additional dissipation is added in time stepping. In
Table 2, the real parts of
are negative for all upwind schemes which introduces a certain amount of spurious dissipation so as to stabilize the scheme. From
Table 1 we can see that
are pure imaginary numbers with all choices of
M, which means that the dispersion is dominant in the numerical error for the central schemes. However, it can be easily observed that the dispersion and dissipation error decays rapidly as
M increases for the proposed schemes.
In (
39), one can see that the numerical error in computing
dominates the error of temporal derivatives, since
and
have explicit expression. However, the error of computing
is inevitably increasing as
M increases, which indicates that there is a limit of accuracy for the presented method. Denote the numerical inverse of
by
, which is obtained by Matlab. It is well known that the error
hinges on the condition number of
. Moreover, a more accurate estimate
can be bounded by the following inequality,
Then
determines the smallest error that RDO-
M can reach through the grid refinement. The results of
and
with
M varying form 4 to 11 are shown in
Table 3. We can see that
and
grow quickly as
M increases. Hence, there is a tradeoff between the higher convergence rate and
. In this work, we suggest that
is ideal for the practical numerical experiments.
4.3. Stability Analysis
Now we investigate the maximum tolerant time step length when using the Runge–Kutta time integration scheme. Let the time step length be
and
c is the Courant number for (
31). The amplification matrix for the proposed schemes using the fourth order Runge Kutta method for time integration (
19) can be simplified as
To ensure the stability, the Courant number
c should be chosen such that the spectral radius of the amplification matrix
, i.e., the largest absolute value of its eigenvalues, is not greater than 1 for all wavenumbers
w. The spectral radius of
for the upwind schemes and the central schemes are presented in
Figure 2 and
Figure 3, respectively, which can also tell the Courant number limitation related to the stability condition.
Both types of schemes have a fairly large permissible range of
c. As more moments are used, the stability condition becomes more strict for both types of RDO schemes. Generally, the Courant number limitation of the central scheme is larger than that of the upwind schemes. For example, the Courant number limit for upwind RDO-5 is about
and this number is approximately
for the central RDO-5, although both schemes are fourth-order accurate. The Courant number limits for RDO schemes are presented in
Table 4.