Quasi-linear analysis of dispersion relation preservation for nonlinear schemes

In numerical simulations of complex flows with discontinuities, it is necessary to use nonlinear schemes. The spectrum of the scheme used has a significant impact on the resolution and stability of the computation. Based on the approximate dispersion relation method, we combine the corresponding spectral property with the dispersion relation preservation proposed by De and Eswaran (J Comput Phys 218:398–416, 2006) and propose a quasi-linear dispersion relation preservation (QL-GRP) analysis method, through which the group velocity of the nonlinear scheme can be determined. In particular, we derive the group velocity property when a high-order Runge–Kutta scheme is used and compare the performance of different time schemes with QL-GRP. The rationality of the QL-GRP method is verified by a numerical simulation and the discrete Fourier transform method. To further evaluate the performance of a nonlinear scheme in finding the group velocity, new hyperbolic equations are designed. The validity of QL-GRP and the group velocity preservation of several schemes are investigated using two examples of the equation for one-dimensional wave propagation and the new hyperbolic equations. The results show that the QL-GRP method integrated with high-order time schemes can determine the group velocity for nonlinear schemes and evaluate their performance reasonably and efficiently.


Introduction
Nonlinear approaches, such as the weighted essentially non-oscillatory (WENO) scheme, are widely used to calculate complex flow fields with discontinuities, such as shock waves.Because of the mathematics of such schemes, truncation errors can affect the convergence rate of the computation and the spectral properties affect the deviation of the Fourier modes of the numerical results from the exact solution.Therefore, studying the spectral properties of nonlinear schemes is important and meaningful when developing high-order nonlinear schemes.
Tam and Webb [1] suggested that it was necessary to consider the spectral properties besides using the standard Taylor series method in the construction of difference schemes.Specifically, by using the Fourier transformation of the governing equation, they proposed a basic method for analyzing the dispersion and dissipation relation of a linear scheme and developed a dispersion relation preservation (DRP) scheme.
Moreover, this method can also be used to evaluate the spectral properties of different linear schemes and to construct new schemes.However, the method is inappropriate for analyzing the spectral properties of nonlinear schemes.To remedy this deficiency, Pirozzoli [2] studied the amplitude evolution of a disturbing mode after a tiny propagation period through a discrete Fourier transform (DFT), through which he proposed an approximate dispersion relation (ADR) analysis.ADR [2] gives an estimate of the total error generated by a nonlinear scheme and provides a method for analyzing its spectral properties.However, there are still some issues.For example, the time scheme can introduce errors, and if the grid number is not chosen carefully, there can be unexpected jumps in the spectral distribution [3] .Observing these problems, Mao et al. [3] proposed an improved ADR without time discretization called ADR-NT.Using the assumption of a tiny time period, they transferred the effect of the temporal derivative into the contribution of the spatial derivative.The researchers reported [3] that ADR-NT not only reduced the computational cost but also avoided the error due to the time discretization, as occurs in ADR.Mao et al. [3] also suggested ways to avoid jumping points when using ADR, i.e., the size of the grid should be twice some large prime or take the average over a large number of results with different initial phases.
The dispersion relation usually reflects the phase velocity of a spatial scheme.For practical problems with wave propagation, it is necessary to consider the group velocity, which is the one the energy propagates at [5] .
De and Eswaran [4] analyzed DRP from the perspective of the group velocity, which we refer to as the DRP-g method or just DRP-g.They pointed out that DRP-g was important when evaluating the spectral properties of linear schemes.In the space of the reduced wave number and frequency, group velocity preservation (GVP) occurs in the region where the ratio of the numerical group velocity to the theoretical group velocity is in [0.95, 1.05].Obviously, the analyses of DRP-g proposed by De and Eswaran [4] are different from those of Tam and Webb [1] and Pirozzoli [2] .However, the DRP-g method can be used only for linear schemes and not for nonlinear schemes.Thus, we derived the modified wave number of nonlinear schemes with the ADR method, combined the aforementioned group velocity analysis with high-order Runge-Kutta schemes, and finally, derived a quasi-linear method to analyze the spectral properties of nonlinear schemes.For brevity, the method is called the QL-DRP method or just QL-DRP.Its rationality was then verified with a numerical test and the DFT method.
As usual, a one-dimensional wave propagation equation [4] is used in this study to explore the GVP of nonlinear schemes by QL-DRP.As shown in [4], for this governing equation, the group velocity is the same as the phase velocity.Therefore, it is insufficient to investigate the GVP of difference schemes.To overcome this, we devised hyperbolic equations such that the solutions are synthetic waves where the group velocity and phase velocity are different.Using the equations, we numerically analyzed the GVP of different difference schemes and also verified the validity of QL-DRP.This paper is arranged as follows.The DRP and ADR methods are reviewed in Section 2. In Section 3, the DRP-g formula is derived for a high-order Runge-Kutta scheme, the QL-DRP method is described, and the comparative results from QL-DRP are given for typical schemes.In Section 4, the group velocities of selected cases are obtained computationally using DFT and compared with those from QL-DRP, through which the rationality of QL-DRP is verified.We construct hyperbolic equations for which the solutions are waves where the group velocity and phase velocity are different.In Section 5, the one-dimensional wave propagation equation and the aforementioned hyperbolic equations are solved and analyzed.The conclusions are drawn in Section 6.

Review of DRP-g and ADR methods
Before we consider the QL-DRP method, the DRP-g and ADR methods are reviewed first.

DRP-g method
Consider the one-dimensional wave propagation equation [1][2][3][4] : For the initial distribution When using a time scheme, the following approximate relation is obtained: where n u is the value at the initial time When the explicit Euler scheme is used, Eq. ( 2) can be written as: where   , n n f t u denotes the contributions from the spatial derivative.When a spatial scheme is used, one can see [4] that , where  represents the modified wave number.It can be further derived from Eq. (3) that: Following the approach of De and Eswaran [4] , the numerical value , where t  is the reduced frequency, kx   is the reduced wave number, and c is the theoretical group/phase velocity.Equation ( 4) is exactly the same as that proposed by De and Eswaran [4] .
When a linear scheme is employed, the spatial derivative at i x can be approximated by: 1 () , where M and N are the numbers of nodes on the right and left sides of i x , and j a are the scheme coefficients.
The modified wave number can be obtained by a Fourier transformation [1] : Once the dispersion and dissipation relations are obtained, , g num V for the first-order explicit Euler scheme can be obtained from Eq. (4).The explicit Euler scheme yields relatively large temporal errors, and therefore, the group velocities obtained by DRP-g are inappropriate for unsteady problems when high-order Runge-Kutta schemes are used.Thus, Sengupta et al. [7]     is the Courant-Friedrichs-Lewy (CFL) number.On checking, the following problems were found with this work: 1.There are no details for the derivation.Therefore, it is difficult to check the correctness of the formula.
2. When  approaches 0, the formula does not regress to that of the first-order explicit Euler scheme in Eq. (4).

 
, the denominator cos( ) t  is 0, and there is a singularity.Therefore, the DRP-g formula needs further investigation for higher-order Runge-Kutta schemes.Thus, we employ the idea of [4] and derive the DRP-g results for high-order Runge-Kutta schemes in Section 3.1.

ADR method
Despite its applicability to linear schemes, the spectral analysis method by Tam and Webb [1] is unsuitable for nonlinear schemes.Thus, Pirozzoli et al. [2] proposed a quasi-linear spectral analysis method for nonlinear schemes or ADR.The process is as follows [2] .
For Eq. ( 1), suppose the initial distribution is 0 ( ,0 and consider a semi-discretized approximation of Eq. ( 1) on grids    with the spacing x  as 0 0, (0) where  is for a specific numerical scheme, e.g., the linear difference scheme: , as before.From [2], Eq. ( 6) has the following exact solution: . Furthermore, Pirozzoli [2] proposed ADR to derive  : 0 ( , ) ( ) ln () In this equation, 0 () n v  is the Fourier coefficient corresponding to the initial mode and ( , ) where x N is the number of grid number.
As indicated in the introduction, ADR depends on the choice of time step and number, and inappropriate choices can lead to temporal errors or unreasonable jump points in the spectral distributions [3] .Mao et al. [3] expanded the solution   j v  of Eq. ( 6) with a Taylor series: Substituting the above formula into Eq.( 8) and the subsequent result into Eq.( 7), Mao et al. [3] proposed ADR without time discretization, known as ADR-NT, to derive () where jl b  are the coefficients of the nonlinear scheme.Compared with ADR, ADR-NT reduces the computation cost and minimizes the influence of temporal errors.

QL-DRP for nonlinear schemes
De and Eswaran [4] first analyzed and derived the numerical group velocity for linear schemes.However, their method can be used only for linear schemes and not for nonlinear ones.Considering the development of ADR [2] from DRP [1] , which extends the analysis from linear to nonlinear, we introduce the influence of the nonlinearity in DRP-g [4] by using ADR to account for the spectral property of the nonlinear scheme.Next, employing a similar procedure as in DRP-g, the so-called QL-DRP is derived for nonlinear schemes.Prior to a further discussion on QL-DRP, we first derive the DRP-g formula for high-order Runge-Kutta schemes.

DRP-g formula for high-order Runge-Kutta schemes
Taking the third-order Runge-Kutta scheme as an example, we derive the DRP-g formula based on the method of [4].For convenience, the nth-order Runge-Kutta scheme is referred to as n RK .
Applying (.) / d dk to the above equation and taking the real part as in [4], the group velocity for RK3 is For a linear scheme,  is found from Eq. ( 5) [4] .Similarly, when an 4 RK scheme is used for time discretization, the DRP-g formula can be derived as: From Eqs. ( 13) and ( 14), it can be seen that the numerical group velocity for a high-order Runge-Kutta scheme depends on  , unlike the first-order explicit Euler scheme in Eq. ( 4).For illustration, the group velocity of a linear fifth-order scheme, UPW5, namely: () , as mentioned before.As indicated by De and Eswaran [4] , the group distribution near the origin is of special concern.Figure 1 shows that when the CFL number is small, the distributions of the first-order explicit Euler scheme, RK3, and Euler scheme, which indicates that higher-order time schemes enhance the GVP.This outcome confirms the potential advantages of higher-order time schemes for unsteady problems.
In summary, in this section the group velocities of RK3 and 4 RK were found by applying the same procedure used for the explicit Euler scheme in [4].The following remarks are given: (1) When the CFL number approaches 0, the group velocity distributions of higher-order Runge-Kutta schemes gradually degenerate into that of the first-order explicit Euler scheme.Further, the results from DRP-g for different time schemes become the same under this condition.(2) No singularity arises in Eqs. ( 13) and ( 14), which indicates their rationality.

QL-DRP
As mentioned in the introduction, QL-DRP is proposed to evaluate the characteristics of the group velocity in a quasi-linear manner.The main difference between QL-DRP and DRP-g is that the modified wave number in the former considers the nonlinearity of schemes.Therefore, the formulas for QL-DRP are the same as those for DRP-g, such as Eqs.( 13) or ( 14), except that different evaluations of  are implemented.
The concrete implementations are summarized as follows: 1. Solve for the modified wave number ()   of nonlinear spatial schemes.
Based on ADR, the dispersion and dissipation relations of the nonlinear scheme, ()   , are solved with Eq. ( 9).
When the spatial scheme is linear, ()   can be obtained as () . Because the coefficient j a is constant, dd   can be derived analytically.However, when the spatial scheme is nonlinear, the coefficient jl b  in Eq. ( 9) correlates with the initial variable distribution nonlinearly and dd   is hard to resolve analytically.To overcome this difficulty, we use the difference method to evaluate it numerically.
Taking the second-order central difference as an example, dd   can be evaluated as: where j is the index of the reduced wave number and When solving for ()   for nonlinear schemes with ADR, t  in the time integration must be small enough [2] , which means the corresponding reduced frequency t   should also be small.Correspondingly, QL-DRP can describe the GVP of nonlinear schemes effectively only for a small reduced frequency.When t   is large, ADR does not accurately predict the spectral property, and therefore, QL-DRP cannot accurately describe GVP either.However, QL-DRP can give some insights into the group velocity of a nonlinear scheme.
Hence, QL-DRP can provide a useful reference for the overall characteristics of the group velocity for a nonlinear scheme.
For demonstration, we used QL-DRP to find the group velocity for three typical spatial schemes: UPW5, WENO5-JS [8] , and WENO5-M [6] .Adopting the aforementioned procedure, the corresponding distributions of  The figure shows that: 1. Comparing Fig. 3(a) with Fig. 1, one can see that the group velocity distributions for QL-DRP with the UPW5 scheme are the same as those for DRP-g.This is because, for linear schemes, ()   in ADR is the same as that in DRP [1] .In other words, QL-DRP degenerates into DRP-g for a linear spatial scheme.
2. The distributions of the group velocity are obviously different for the three schemes, which indicates that QL-DRP can be used to distinguish between the GVP of different spatial schemes under the same time scheme.Specifically, the area of the GVP region (green) for the WENO5-JS scheme is smaller than that for the UPW5 scheme, which indicates that the group velocities of the nonlinear scheme differ from that of scheme.Moreover, the GVP area of the WENO5-M scheme is larger than that of WENO5-JS, which shows that the nonlinear optimization of the former enhances the GVP.
3. The GVP area near the origin, especially along the vertical axis, is the largest for UPW5, followed by WENO5-M and then WENO5-JS.
In short, QL-DRP can be applied not only to analyze the group velocity of a linear scheme but also of a nonlinear scheme.Moreover, the method can clearly differentiate spatial schemes by considering the size of the GVP region.

Numerical analysis and validation of QL-DRP
The previous section introduced the derivation and implementation of QL-DRP.In this section, the numerical group velocity is obtained numerically and the results are compared with those from QL-DRP.Thus, the rationality of the method was analyzed and verified.Moreover, we devise hyperbolic equations and a corresponding example to determine the numerical group velocity of the scheme directly, which further verifies the approach.

Numerical analysis and validation
In this section, the group velocity of the different schemes is numerically obtained by solving the onedimensional wave propagation equation and by using DFT.By comparing with the analytic results from QL-DRP, the numerical error is acquired and the rationality of QL-DRP is evaluated.Next, the derivation of numerical group velocity will be discussed.The solution of the semi-discrete equation ( 6) is where ()      [10]   .Compared with the exact solution , the numerical frequency can be expressed as / cx    , which depends on  .Hence, the reduced frequency at Pi can be expressed with the modified wave number as: where L is the length of the computational domain.We consider Eq. ( 1) with 1 c  and the initial distribution . From [2],  can be numerically found with ADR [2] for a spatial scheme and time scheme using DFT.Here, we use WENO5-JS and RK4 for the spatial and time discretization.
The number of grid cells and time step are given in Table 1.As shown in the table, four cases with respective (Nx, t  ) are chosen, and these correspond to Pi.Then,  can be acquired with respect to  , as illustrated in Fig. 4.Besides P1 to P4, we also investigated the case when the reduced frequency was larger.There were oscillations in the results.Hence, the frequency region for investigation was limited to  Thus far, the numerical solution of the group velocity can be found from: Re( ) Re( ) ( ) where k1 and k2 represent two wave numbers with a small difference, 1  and 2  are the corresponding numerical frequencies, and Re denotes taking the real part.Substituting Eq. ( 16) into Eq.( 17): Re( ( ) ( )) Under the assumption that at Pi can be obtained from Eq. ( 18).These are compared in Table 2. Table 2 shows that when t  is small, as for P1, P2, and P3, then the group velocity given by QL-DRP is quite close to the numerical value.This indicates that QL-DRP can reasonably predict the group velocity at low reduced frequency.With an increase of t  , the error between , g anal V and , g num V becomes large but was still within 21% when 1 t  .Therefore, QL-DRP can provide important insights for the group velocity at medium reduced frequency.
Note that QL-DRP can provide an overview of the group velocity in the kx  vs. t  space, which is difficult for the numerical approach just mentioned.In addition, QL-DRP avoids the tedious process of finding a numerical solution and is more convenient for practical analyses.  , respectively.One can see that in such situations, the group velocity and phase velocity are indistinguishable, which is unfavorable for investigating the GVP property of the numerical scheme.
In textbooks, the combination wave is usually used to explain the concept of group velocity, which can also be written as: The envelope of Eq. ( 19) is Although the combination wave is used to explain the group velocity and phase velocity, the distribution does satisfy Eq. ( 1).Therefore, this equation cannot be used to simulate the wave.We devised the following hyperbolic equations so that we could numerically study the GVP of different schemes:

Numerical examples
In the following, to compute the group velocity, the envelope of the wave is derived by a Hilbert transform.

One-dimensional wave propagation
0.8259, 0.9592, and 0.9950, respectively.These values are nearly the same as those obtained numerically.
Hence, QL-DRP has been verified quantitatively, and its usefulness in analyses of the group velocity has been demonstrated.

Wave propagation with different group and phase velocities
In this section, the case described in Section 4.2 with different group and phase velocities is evaluated for the different spatial schemes, namely UPW5, WENO5-M, and WENO5-JS, which further demonstrates the validity QL-DRP.
For the distribution in Eq. ( 19), f multiple waves were included in one wave packet in purpose of better illustration.The initial condition is in the domain [ 3 ,3 ]   .We set 11 6 k   , 2 8 k  , and 2 12   .The exact group and phase velocities were The corresponding envelopes, which contain information about the group velocity, were derived, as shown in Fig. 9.By convention, the envelopes were drawn after taking the absolute value.Figure 9 shows that UPW5 yields an envelope such that the numerical group velocity is closest to the theoretical one, followed by WENO5-M and then WENO5-JS.The amplitudes of the envelope show that was investigated for different CFL numbers.The velocity distributions are shown in Figs. 1 and 2 for 0.01   and 0.1, respectively.The effect of  is apparent through the variation of the GVP region, for which

/ g Vc with the 4 RK
scheme are shown in Fig.3.

First, we select
several points Pi in the kx  vs. t  plane, then numerically compute the group velocity using DFT.The coordinates of Pi are the reduced frequency and wave number.Without loss of generality, we set kx  as 1.Because kc are determined, the coordinates of Pi are defined.

Fig. 4 .
Fig. 4. Illustrations of Pi on the distributions of Re( )  for different sizes of grid and different time steps with WENO5-JS and 4 RK .

1 t
 .For illustration, the positions of Pi (i = 1, 2, 3, 4) are shown in the kx  vs. t  plane in Fig. 5.The colors indicate the group velocity calculated by QL-DRP with WENO5-JS and RK4.

Fig. 5 .
Fig. 5. Positions of Pi in the kx  vs. t  plane.The colors indicate the group velocity calculated by QL-DRP with the group velocity is different from the phase velocity.Hence, a measure is provided to study the combination wave with the cooccurrence of different group and phase velocities, which also provides a powerful way to verify the outcome of QL-DRP.
are shown in Fig.7.

Fig. 9 .
Fig. 9. Envelopes for T = 1 for the hyperbolic equations under the conditions Nx = 120 and assumed that 44

Table 1 .
Size of computational grid, time step, and coordinates of Pi.

Table 2 .
Both group velocities and errors of Pi for WENO5-JS and RK4.