New approach to the Lax‐Wendroff modified differential equation for linear and nonlinear advection

The paper presents an enhanced analysis of the Lax‐Wendroff difference scheme—up to the eighth‐order with respect to time and space derivatives—of the modified‐partial differential equation (MDE) of the constant‐wind‐speed advection equation. The modified equation has been so far derived mainly as a fourth‐order equation. The Π‐form of the first differential approximation (differential approximation or equivalent equation) derived by expressing the time derivatives in terms of the space derivatives is used for presenting the MDE. The obtained coefficients at higher order derivatives are analyzed for indications of the character of the dissipative and dispersive errors. The authors included a part of the stencil applied for determining the modified differential equation up to the eighth‐order of the analyzed modified differential equation for the second‐order Lax‐Wendroff scheme. Neither the derived coefficients at the space derivatives of order p ∈ (7 – 8) in the modified differential equation for the Lax‐Wendroff difference scheme nor the results of analyses on the basis of these coefficients of the group velocity, phase shift errors, or dispersive and dissipative features of the scheme have been published. The MDEs for 2 two‐step variants of the Lax‐Wendroff type difference schemes and the MacCormack predictor–corrector scheme (see MacCormack's study) constructed for the scalar hyperbolic conservation laws are also presented in this paper. The analysis of the inviscid Burgers equation solution with the initial condition in a form of a shock wave has been discussed on their basis. The inviscid Burgers equation with the source is also presented. The theory of MDE started to develop after the paper of C. W. Hirt was published in 1968.

predictor-corrector scheme (see MacCormack's study) constructed for the scalar hyperbolic conservation laws are also presented in this paper. The analysis of the inviscid Burgers equation solution with the initial condition in a form of a shock wave has been discussed on their basis. The inviscid Burgers equation with the source is also presented. The theory of MDE started to develop after the paper of C. W. Hirt was published in 1968.

KEYWORDS
Π-form of the first differential approximation, difference scheme dispersive and dissipative features, group velocity, higher order modified differential equation, phase shift errors, Lax-Wendroff method

INTRODUCTION
The explicit difference schemes for the constant-wind-speed advection equation, especially the Lax-Wendroff second-order method, have been commonly dealt with in the recent 50 years-there are many books and papers easily found in the literature concerning the numerical methods for hyperbolic partial differential equations. The history of the MacCormack difference method application in computational fluid dynamics problems may be found in [1]. The Lax and Wendroff [2][3][4], Lax and Richtmyer [5] and all editions of Richtmyer and Morton [6] publications are referred to in each of them. These papers and books are a crucial part of the history of nondispersive and nondissipative numerical methods development. However, the modified partial differential equation of order higher than 6 was derived in none of them.
Warming and Hyett [7] prove that the form of the modified differential equation with the highest order space derivatives is the best scientific basis for the analysis of the difference schemes' features. However, Warming and Hyett, and similarly Duran [8, p. 104], Morton and Mayers [9, p. 172] and Pletcher et al. [10, p. 118], limited the Lax-Wendroff modified differential equation only to the coefficient at the fourth-order spatial operator. In both his monographs LeVeque [11, p. 221] and [12, p. 154] presented the MDE for the Lax-Wendroff scheme only up to third-order.
Authors of this paper found only one publication which contains the modified differential equation for the Lax-Wendroff scheme up to sixth-order. It is Ahmed [13] doctoral dissertation. Unfortunately, the coefficient (5) (we use the symbols introduced by Warming and Hyett in [7, p. 45, eq. (3.34)] at the fifth-order derivative is incorrect. We present its correct form in Section 2, Equation (5). Therefore this paper presents an excerpt of the stencil for determining the coefficient at the fifth-order derivative ( Table 1).
We also used the procedure present in the Warming and Hyett paper [7] and the Pletcher, Tannehill and Anderson [10] and the LeVeque [11,12] monographs to determine the MDE for the Lax-Wendroff difference scheme up to the eighth-order space derivative. We found that for irregular initial conditions (the larger wave-numbers k) it can turn out that the terms of the MDE with higher-order derivatives may play an equally important role.
Shokin et al. [14][15][16][17][18][19] also present a very precise approach to the problem. These scientists call the modified partial differential equation: Γ-and Π-form of the first differential approximation or Γand Π-form of the differential representation of the difference scheme. In the case of hyperbolic type  1 An excerpt of the stencil for determining the modified differential equation up to the eighth-order for the second-order Lax-Wendroff difference scheme (a similar algorithm may also be automated by means of the Maple system-see [13]) equations, the Π-form indicates the Courant number value of the separated diffusive or dispersive difference schemes' features. It also leads to the difference method stability criterion determination. Their monographs [16,17,19] present in the tabular form the stability criteria and examples of the shorted form of the MDEs of the difference schemes, most popular in the middle of the 1980s. This paper was partly prepared to supplement of the Li et al. articles [20][21][22] for the linear advection process and of the Li, Du, et al. [23,24] for nonlinear. Where a theoretical basis for the von Neumann analysis and the modified equation approach for the finite difference schemes are discussed. This paper presents several examples for their theory.
The comeback of the Lax-Wendroff type difference schemes constructed for the scalar hyperbolic conservation laws and characterized by high order time-accurate and high order numerical boundary condition has been observed in a few recent years. Li and Du [23] presented a proposition of the novel two-stage fourth order time-accurate discretization based on the Lax-Wendroff idea. Du and Li [24] derived the fourth order numerical boundary condition which satisfies the hyperbolic conservation laws. So, in our opinion, Section 7 of this paper may be treated as a complement to their papers.
The modified differential equations for the hyperbolic linear and nonlinear problem were developed in the late 1960s by Hirt [25] and Yanenko and Shokin [26]. Next in the early 1970s this study was continued by Shokin [27] and Lerat with Peyret [28][29][30][31] and also Peyret and Taylor [32, section 2.7]. The French scientists called the MDE the equivalent equation.
On the basis of the analysis presented here we come to the conclusions that they may be treated as a supplement of the theory of the MDE for nonlinear problems developed by Lerat and Peyret, especially of their paper [28] published in 1973. We shall discuss it in Section 7.
This paper contains seven sections. We derive the eighth-order modified differential equation for the Lax-Wendroff difference scheme in Section 2. The necessary and sufficient stability conditions are presented in the following section. The dispersion, leading and lagging phase shift errors [11,12,33] and phase and group velocity [34] related with the Lax-Wendroff difference scheme are discussed in Sections 4 and 5. On the basis of the MDE for the Beam-Warming difference scheme (up to the eighth order space derivatives) the equations for the phase shift error and phase and group velocities are also presented in this section. These equations have not been published so far. Section 6 concerns the discussion on the Lax-Wendroff and partly the Beam-Warming difference schemes dissipation.
Section 7 may also be treated as an appendix in which the authors derived in the explicit form the modified differential equations (Γ-and Π-) for the predictor-corrector MacCormack difference scheme [35] and two difference schemes of the Lax-Wendroff type (see also  and Peyret and Taylor [32]) for the scalar hyperbolic conservation laws (inviscid Burgers equation). These derived equations are presented in the divergence form.
In this section we also present the graphs of the nonlinear dissipation and dispersion coefficients and we discuss the new term which appears as an effect of nonlinearity. This analysis enables to minimize the dispersion and dissipation features of the studied methods.

THE LAX-WENDROFF SCHEME MODIFIED DIFFERENTIAL EQUATION
In the case of the homogeneous advection equation the Lax-Wendroff difference scheme (derived in 1960 [2]) is presented below: where h = is the Courant number, h is the spatial step, is the time step, c is the constant-wind-speed. The scheme (2) is second-order accurate ( 2 , h 2 ). Substituting Taylor series expansions for all terms in (2), we can obtain the mixed MDE (so-called Γ-form of the first differential approximation of the difference equation [16,17,19]): Time derivatives of order higher than the first must be removed from (3) by differentiation of the modified partial differential equation to enable discussion of the difference schemes features, for example, the dispersion and diffusion. This approach eliminates the time derivatives, however, it introduces the mixed time and space derivatives which also need removal [10][11][12][14][15][16][17]. Then the modified differential equation's final form (the so-called Π-form of the first differential approximation) for the Lax-Wendroff difference scheme (2) is (using Warming and Hyett's symbols; see [7, p. 166, eq. (3.1)]): where (3), (5) and (7) are the coefficients which are responsible for the dispersion (phase shift errors) and (2), (4), (6) and (8) are the coefficients which are responsible for the dissipation (amplitude errors). The full expansion of these coefficients is presented below: All coefficients (p) (p = 2 -8) were derived using the technique described, for example, in [7,10,14,16,17] (for (5) see Table 1). In the MDE for the Lax-Wendroff difference scheme the coefficient (2) = 0, so the first nonvanishing coefficient at even-order derivative is (4). Usually it is the most important term which is responsible for the dissipative feature of the difference schemes.
Taking the partial differential equations theory and the modified partial differential equation (6) under consideration, we can notice that the truncation error (TE) of the Lax-Wendroff method can be written as:

THE NECESSARY AND SUFFICIENT CONDITIONS OF THE LAX-WENDROFF METHOD STABILITY
The stability condition of the Lax-Wendroff difference scheme is well known. Omitting the transformations, we remind the stability criterion in the form: where k is the wave-number, k = 2 /L, L is the length of the disturbance. So: For the waves which length satisfies the condition kh = , Equation (9) can be rewritten in the following shorter form: Let us refer to the discussion presented in [7]. On its basis, we can derive the necessary and sufficient stability conditions of the two-level explicit difference schemes. The heuristic necessary condition for the stability of these schemes has the form of: where 2r is the lowest even-order nonvanishing coefficient (2r) in the MDE.
For the Lax-Wendroff modified partial differential equation (4) and (5), the even-order nonvanishing coefficient is equal to 4 = 2r, so r = 2. Thus, in the case of the necessary stability condition, the below inequality must be satisfied: Warming and Hyett [7] derived also the sufficient condition of the two-level explicit difference schemes. For r = 2 it has the form: Taking coefficients (5) into consideration, the condition (13) for the Lax-Wendroff modified partial differential equation is: Comparing (12) and (14) we can conclude that in this case the necessary and sufficient conditions are the same and they are equal to 0 ≤ ≤ 1.

THE DISPERSION AND PHASE SHIFT ERRORS
It is well known that the finite difference schemes introduce nonphysical effects to the numerical solutions-usually as dissipation or nonphysical oscillations (called dispersion). It is especially noticed in numerical solutions of the differential problems with initial conditions that include discontinuities or shocks (e.g., a square pulse). Let us split the MDE (6) up to eighth-order derivative: into two parts, named by Durran: prototypical equations. The first, which is responsible for dispersion: x 2p+1 (16) and the second, which is responsible for dissipation: The issue of dissipation will be discussed in Section 6. Comparison of the so called phase shift of the numeric and the phase shift exact e amplification factors may be used to determine the quantity of the dispersion. It enables to compute the relative phase shift error e per one time-increment. For an explicit difference scheme, the phase shift error is a function of the odd-order coefficients of the modified differential equation. It is described by [7]: and the exact phase shift e per one time step is: Considering (5) and (18) one can derive the formula for (k, h, ) for p = 3: So the relative phase shift error 7 / e for the Lax-Wendroff difference scheme takes the form: 6 5,040 The graphs of the function 7 / e (kh) for several selected values of are presented in Figure 1a. This function achieves its minimum for disp = 0.55979 in the range of ∈ (0 -1) and kh = . On the basis of (21) one can say that from the dispersion and irregular initial conditions points of view (e.g., the square pulse) the value of = 0.55979 is the worst choice to do calculations. The graph of 7 / e ( ) for kh = is plotted in Figure 1b. Its minimum is equal to An other important problem is to show how the functions 3 / e , 5 / e and 7 / e defined by (21) in the range of kh ∈ (0 -) for different values of ∈ (0 -1) change separately. In a few papers and books [7,11,12,36] there is information which leads to the conclusion that in the first approximation for small wave-number k (k → 0, so kh → 0) we can neglect all the odd-order terms in the MDE except for the lowest-order nonzero coefficient. In the case of the Lax-Wendroff method this coefficient is (3) and then the phase shift error is equal to 3 ∕ e = 1+( 2 −1) ( ℎ) 2 6 . In general the conclusion mentioned above is true but if we want to provide a detailed analysis, we should also take the coefficients (2p + 1) at the higher-order derivatives under consideration. The expression (21) for 3 / e corresponds with the formula for dispersion of the Lax-Wendroff difference scheme presented by Strikwerda [37].
If we analyze the phase shift error on the basis of 5 we notice that the shortest wavelength (of order 2h) are performed quite well and for = 0.41361 and kh = we obtain 5 = e (Figure 2c-dot-dashed line). It could be meant that for ≥ 0.41361 the Lax-Wendroff difference scheme performs the irregular initial conditions without large deformations. This ambiguity is partly resolved by the graph of the function 7 / e (see Figure 2-dotted line). The graphs of 3 / e , 5 / e and 7 / e defined by Equation (21) are presented in Figure 2. It is clearly seen that at the point ( ℎ) 1 = √ 10 3(6 2 +1) the functions 3 / e and 5 / e start to split and from the point ( ℎ) 2 = √ 10 6 2 +1 they have opposite directions of changes (the function 3 / e constantly decreases and the function 5 / e for √ 10 6 2 +1 < ℎ < increases). The function 7 / e also constantly decreases. Therefore, on the basis of Figure 2 we can say that for the Lax-Wendroff difference scheme the tendency of the successive functions 2p + 1 / e (p = 1, 2, …) changes alternately. It can also prove the correctness of the previous statement that in the first approximation for small wave-numbers we can neglect all the odd-order terms in the modified differential equation but the lowest-order nonzero coefficient [7,11,12,36]. So one can say that − (3)∕c = h 3 6 ( − 3 )∕c and 3 ∕ e = 1 + ( 2 − 1) ( ℎ) 2 6 are the predominant terms which indicate the kind of the phase shift (all subplots in Figure 2-solid lines) and they are amplified by the terms with coefficient (7). The case of the relative phase shift error / e < 1 is called lagging shift and then the Fourier modes will have a lower wave speed than the exact wave speed [7, 10-12, 33, 36]. The opposite case, with relative shift error / e > 1 is called leading shift and then the corresponding Fourier mode will have a wave speed exceeding the wave speed of the exact solution. The Lax-Wendroff difference scheme has only a lagging phase error for 0 < < 1 because for all from this range the phase shifts 3 / e < 1 and 7 / e < 1 even though for ≥ 0.41361 the term 5 / e exceeds 1 for large wave-numbers. For ≥ 0.41361 the inequality ℎ > This difference scheme is unconditionally unstable for > 1, so in the Lax-Wendroff method the leading shift does not appear.
The best example for explaining the differences between leading and lagging shift phase errors is the solution of the advection equation with the periodic boundary conditions and with the regular shape initial condition (sine pulse- Figure 3-top and triangular chapeau pulse-middle). Unfortunately, it is not so obvious for the irregular initial condition (square pulse- Figure 3-bottom). It is clearly shown that in the case of the regular shape initial conditions and / e < 1 the oscillations fall behind the main pulse [11,12]. In both cases (the sine and chapeau pulses), the oscillations amplitudes are rather small. For = 0.55979 and = 0.88487, the values of the first oscillations amplitudes directly at the sine pulse bottom are equal to 0.058 and 0.038, respectively. In the case of the chapeau pulse the values are smaller and they are equal to 0.040 and 0.025.
As it was mentioned above, the explanation of the shift phase errors in the case of the irregular initial condition is not so unambiguous. In the bottom of Figures 3 and 4 (the dot-dashed line) we can observe that the oscillations appear both in front of and behind the discontinuities. We can ask-which oscillations-lower behind or upper in front of-decide about the leading or lagging shift phase errors? It means that the interpretation of the solutions for the strong irregularities is not so clear if we can only see the figure without knowing the values of / e . Furthermore, we can also notice: the main differences between regular and irregular initial conditions are the amplitude values of the oscillations. For the same = 0.55979 and = 0.88487 the irregularities (square pulse initial condition) introduce considerably larger amplitudes of oscillations. They are equal to 0.20 (above the initial value equal to 1) and to 0.153, respectively. Hence they are much bigger than for the sine and chapeau pulses. For much less than = 0.55979 the amplitudes of oscillations are bigger and bigger (e.g., for = 0.287 their magnitude reaches the value of 0.25). In Figure 4 we present the graphs of the solutions of the advection equation computed with the Beam-Warming second-order and the Lax-Wendroff difference schemes. The main differences between them are the locations of the oscillations. For the same value of they there are at the opposite edge of the square pulse. In Section 5 we shall show that they are the consequences of the different features of the both schemes. In the range of ∈ (0 -1) the phase shifts 3 / e and 7 / e for the Beam-Warming methods are equal to or greater than 1 ( Figure 5) while the 3 / e and 7 / e for the Lax-Wendroff method are less or equal to 1. The function 7 / e for the Beam-Warming difference scheme has the form: The graphs of the functions (21) and (22) are plotted in Figure 5. This figure also considers the phase shifts 3 / e and 5 / e for the Crowley difference scheme [38] (Figure 5a,b,d,e). Both methods,  Figure 5 do not present its graphs.
The phase shift 5 / e for the Crowley difference scheme has the form [7,38]: (the coefficient (7) for the Crowley method has not been published so far).
On the basis of Figure 4 we can also state that the amplitudes of the oscillations are a little bigger in the Beam-Warming difference scheme (0.23) than in the Lax-Wendroff method (0.20).
It is worth to notice that for kh = both the phase shifts 5 / e and 7 / e for the Lax-Wendroff (21) as well as for the  difference schemes are equal to 1 for 1 = 1 while the phase shift in the Crowley method ( 5 / e ) Crow for 1 = 1 is greater than 1. For both the Beam-Warming and the Crowley methods the phase shifts ( 5 / e ) Crow = ( 5 / e ) B-W = 1 for 2 = 2 and for this value of 2 the Beam-Warming difference scheme phase shift ( 7 / e ) B-W is also equal to 1.
The coefficient (7) in the modified differential equation for the Lax-Wendroff difference scheme has not been published so far.

THE GROUP VELOCITY OF THE LAX-WENDROFF DIFFERENCE SCHEME
From the very beginning of the numerical methods development it was state that a finite difference method must be stable and accurate [11,12] but various difference methods gave different solutions for the same initial conditions. It was particularly noticeable if the initial conditions were irregular. It soon turned out that the stability and accuracy was not enough. Roache [39] and Mueller [40] proved that the correct transport of any disturbance in the direction of the main flow, so-called transportive property, played equally important role. This feature is characteristic of only the up-wind difference schemes (well known methods of the first-order ( , h) accurate and the second-order ( 2 , h 2 ) accurate-check the Beam-Warming method).
Scientists started to look for new difference methods and new features of the methods already known and used in practice. Group velocity is such a new characteristic. The significance of the group velocity in finite difference schemes is clearly described by Trefethen [34], LeVeque [11,12] and Strikwerda [37].
In this section the equation for the phase and group velocities of the Lax-Wendroff difference scheme is derived on the basis of coefficients (5). We also present the forms of these functions for the Beam-Warming method.
Let us rewrite the MDE for the Lax-Wendroff difference scheme (4) only with the odd-order derivatives and let us assume that this equation presents only its dispersive features [7,8,41]: Equation (24) with (5) has the form: ( Replacing (x, t) by the single Fourier mode we obtain its spectral representation: which directly leads to the phase c ph = (k) k and group c gr = (k) velocities definitions for the Lax-Wendroff difference scheme. On the basis of (26) the phase velocity has the following form: The expansion of (27) with (5) leads to the explicit equation for c ph : 6 5, 040 The phase velocity c ph (28) is equivalent to the relative phase shift error e definition and for the MDE (6) for 2p + 1 = 7 its form is presented below: 6 5, 040 = 7 e (29) The first derivative of the frequency being a function of the wave-number k: is called the group velocity and its form for the other difference schemes can be written as: Its explicit form for the modified differential equation (6) is as follows: 4 24 In Section 4 we discussed the leading and the lagging phase shift errors. Now we can see that the phase velocity c ph /c as well as the group velocity c gr /c for the Lax-Wendroff difference scheme change in the same direction as the relative phase shift errors / e . If 0 < < 1, the values of both the phase and the group velocities are smaller than c (see Figure 6 and Figure 2-dotted lines and also e.g., [10][11][12]42]). In general, the leading and the lagging phase shifts should be associated with the group velocities defined on the basis of the known coefficients of the MDE. The group velocity is the function responsible for moving the oscillations in front of ( for ) the main peak of the solutions.
Analyzing Figure 6 we can understand the differences in the solutions of the linear advection equation computed with the Lax-Wendroff and the Beam-Warming methods (presented in Figure 4). Their best explanation are the graphs of the group velocity functions: Equation (32)-the Lax-Wendroff and Equation (34)-the Beam-Warming plotted next to the solutions of (1) with the initial condition in a sine pulse form (Figure 7).
In the range of ∈ ⟨0 -1⟩ and for kh ∈ ⟨0 -⟩ the group velocity for the Lax-Wendroff difference scheme (32) is less than c (Figure 7a) and then the oscillations move behind the main peak (Figure 7b). For the Beam-Warming difference scheme the group velocity c gr (34) for the same is greater than c (Figure 7c) and the oscillations move in front of the main peak (Figure 7d).
The solid lines in Figure 6 are described below:

THE LAX-WENDROFF METHOD DISSIPATION
In this section we discuss the Lax-Wendroff difference scheme dissipative features. In the modified differential equation (6) there are even-order derivatives with nonzero coefficients. It means that the effects of diffusion or backward diffusion can appear in the solutions. Let us rewrite the equation only with the even-order space derivatives [7,8,41]: and let us assume that the equation derived this way presents only its dissipative features. Equation (35) describes the effects of dissipation (or anti-dissipation, so called backward diffusion or amplification) and with (5) it has the form: ( Replacing again (x, t) by the single Fourier mode: the spectral representation of (35) can be written:  Figure 9b). The opposite case presents the inequality ( (4)k 4 − (6)k 6 + (8)k 8 ) > 0. It follows the exponential growth of the amplitude of the elementary solution (39) and presents the amplification features of the scheme. This fact is partly seen on the graphs presented in Figure 10 by dot-dashed lines illustrating the function Am 6 .
Equation (39) corresponds to results presented by Warming and Hyett [7]. However, the most important difference is that the authors of [7] have not presented the explicit form of the coefficients (2p) for p ∈ (2 -4) for any other difference scheme but the Crowley scheme. In this paper, these forms are also presented for the Lax-Wendroff scheme (see (5)).
Our conclusion here is that it is sufficient to know the lowest even-order nonvanishing coefficient (2r) in the modified differential equation to assess that the Lax-Wendroff difference scheme is strongly dissipative because for 0 ≤ kh ≤ and for the ∈ (0, 1) the coefficient (4) ≤ 0. Due to strong dissipativity, the difference scheme damps all Fourier modes that correspond to nonzero wave-numbers (see also [21, 22, p. 617]). The fourth-order derivatives' coefficients in the MDEs are well-known (see references) but we provide the next two for the even-order derivatives for the Lax-Wendroff scheme. The analysis of the sum Am 6 = ( (4)k 4 − (6)k 6 ) , indicates that for every in the range of (0, 1) and 0 ≤ ℎ < √ 6 this method might be weakly dissipative ( (4)k 4 − (6)k 6 ) < 0 and for √ 6 < ℎ ≤ it is anti-dissipative (( (4)k 4 − (6)k 6 ) > 0). If we plot the amplitude e Am 6 it occurs that for √ 6 < ℎ ≤ the Lax-Wendroff difference scheme amplifies nonphysical solutions. The curves in Figure 9b increase instead of decreasing. This anti-dissipative feature indicates that for these values of kh and ∈ (0, 1) this difference scheme is unconditionally unstable.
However, such analyses result in wrong conclusions. In the modified differential equation, the lowest even-order nonvanishing coefficient (2r) has the strongest influence on the dissipative features of the difference scheme. Here it is (4). So the Lax-Wendroff difference scheme is strongly dissipative of order 4 (see also Strikwerda [37] and Trangenstein [42]).
The Am minimum value of Equation (41) is equal to −11.24 for kh = and for diss = √ 2∕2 (Figure 9a). So we can say that from the dissipative point of view this lambda value is the worst choice to do calculations (see Figure 8b). Moreover, the coefficient (8) has the greatest influence on the Lax-Wendroff difference scheme's dissipative features. Our analysis also leads to a wider conclusion that (8) amplifies the influence of the coefficient (4) rather than decides about the dissipative features of the Lax-Wendroff difference scheme. Therefore, on the basis of MDE (6) one can agree that the coefficient (8) is a supplementary coefficient of the dissipation process explanation and (4) remains the leading coefficient.
The graphs of Am presented in Figure 10 (dotted line) show the magnitude of the Lax-Wendroff method dissipation and explain the numerical solutions of the constant-wind-speed advection equation for k = 126- Figure 12 (subplots to the right).
The coefficients (6) and (8) in the modified differential equation (see Equation (6) for the Lax-Wendroff difference scheme have not been published so far.
The exponent Am for the Beam-Warming second order difference scheme has the following form: 8 It has two extrema (minima) for 1 = 1 − √ 2∕2 and 2 = 1 + √ 2∕2 and they are symmetrically placed relative to a straight line for = 1. The graphs of Am B -W for selected are presented in Figure 11. On the basis of the subplot in Figure 11b we can see that for = 0.5 the exponents Am and Am B -W for both the Lax-Wendroff and the Beam-Warming methods are equal. For < 0.5 the damping of the Beam-Warming scheme increases and for = 1 − √ 2∕2 achieves its maximum (see Figure 11a). It is the worst value of to do calculations on the basis of this scheme. For → 1 the exponent Am B -W rapidly achieves 0 and the amplitude of e Am B-W does not change in time and the damping of the solutions is equal to zero.
For > 0.5 the Lax-Wendroff method dissipative feature starts to grow and for = √ 2∕2 achieves its maximum (see Figure 11a). In the range of ∈ (0.5, 1.0) this feature slowly gets weaker, the exponent Am approaches zero and amplitude of e Am does not change in time. For = 1 it is equal to one.
The solutions of the linear advection equation computed by the Lax-Wendroff difference scheme for various wave-numbers (k = 9, k = 25 and k = 126) and for different types of initial conditions: sine pulse, triangular chapeau pulse and square pulse are presented in Figure 12. The successive three values of k define various kinds of initial disturbances. In our calculations, the spatial step h = 2 N x = 0.024933, so for k = 126 we obtain kh = and k = 126 is the largest wave-number that, in this case, can be solved by the finite difference schemes. For k = 9 the initial disturbance spreads on 29 spatial points, for k = 25 on 11 points and for k = 126 only on 3 spatial points.
The damping of the irregularities is the biggest for the wave-number k = 126 but the range of spreading nonphysical oscillations is much the same for other wave-numbers k (Figure 12). Analyzing all subplots in Figure 12 we can understand the wave-number k value influence on the differences in the solutions.
For large values of k (k = 126) the damping of the solutions for all initial condition types is almost the same.
For middle values of wave-number k (k = 25) the damping clearly depends on the initial condition type-it is the biggest for the chapeau type pulse.
For small values of wave-number k (k = 9) and regular initial condition (sine pulse) the damping of the solution is almost equal to zero. In the case of the chapeau pulse the damping starts to develop. For irregular initial conditions (square pulse) the damping of the solution is the smallest but the magnitude of the oscillations is the biggest. Figure 13 presents the graphics of the solutions of the advection Equation (1) with irregular initial condition (square type pulse) for k = 126 and for selected (Top). The initial  Figure 13 by the envelope of the peaks max .

THE ANALYSIS FOR ONE-DIMENSIONAL SCALAR CONSERVATION LAWS
Let us consider the scalar hyperbolic conservation laws in the following form: where is a conserved variable, and f( ) is a known flux function. For 2 2 Equation (43) becomes the inviscid Burgers equation.
Applying two Lax-Wendroff type difference schemes and the MacCormack method (McC) we shall solve this equation with the initial condition in the form of the right-moving discontinuity: But at first we shall derive their modified partial differential equations and next on their basis we shall discuss the differences between the obtained solutions.
Let us start our analysis from the method which is a modification of the most popular two-step variant of the L-W type scheme (LW-1): where = /h. The mixed (with time and space derivatives) MDE for the LW-1 method (45) is as follows: Almost all terms in Equation (46) depend on the space derivatives. But one-in box-is the second order time derivative. The construction of the modified differential equation requires to eliminate this derivative (see also Section 2). Several mixed derivatives, for example, tx 2x , x t2x , t 2x and so on and higher order space derivatives (p > 3) are hidden in the TE (h 3 ). The underlined terms play the most important role. The successive steps of the transformations introduce a new term, 3 3x , which is also an important part of the Π-form of the first differential approximation for the Lax-Wendroff difference scheme (45).
After several complicated transformations we can present the final form of the modified differential equation for LW-1 (46): (47) When we gather the coefficients of Equation (47) at the diffusion 2 x 2 and dispersion 3 x 3 parts we obtain the simpler MDE and we can describe its form explicitly: We have a few new expressions in the brackets of (48). The first is a nonlinear numerical dissipation with coefficient − 1 4 (5 2 2 − 1) x . The second is a nonlinear numerical dispersion with coefficient − 6 (− 2 2 − 1). And the third, − 3 , is a function of the first derivative values on the shock wave slope and it can be treated as values of the shear on the slope. Let us call it the shear.
The analysis of these coefficients will be presented for the two difference schemes of the Lax-Wendroff type and for the McC.
The terms of the modified differential equation (47) may come down to identities: which enable us to present the final form of the MDE (47) in the scalar hyperbolic conservation laws: (50) where N is the nonlinear operator for MDE of (43) in the conservative form.
The next difference scheme of the Lax-Wendroff type has the form of LW-2: Its mixed modified differential equation is presented below: Omitting the indirect transformations we can present its final form: (53) The second and the third expressions in brackets can be transformed to the conservative form: These identities are noticeably less complicated than those in (49) for the Lax-Wendroff LW-1 difference scheme (45). The full Π-form of the MDE for the LW-2 difference scheme in the conservation laws can be written as: Expansion of Equation (53) into a form similar to (48) leads to: The two-step difference scheme, well-known as the McC, was for the first time introduced to the computational fluid dynamics by MacCormack [35]: which leads to the mixed modified differential equation: and its final Π-form of the first differential approximation: After several transformations we obtain the modified differential equation for the MacCormack difference scheme in the scalar hyperbolic conservation laws: Gathering the coefficients at the second and third derivatives in (59) leads to: (61) On the basis of this discussion one can collect all the coefficients at the terms of the MDE for future analysis of the presented difference methods ( Table 2). These coefficients enable us to explain the differences in the observed amplitudes values of the oscillations and their range of spreading.
It is important to state that the coefficients presented in Table 2 were published for the first time by Lerat and Peyret [28, p. 760]. In their paper the authors described in a coherent way the influence of E 1 , E ′ 2 and E 3 on the solutions of (43) with the initial condition (44) for two Lax-Wendroff difference schemes and two McCs. In this paper we derive step by step all forms of the MDEs (Γ-and Π-), their coefficients and we present them in the divergence forms. In the next section we shall also show the graphs of the coefficients E 1 and E ′ 2 separately, estimate the values of the function E 3 ( x ) 3 and Note: They describe three processes: dispersion, dissipation and nonlinearity.
present graphically the sum of the terms in the brackets in the equation over Table 2. Thus this section may also be treated as an appendix and expansion of the Lerat and Peyret paper [28].

The analysis of the numerical test
Knowing the general form of the modified differential equations: let us solve the inviscid Burgers equation (43) with the initial condition (44) and estimate the influence of the specific terms of Equation (62) on the solutions. In Figure 14 we present 10 graphs for = 0.1 + nΔ (n = ⟨1 -9⟩) and Δ = 0.1. They were computed for L = 1, R = 0. On the basis of Figure 14 we can state that nonphysical oscillations appear in all solutions. Their amplitudes and range of spreading are different-the biggest occur for the LW-1 scheme and the smallest for the MacCormack scheme.
Let us notice that the coefficients E 1 ( , ) = 6 ( 2 2 − 1) at the third space derivatives are the same in all derived modified differential equations (see (48), (56), (61) and Table 2) although they are different in the mixed forms of their first differential approximations (46), (52) and (58): (63) Besides, Lerat and Peyret showed in [28] that the other difference scheme of the MacCormack type (not presented here) characterizes the same coefficient at the dispersion term.
Taking these facts under consideration we can conclude that both methods of the Lax-Wendroff type and the MacCormack scheme introduce into solutions the same values of the dispersion equal to −E 1 ( , )   The presented analysis leads also to the conclusion that independently of the different space-steps h one can find such a value of that the coefficient of the dispersion E 1 ( , ) = 0. Then the oscillations are evidently smaller but in the Lax-Wendroff methods they still do not disappear. It can mean that the differences in amplitudes values shown in Figure 14 Section 6) or by the third term which is responsible for the shear (nonlinearity) 3 or by both these effects simultaneously. Figure 15b presents the graphs of the functions being the coefficients of −E ′ 2 ( , ) at the second space derivatives in the derived modified differential equations. Their analysis can help us to explain the differences in the values of the oscillations amplitudes.
For several initial s all coefficients −E ′ 2 ( , ) are negative, so they describe typical cases of the dissipation (see the equation at the top of Table 2). For greater s the coefficients −E ′ 2 ( , ) in both Lax-Wendroff methods become positive and the origin of the oscillations is mainly related to the anti-dissipative features of these methods and in our opinion it is the main reason that the oscillations remain.
The smallest damping of the oscillations characterizes the LW-1 difference scheme: −E ′ 2 ( , ) = −0.237 for = 0.1. It is noticeable in Figure 14 Figure 15b). It means that for > 0.4 the LW-1 scheme amplifies the unwanted numerical effects and becomes the anti-dissipative scheme. For = 1 the basic dispersion is equal to zero, E 1 ( , ) = 0, but the oscillations still develop, though their amplitudes are evidently smaller than those for = 0.1 and equal to max = 1.127 (see Figure 14 on right).
The solutions computed by the LW-2 difference scheme are placed between the results of LW-1 and the McC methods (see Figure 14-middle). For = 0.1 the LW-2 scheme damping is equal to −0.485 and it is greater than in the LW-1 method (see Figure 15b). The oscillations amplitudes in LW-2 scheme are nearly the same (1.397) as in the McC (1.366). This method starts to introduce into solutions the anti-dissipative features for > 0.6 (n = 5). This conclusion is graphically confirmed in Figures 14 and 15b (green curve).
From the point of view the damping and amplifying of the nonphysical effects, the MacCormack scheme is the best one. Admittedly this method does not become anti-dissipative in the whole range of ∈ ⟨0 -1⟩. For every this method strongly damps the initial oscillations and for = 1 it hardly amplifies its dispersive and dissipative features. On the basis of Figure 14 (bottom) one can conclude that for > 0.8 the MacCormack difference scheme is neither dispersive nor dissipative (the oscillations disappear).
For = 1 the coefficients −E ′ 2 ( , ) = 1 for both LW-1 and LW-2 methods (Figure 15b-green and brown lines). Moreover, the coefficient of the dispersion −E 1 ( , ) = 0. So one can ask: why these methods with the same coefficients values of the dispersion and dissipation in the modified differential equations compute such different solutions?
In our opinion the explanation of the above described effects should be connected with the non- 3 . In Table 2  should not be interpreted only in the category: dispersion-dissipation. Its interpretation should be a little wider, for example: dispersion-dissipation-nonlinearity or: dispersion-dissipation and at the right-hand side-the source. Because this term does not appear in the MDE of the linear advection Equation (6), one can say that it is the effect of the presence of the conservative term ( 2 ∕2) x in Equation (43). The graphs of the terms which are responsible for the nonlinearity are drawn in Figure 16. For both Lax-Wendroff schemes they are always positive and for → 1 their values are comparable with the Comparing the subplots Figures 17b and 16 we come to the conclusion that the share of the nonlinearity in the solutions computed by means of the analyzed difference schemes is evidently significant.
LeVeque [14, section 2.5 and Chapter 17] considers the nonhomogeneous conservation law (also often called the balance law) in the form: where ( ), generally called the source, represents a sink if ( ) < 0 or a source if ( ) > 0. In the case of the modified differential equations derived here, the right-hand side of (64) should rather be treated as a sink. Let us rewrite the MDEs (62) to the form: Analyzing Figure 17c one can conclude that for LW-1 and LW-2 the right-hand side of Equation (65): is always negative for all s. In the case of the McC difference scheme it is negative for > 0.5. In our opinion the role of this function is important. Its property weakens the anti-dissipative oscillations. The explanation of this statement is seen in Figure 17b,c. For = 1 both functions −E ′ 2 ( , ) 2 x 2 in the MDEs of the Lax-Wendroff difference schemes are equal to 1. But the oscillations in the solutions computed by the LW-2 scheme are smaller than those in LW-1 (see Figure 14). It means that the damping by the nonlinearity term is greater in the LW-2 scheme than in the LW-1 (see Figure 17c-green line).
For the assumption that ( ) = h 2 E 3 ( , ) ( x ) 3 the nonlinearity term cannot be neglected in such analysis. It confirms our previous conclusions.

THE FINAL CONCLUSIONS
1 The Lax-Wendroff method is dissipative over the intervals of ∈ (0, 1). It is never anti-dissipative. Thus, it means that this difference scheme does not amplify nonphysical solutions. 2 In the modified partial differential equation (6) of the Lax-Wendroff method, no numerical dissipation is described by the second-order space derivative 2 x 2 , typical of one-sided finite-difference approximation. 3 In case when approaches 0, the method's properties are determined by anti-dispersion and dissipation (see Figures 2 and 8). This conclusion is satisfied for kh > (kh) 0 . The value (kh) 0 must be experimentally computed for every on the basis of assumption that for kh = (kh) 0 the relation 7 e ( ℎ) = 0 is satisfied. It means that for kh > (kh) 0 the function 7 e < 0. It also means that the short-wave disturbances (large-wave-number k) are reproduced worse and worse (see Figure 13). 4 In case when approaches 1, the Lax-Wendroff method's properties are determined by dissipation and dispersion (see Figures 2 and 8). If changes in the range (0.88458 -1) this method is dispersive for all kh ∈ (0, ) and if changes in the range (0 -0.88458) it is dispersive only for kh > (kh) 0 . The short-wave disturbances are also poorly reproduced (see also Figure 13 for wave-number k = 126).