NONLINEAR MONOTONIZATION OF THE BABENKO SCHEME FOR THE QUASI-LINEAR ADVECTION EQUATION 1

. The paper is devoted to construction and development of new method for numerical solution of hyperbolic type equations [14, 17]. In the previous papers [4, 5, 6, 7, 8, 9] authors have investigated theoretically and tested experimentally 26 different ﬁnite-difference schemes on 4 point patterns for the simplest hyperbolic equation: linear advection equation. This equation has the main features of every hyperbolic equation and is the important part of many mathematical models. In other cases the advection operator is the important part of the full operator of the problem. All 26 schemes have been compared experimentally on the special representative set of tests. Nevertheless to simplicity of the equation, almost all schemes have different disadvantages. They are dis-cussed in detail in the cited papers. So, the investigation of new schemes for this equation is still an important task. monotone Babenko ("square") generalize to more difﬁcult The a quasi–linear advection

Received October 21, 2004; revised March 4, 2005 Abstract. The paper is devoted to construction and development of new method for numerical solution of hyperbolic type equations [14,17].
In the previous papers [4,5,6,7,8,9] authors have investigated theoretically and tested experimentally 26 different finite-difference schemes on 4 point patterns for the simplest hyperbolic equation: linear advection equation. This equation has the main features of every hyperbolic equation and is the important part of many mathematical models. In other cases the advection operator is the important part of the full operator of the problem. All 26 schemes have been compared experimentally on the special representative set of tests. Nevertheless to simplicity of the equation, almost all schemes have different disadvantages. They are discussed in detail in the cited papers. So, the investigation of new schemes for this equation is still an important task.
In [4,5,6,7,8,9] some new schemes were constructed for solving this advection equation. The nonlinear monotone Babenko scheme ("square") proved to be the best among all 26 schemes. So, it is a big interest to generalize this scheme to more difficult equations. The important example is a quasi-linear advection equation.
In this paper our basic aim is to construct a quasi-monotone nonlinear Babenko scheme for solving the quasi-linear advection equation and to test it experimentally. The monotonisation of the scheme is done by adding the artificial diffusion with limiters. We also present advanced results of comparative analysis of the new scheme with other known schemes. We have considered explicit and implicit upwind approximation schemes [4,6,13,16] which is firstorder accurate in time and space, the Lax-Wendroff scheme [4] which is the first order accurate in time and second order accurate in space. We also analyze the monotonised "Cabaret" scheme proposed in [10,11]. It is second order accurate in time and space, and its monotonisation is based on apriori knowledge of the dependence region of the exact solution. The authors of this scheme called it by "jumping advection". The considered schemes are compared numerically by using a set of tests, which is similar to one used in [4,5,6,8].

Problem statement
Let us consider the Cauchy problem for the quasi-linear advection equation with finite (mostly) initial condition u(x, 0) = u 0 (x). We are going to construct the Babenko scheme for this equation and to test it on the special set of test solutions. The initial functions for this set have following forms: They have a form of triangle, rectangle, left-hand triangle, right-hand triangle, stepdown and step-up functions, respectively.
Here are the exact solutions for initial conditions defined above: 1.
To compare numerical and exact solutions we use finite-dimensional analogs of norms in spaces C, L 1 , L 2 .

Nonlinear Monotonisation of the Babenko Scheme
For the numerical solution let us introduce the uniform (for simplicity) space-time gridω hτ =ω h ×ω τ : where h is the space grid step and τ is the time grid step.
Further we use the following standard notation of grid functions and operators [15,16] where the superscript is a time layer number and the subscript is a number of point on x. We assume that the solution y n is known. For equation (1.1) the Babenko scheme has the following form: It can be rewritten as Suppose a = 1 2 (ŷ + y), then equation (2.1) takes the form: where γ = aτ /h is the Courant number. In the rest of the paper only difference derivatives are used. Therefore we use the notations for difference derivatives similar to the partial ones.
It is well known [5], that for smooth solutions the Babenko scheme [2, 3] is second order accurate in x and t, but it is non-monotone. This is an implicit scheme, but it does not prevent us from calculating the solution on a next time layer by an explicit formula. For the linear advection equation the Babenko scheme gives the exact solution on some grids, when the corresponding Courant number is equal to 1. Let us notice that first two terms in (2.2) give the usual upwind scheme.
We monotonize the scheme by adding to (2.2) two terms where µ is the artificial diffusion. Since time derivatives u t , u t,−1 of the exact so- , the introduced items produce diffusive terms. Hence the modified scheme takes the following form: Let us note that for µ ≡ 0 we get the original Babenko scheme and for µ ≡ 1 we get the upwind scheme. We now rewrite (2.3) in the following form: .
The scheme (2.3) is constructed taking into account the assumption that y ≥ 0 in the whole space-time region. Then it satisfies the maximum principle [15] if the following conditions are valid Notice that quantityã = 1 2 (y + y −1 ) is used in the literal condition of maximum principle fulfillment. So, it is different from the second inequality in (2.4). But the desirable condition follows from the second inequality in (2.4) ifŷ ≥ 0. Let us introduce functions . Inequalities (2.4) are valid if function µ = µ(R, γ) takes the form [5]: Function µ provides the monotonisation of the considered scheme. Note that the third line in (2.6) corresponds to the line segment, where the scheme preserves the second order of approximation for linear and closed to linear profiles [18,12]. The fourth line in (2.6) provides change from zero to negative values of µ. In this case anti-diffusion items in the scheme appear. They decrease the numerical diffusion at discontinuity points of the solution.
In this work parameter R * is selected experimentally and R * = 1.2 is used in all numerical experiments [5].
Equation (2.3) is nonlinear due to dependence of R onŷ. The unknownŷ value can be calculated using the known value of R by: We can find R using relation µ = µ(R,γ). First we rewrite equation (2.3) in the form: According to (2.5) this equation takes the following form: . (2.9) We denote by b the right-hand side of equation (2.9): .
(2.10) Equation (2.9) with a known right-hand side is a nonlinear equation with respect to R. The left-hand side is a monotone function of R. Considering the form of function µ = µ(R, γ), we get the following relations: (2.11) Remark 1. The value of b is indefinite if y 2 +1 − y 2 = 0, but in this case R = 0 (as it follows from (2.5)). Also, we get from (2.11) that the value of R is indefinite, if b = 0, but in this case the equalityŷ = y can be obtained directly from (2.10) and (2.8).
So, using equations (2.7), (2.10), (2.11) we can find y . Note that γ depends on y , therefore the value of y is calculated numerically using iteration method. In this work at first iteration step we count the value of y as 0.5( y −1 + y +1 ). For this value we calculate γ, then we define R basing on the known value of b and calculate desired value of y using (2.7).

Numerical Experiments
As it was done in [4,5,8] the numerical solution of problem (1.1) is obtained for all stated above initial data for following parameters: The obtained results make possible to analyze qualitative and quantitative features of the new developed scheme and to compare it with the other schemes considered in this work.

Explicit upwind scheme
The solution obtained by this method is very rough, diffusion of solution discontinuity takes place at few steps (from 4 to 10 steps). If γ is smaller, the solution diffusion is greater. The solutions of depression wave type are diffused on more steps compared with shock waves. This solution drawback grows in time. For initial profiles 1, 2, 4 on time intervals where the exact solution evolutes without amplitude change the amplitude of numerical solution diminishes permanently. It is difficult to trace time moments when shock waves arise. The speed of propagation of shock waves for the numerical solution is the same as for exact solution and is independent of γ.

Implicit upwind scheme
The solution obtained using this scheme is even less similar to exact solution than the solution obtained by using the previous one. For this scheme even greater solution diffusion is typical (from 7 to 15 space steps). This diffusion becomes stronger as the Courant number value increases. As shown in Figures 1-6, at initial time moments the bending of back front of the wave is typical for the numerical solution. In process of time the dependence of numerical solution on Courant number becomes less expressed. For initial profiles 1, 2, 3 and 4 the numerical solution obtained by using explicit and implicit upwind schemes evolutes to a wave of triangle profile with diminishing amplitude and shock wave spreading to the right side. It conforms to the exact solution, but the maximum of numerical solution amplitude goes behind the maximum of the exact solution. So, the upwind schemes badly reproduce such features of solution as appearance and propagation of shock and depression waves.

The Lax-Wendroff scheme
The phenomenon of normal dispersion is typical for this scheme. The appearance of shock waves is accompanied by oscillations behind their front. Appearance of oscillations on the left front of generated depression wave is typical for initial profiles 2, 4 and 6. These oscillations are progressing in time. It leads to growth of errors especially in the C norm. For instance at γ = 0.5 for initial profile 2 at time moment t = 10 the relative error in the C norm is equal to 0.17 and at time moment t = T the error is equal to 1.076. The spreading of solution is less for the Lax-Wendroff scheme compared with implicit and explicit upwind schemes. Solution discontinuities are reproduced better by this scheme, and discontinuity diffusion takes place on 3 to 5 steps.

Jumping advection method
Compared with schemes given above the numerical solution obtained using the jumping advection method represents the exact solution better. The main advantage of this method is the ability to represent discontinuities such as shock waves precisely. But due to features of the method the numerical solution is changing spasmodically and the exact position of shock waves can be obtained at some time moments only. For instance if γ = 0.5 and initial profile 5 is used, then the solution jumps to the right side by one step for every fourth time layer, overtaking the exact solution. Thus at these moments there is no computational error.
However, when the Courant number increases, defects of numerical solution such as "steps" and wrong speed of shock waves appear. It can be seen in Figure 9 that "steps" appear for the compression wave for initial profile 4 with any Courant number, and figures 8 and 11 show that when γ = 0.9, then shock wave of the numerical solution move slower than the exact one. In Figures 7-11

Monotonized Babenko scheme
The numerical solution which is obtained using this scheme represents the exact one in the best way compared with schemes examined above in all explored range of Courant numbers. This conclusion follows from the error analysis. This scheme t = 5 has no normal and anomalous dispersion, the solution spreads slightly, discontinuity spreading takes place at smaller number of steps (1-3 steps). For initial profiles 2, 4, 6 little right shifting of a left front of the resulting depression wave can be noticed.
As it is shown on Figures 12-17, numerical solutions are almost the same for different Courant numbers. And for all tests γ = 0.1 corresponds to the most accurate numerical solution.

Conclusion
Nonlinear monotonised Babenko scheme, considered in this work, gives the best solution of quasi-linear advection equation in all explored range of Courant numbers. The obvious advantage of this scheme is its ability to present solution discontinuities in a most accurate way compared with other schemes. We notice that the jumping advection method [11] gives worse solution at some initial conditions and for some Courant numbers, despite of its ability to represent the exact solution at some time moments. It is connected to the appearance of "steps" leading to strong solution distortion. Some additional details can be found in [1].