Gr\"obner Bases and Generation of Difference Schemes for Partial Differential Equations

In this paper we present an algorithmic approach to the generation of fully conservative difference schemes for linear partial differential equations. The approach is based on enlargement of the equations in their integral conservation law form by extra integral relations between unknown functions and their derivatives, and on discretization of the obtained system. The structure of the discrete system depends on numerical approximation methods for the integrals occurring in the enlarged system. As a result of the discretization, a system of linear polynomial difference equations is derived for the unknown functions and their partial derivatives. A difference scheme is constructed by elimination of all the partial derivatives. The elimination can be achieved by selecting a proper elimination ranking and by computing a Gr\"obner basis of the linear difference ideal generated by the polynomials in the discrete system. For these purposes we use the difference form of Janet-like Gr\"obner bases and their implementation in Maple. As illustration of the described methods and algorithms, we construct a number of difference schemes for Burgers and Falkowich-Karman equations and discuss their numerical properties.

Mathematical operations used in the construction of difference schemes for PDEs are substantially symbolic. Thereby, it is a challenge for computer algebra to provide an algorithmic tool for automatization of the difference schemes constructing as well as for the investigation of properties of the difference schemes. One of the most fundamental requirements for a difference scheme is its stability which can be analyzed with the use of computer algebra methods and software [9].
Furthermore, if PDEs admit a conservation law form or/and have some symmetries, it is worthwhile to preserve these features at the level of difference schemes too. In particular, a tool for automatic construction of difference schemes should produce conservative schemes whenever the original PDEs can be written in the integral conservation law form. One of such tools GRIDOP written in Reduce [10,11] is based on symbolic operator methods and generates conservative finite-difference schemes on rectangular domains in an arbitrary number of independent variables. However, the generation is not entirely automatic. A user of GRIDOP has to specify function spaces together with associated scalar products and define grid operators as finite-difference schemes. Then the user may provide partial differential equations in terms of the defined grid operators or the adjoints of those operators. Under these conditions the package returns the finite-difference equations for the dependent variables.
Besides, a few other applications of computer algebra are known to construct finite-difference schemes [12,13] which, being also not completely automatic, are applicable to PDEs of a certain form.
In this paper we describe a universal algorithmic approach to the automatic generation of conservative difference schemes for linear PDEs with two independent variables admitting the conservation law form. This approach generalizes and extends the observations of paper [14] where it was noticed that a conservative difference scheme can be derived as a compatibility condition for a system of difference equations. The system is composed of a discrete form of the original PDEs taken in the integral conservation law form and of a number of natural integral relations between functions and their partial derivatives. The finite-difference scheme is obtained by elimination of all the partial derivatives from the system. We also show, by the example of Burgers equation, that one can also apply the difference elimination approach to generate of difference schemes without use of conservation law form.
To perform the difference elimination we apply the Gröbner bases method invented 40 years ago by Buchberger [15] for polynomial ideals. This method has become the most universal algorithmic tool in commutative algebra and algebraic geometry and found also numerous fruitful applications for computations in certain noncommutative polynomial rings as well as in rings of linear differential operators and differential polynomials [16]. Nowadays, all modern generalpurpose computer algebra systems, for example, Maple [17] and Mathematica [18], have special built-in modules implementing algorithms for computing Gröbner bases. However, the fastest implementation of these algorithms for commutative polynomial algebra is done in the specialpurpose systems Singular [19] and Magma [20]. As to the difference algebra [21], in spite of known for long time (see [22] and references therein) extensions of Buchberger's algorithm [23] to difference polynomial rings, there are only a few implementations of the algorithm specialized to shift Ore algebra: in the Ore algebra library package of Maple [24], in the library OreModules [25] developed using the latter package and in Singular (Plural) [26]. These packages can be used for computing Gröbner bases of linear difference ideals and modules, and, in particular, for those linear systems which are considered below.
In the given paper we present, however, another algorithm for computing difference Gröbner bases. This algorithm is superior over our Janet division algorithm whose polynomial version [27] in most cases is computationally more efficient than Buchberger's algorithm [28]. In addition, unlike the above mentioned implementations of Buchberger's algorithm for the shift Ore algebra, the algorithm described below and its recent implementation in Maple [29] admit a natural extension to nonlinear difference systems exactly in the same way as differential involutive algorithms [30,31,32]. The algorithm improves our Janet-like division algorithm [33] adapted to linear difference ideals [34]. The improvement includes, in particular, the difference form of the involutive criteria [27] modified for Janet-like reductions. These criteria allow to avoid some useless reductions, and thereby accelerate the computation.
The structure of the paper is as follows. In Section 2 we describe the basic idea of our approach to the generation of finite-difference schemes for two-dimensional PDEs. Section 3 contains definitions and notions of difference algebra which are used in the sequel. Here we define Gröbner bases for linear difference ideals and their special form called Janet-like bases.
In Section 4 we present an improved version of the algorithm in paper [34] and briefly discuss some relevant computational aspects. Section 5 illustrates our approach to construction of difference schemes by simplest second-order equations -Laplace's equation, the wave equation, the heat equation, and by the first-order advection equation. In Section 6 we generate several difference schemes for Burgers equation. But for all that, to avoid problems arising in computing of nonlinear Gröbner bases, we denote the square of the dependent variable by an extra function. Besides, we characterize some of the constructed schemes by the modified equation method. In Section 7 we consider the two-dimensional quadratically nonlinear Falkowich-Karman equation describing transonic flow in gas dynamics. Here, we succeeded in computing of the nonlinear Gröbner basis by hand, and in that way generated the cubic nonlinear difference scheme which possesses some attractive properties. These properties as well as those of the schemes generated for Burgers equation are illustrated by some numerical experiments in Section 8. We conclude in Section 9.
It follows that eliminating from (6) all the proper grid partial derivatives gives equations containing only independent (vector) function u, and, hence composing a finite-difference scheme. If system (6) is linear, then this difference elimination can always be algorithmically achieved by the Gröbner bases method considered in the next section.
It should be noted that generation of finite-difference schemes on grid (3) by the elimination can be also applied to PDEs irrelative to their conservation law properties. Again, one has to add to the initial differential equations, written in terms of grid variables (4), the corresponding number of integral relations (5) and approximate them by numerical quadrature formulas. Such an approach may give more flexibility in generation of distinct difference schemes, and we apply it in Section 4 to the first-order advection equation, and in Section 5 to Burgers equation. Generally, however, for the second-(and higher-) order PDEs admitting the integral conservation law form, the difference scheme obtained directly from the differential form may not be conservative. Besides, the difference elimination based on the integral form is usually more efficient than that based on the differential form. This is because the number of partial derivatives to be eliminated in the former case is smaller than in the latter case. Indeed, the integrand in (2) has the differential order smaller by one than that in (1) whereas computational complexity of the elimination is at least exponential in the number of eliminated variables [35].

Difference Gröbner bases
A difference ring R is a commutative ring with a unity together with a finite set of mutually commuting injective endomorphisms θ 1 , . . . , θ n of R. Similarly, one defines a difference field. Elements {y 1 , . . . , y m } in a difference ring containing R are said to be difference indeterminates over R if the set θ k 1 1 · · · θ k 1 1 y j | {k 1 , . . . , k n } ∈ Z n ≥0 , 1 ≤ j ≤ m is algebraically independent over R.
The field Q(x 1 , . . . , x n ) of rational functions in {x 1 , . . . , x n } whose coefficients are rational numbers is an example of difference field, and we shall assume in the next sections that the coefficients of PDEs belong to this field.
Let K be a difference field, and R := K{y 1 , . . . , y m } be the difference ring of polynomials over K in variables {θ µ • y k | µ ∈ Z n ≥0 , k = 1, . . . , m}. Hereafter, we denote by R L the set of linear polynomials in R and use the notations: A difference ideal is an ideal I ⊆ R closed under the action of any operator from Θ. If F := {f 1 , . . . , f k } ⊂ R is a finite set, then the smallest difference ideal containing F will be denoted by Id(F ). If for an ideal I there is F ⊂ R L such that I = Id(F ), then I is a linear difference ideal.
A total ordering ≻ on the set of θ µ • y j is a ranking if ∀ i, j, k, µ, ν the following hold: Given a ranking ≻, a linear polynomial f ∈ R L \ {0} has the leading term of the form a ϑ • y j , ϑ ∈ Θ, where ϑ • y j is maximal w.r.t. ≻ among all θ µ • y k which appear with nonzero coefficient in f . lc(f ) := a ∈ K \ {0} is the leading coefficient and lm(f ) := ϑ • y j is the leading (head) monomial.
A ranking acts in R L as a monomial order. If F ⊆ R L \ {0}, then lm(F ) will denote the set of the leading monomials and lm j (F ) will denote its subset for indeterminate y j . Thus, Given a nonzero linear difference ideal I = Id(G) and a ranking ≻, the ideal generating set G = {g 1 , . . . , g s } ⊂ R L is a Gröbner basis [22,23] It follows that the head monomial of f ∈ I \ {0}, as well as the polynomial f itself, is reducible modulo G and yields the head reduction: If f ′ = 0, then its leading monomial is again reducible modulo G, and, by repeating the reduction finitely many times [16,22,23] we obtain f − → G 0. Generally, if a polynomial h ∈ R L contains a term with monomial u and coefficient c = 0 such that u = ϑ • lm(f ) for some ϑ ∈ Θ and f ∈ F ⊂ R L \ {0}, then h can be reduced: By applying the reduction finitely many times, one obtains a polynomialh which is either zero or such that all its (nonzero) terms have monomials irreducible modulo the set F ⊂ R L . In both casesh is said to be in the normal form modulo F (denotation: In our algorithmic construction of reduced Gröbner bases we shall use a restricted set of reductions called Janet-like (cf. [33,34]) and defined as follows.
For a finite set F ⊆ R L and a ranking ≻, we partition every lm k (F ) into groups labeled Now we characterize a monomial u ∈ lm k (F ) by the nonnegative integer λ i : Let DP (f, F ) denotes the set of all difference powers for f ∈ F . Now we define the partition of the set Θ into two disjoint subsets which is similar to the partition of monomials into nonmultiplicative and multiplicative ones in the involutive approach [27]. A finite basis G of I = Id(G) is called Janet-like [33,34] if In full analogy with (9) a J -reduction is defined as Apparently, any element in the ideal I = Id(G) is J -head reduced to zero by the finite sequence of J -head reductions by elements g ∈ G in the Janet-like basis G: If the leading monomial of p ∈ R \ {0} is not J -reducible modulo a finite subset F ⊂ R \ {0} we say that p is in the J -head normal form modulo F and write p = HN F J (p, F ). If none of monomials in p is J -reducible modulo F we say that p is in the (full) normal form modulo F and write p = N F J (p, F ).
Since J -reducibility implies the Gröbner reducibility (9), a Janet-like basis satisfying (10) is a Gröbner basis. The converse is generally not true, that is, not every Gröbner basis is Janet-like.
The properties of a Janet-like basis are very similar to those of a Janet basis [36], but the former is generally more compact than the latter. For all that we consider hereafter only minimal bases.
Let GB be a reduced Gröbner basis, satisfying [23]: Let now JB be a Janet basis, and JLB be a Janet-like basis of the same ideal and for the same ranking. Then for their cardinalities the inclusion Card(GB) ≤ Card(JLB) ≤ Card(JB) holds [27,33]. Here Card abbreviates cardinality, that is, the number of elements. Whereas the algorithmic characterization of a Gröbner basis is zero redundancy of all its S-polynomials [15,23], the algorithmic characterization of a Janet-like basis G has the following form (cf. [33]): These conditions are at the root of the algorithmic construction of Janet-like bases as described in the following section.

Algorithm
In this section we present an algorithm for constructing a reduced Gröbner basis (8) of the ideal generated by an input set of linear difference polynomials. The algorithm is an improved version of the algorithm in paper [34] and translates to the difference form of the polynomial involutive algorithm [27] modified for the Janet-like reductions. To apply the difference form of criteria to avoid some unnecessary reductions we need the following definition.
An ancestor of a difference polynomial f ∈ F ⊂ R L \{0} is a polynomial g ∈ F of the smallest deg(lm(g)) among those satisfying f = θ • g modulo Id(F \ {f }) with θ ∈ Θ. If for all that deg(lm(g)) < deg(lm(f )), then the ancestor g of f is called proper.
If an intermediate polynomial h that arose in the course of the below algorithm has a proper ancestor g in the intermediate basis G, then h has been obtained from g via a sequence of shift operations of the form ϑ • g where ϑ ∈ DP (g, G) with lm(ϑ • g) J -irreducible modulo G. For the ancestor g itself the equality lm(anc(g)) = lm(g) holds.
In the main algorithm GröbnerBasis and its subalgorithms we endow every element f ∈ G in the intermediate set of difference polynomials G (occuring in the set T ) with a triple structure of the form: where pol(p) := f is the polynomial f itself, anc(p) := g is an ancestor of f in G, dp(p) := dpow is a (possibly empty) subset of DP (f, G).
The set dpow associated with the polynomial f accumulates all the difference powers for f which have been already applied to f in the course of the algorithm. Keeping this information serves to avoid useless repeated applications of the difference power operators. Knowledge of ancestors for elements in the intermediate basis helps to avoid some unnecessary reductions by applying Buchberger's chain criterion [23] adapted to Janet-like reductions.
Algorithm: GröbnerBasis(F, ≺) Input: F ∈ R L \ {0}, a finite set; ≺, a ranking Output: G, a reduced Gröbner basis of Id(F ) choose p ∈ Q with the lowest lm(pol(p)) w.r.t. ≻ 7: if pol(p) = anc(p) then 9: T := T ∪ {h, anc(p), dp(p)} 15: for all q ∈ T and ϑ ∈ DP (q, T ) \ dp(q) do 16: In the above main algorithm GröbnerBasis and its subalgorithms presented below, where no confusion can arise, we simply refer to the triple set T as the second argument in DP , N F J , and HN F J instead of the polynomial set {g = pol(t) | t ∈ T }. Sometimes we also refer to the triple p instead of pol(p). Besides, when we speak of reduction of the triple set Q modulo triple set T we mean reduction of the polynomial set Correctness and termination of algorithm GröbnerBasis can be shown exactly as in the polynomial case [27,33]. Here we only elucidate some related features of the algorithm.
At steps 4 and 19 the J -head reduction is performed for the difference polynomials in Q modulo those in T . Then the remaining tail reduction is done in line 13 to obtain the (full) J -normal form. Thereby, the main while-loop 5-20 terminates when the conditions (14) hold for the difference polynomial set G composed from the first elements of triples in T and the set Q is empty. The upper for-loop 9-11 provides minimality of the output Janet-like basis contained in T [27]. Another for-loop 15-18 constructs new conditions (14) which have to be further examined because of the insertion of a new element in T at step 14. Besides, the set of difference powers is upgraded at step 17.
Furthermore, the main algorithm GröbnerBasis together with its subalgorithms presented below ensures that every element in the output Janet-like basis composed from the first elements in the triple set T has one and only one ancestor. This ancestor is apparently irreducible, in the Gröbner sense (9), by other elements in the basis. Thereby, those elements in the output basis that have no proper ancestors constitute the reduced Gröbner basis (13) that is returned by the main algorithm at the last step 21.
The algorithm HeadReduce invoked in lines 4 and 19 of the main algorithm returns the set Q which, if nonempty, contains part of the intermediate basis J -head reduced modulo T . The reductions are performed by its subalgorithm HeadNormalForm that is invoked at step 6 of the algorithm.
If algorithm HeadNormalForm returns h = 0, then lm(h) does not belong to the initial ideal generated by {lm(pol(f )) | f ∈ Q ∪ T } [27,36]. In this case the triple {h, h, ∅} for h is inserted (line 9) into the output set Q. Otherwise, the output set Q retains the triple p as it is in the input.
In the case when h = 0 and pol(p) has no proper ancestors that is checked at step 14, all the descendant triples for p, if any, are deleted from the intermediate set S at step 16. Such descendants cannot occur in T owing to the choice conditions at steps 1, 6 and to the displacement condition of step 9 in the main algorithm GröbnerBasis. Steps 14-18 serve for the memory saving and can be ignored if the memory restrictions are not very critical for a given problem. In this case all those descendants will be casted away by the criteria checked in the below algorithm HeadNormalForm.
Algorithm HeadNormalForm performs verification (step 3) of J -head reducibility of the input polynomial h modulo the polynomial set (15). This verification consists in searching a difference polynomial (reductor) in the set G defined in (15) such that G yields the reduction (11). If the search fails, that is, there is no J -reductor, the algorithm returns at step 4 the input polynomial.
In addition, if all difference polynomials in the input set F for the main algorithm Gröbner-Basis have constant coefficients, then the set of criteria (16) can be enlarged with one more criterion C 4 : C 4 (p, g) is true for lm(pol(p)) = θ • y k , lm(pol(g)) = ϑ • y k ⇐⇒ lcm(θ, ϑ) = θ ϑ .
The last algorithm TailNormalForm completes J -reduction of the J -head reduced polynomial in the input triple by performing its J -tail reduction. This algorithm is invoked at step 13 of the main algorithm GröbnerBasis. The tail reduction is performed in the while-loop as a sequence of elementary reductions (11).
Algorithm: TailNormalForm(p, T, ≺) Input: p, a triple such that pol(p) = HN F J (p, T ); T , a set of triples; ≺, a ranking Output: h = N F J (p, T ), the (full) J -normal form of pol(p) modulo T 1: Because of the lack of an appropriate collection of benchmarks for linear finite-difference polynomial systems, the algorithmic efficiency of algorithm GröbnerBasis can be indirectly analyzed by running its polynomial (non-difference) counterpart [27,33] for the extensive benchmarks collection in [38,39]. Some timings for our polynomial implementation can be found on the Web page [28].
Recently, the algorithm in its difference version was implemented in Maple [29]. Just this implementation was used for generation of linear finite-difference schemes as described in the next sections. Though one needs special and intensive benchmarking for linear difference systems, our first experimenting with the Maple implementation and with that for commutative polynomials gives us a good reason to expect that the following merits revealed for the pure polynomial version [27] hold also for the difference one: • automatic avoidance of some useless reductions; • weakened role of the criteria: even without applying any criteria the algorithm is reasonably fast; • smoothed growth of intermediate coefficients; • fast search of a reductor which provides the elementary Janet-like reduction (11) of a given term. It should be noted that there can be at most one reductor [27]; • natural and effective parallelism.

Laplace equation
In this section we illustrate the approach of Section 2 to the automatic generation of difference schemes by simplest elliptic, parabolic and hyperbolic equations. To compute Gröbner bases providing the elimination of the partial derivatives to construct difference schemes we used the Maple package [29] implementing the algorithms described in the previous section. We start with the Laplace equation [3,4,5,6,7] u xx + u yy = 0 (17) and rewrite it as the conservation law (1) Now we add the relations (5) for the partial derivatives u x and u y x j+2 x j Thus, we obtain the system of three integral relations (18), (19) for three functions u(x, y), u x (x, y), u y (x, y).
To discretize this system we choose the rectangular contour of Fig. 1 on the orthogonal and uniform grid (3) with and use the midpoint integration method for both (18) and (19). This yields the system: Rewritten in terms of difference polynomials in the ring Q{u, u x , u y } (see Section 2) it reads: Computation of the Gröbner basis for the elimination ranking (Section 2) with u x ≻ u y ≻ u and θ x ≻ θ y gives: The latter equation with eliminated u x and u y is the standard difference scheme with the central approximation of the second-order derivatives in (17) written in double nodes Similarly, the trapezoidal integration rule for relations (19) generates the same difference scheme but written in ordinary nodes

Heat equation
Consider now the heat equation [3,4,5,6,7,8] u t + αu xx = 0 in its conservation law form The integrand in (21) contains the only partial derivative u x . Hence we add the single integral relation Again we discretize u(x, t) and u x (x, t) on the orthogonal and uniform grid with the spatial mesh step h and the temporal mesh step τ , and choose the contour shown in Fig. 2. Then, applying the midpoint rule for the contour integral and the trapezoidal rule for the relation integral we find two difference equations for two indeterminates u, u x in the form By elimination of u x by means of the Gröbner basis with u x ≻ u we obtain the famous Crank-Nicholson scheme [3,4,5,6,7,8,9] The same scheme is obtained for the midpoint integration method applied to (22).

Wave equation
The wave equation [3,4,5,6,7,8] u tt − u xx = 0 in the conservation law form is given by Choosing the same grid with the mesh steps (20), the contour of Fig. 1 and integral relations (19) as are used in Section 5.2 for the Laplace equation (18) and applying the midpoint rule for the contour integral and the trapezoidal rule for the integral relations we obtain the operator equations The Gröbner basis method yields the standard difference scheme
Being of first order, the equation (23) has already the conservation law form (1). By this reason, to generate a difference scheme we shall not convert the equation into the integral form (2). In the latter case one has nothing to eliminate. Instead, we consider equation (23) together with the integral relations: Discretization of u, u t and u x on the orthogonal and uniform grid with the mesh steps h and τ and the explicit integration formula for the upper integral relation in (24) together with the midpoint integration rule for the lower relation give the difference system: Let us apply the operator θ x to both sides of the middle equation in (25) and then use the Lax method, that is, replace θ x with (θ 2 x + 1)/2 in the second term of the right-hand side. This replacement yields • u = 0. Its last element gives the scheme: 6 Burgers equation 6

.1 Conservation law form
Consider Burgers equation [4,5,8,9] in the form where we replaced u 2 by the flux function f in order to avoid computation of nonlinear difference Gröbner bases. ν is called the viscosity. This equation exhibits some difficult features from the point of view of simple finite difference schemes due to the term f = u 2 . Let us, first, convert equation (26) into the conservation law form Then choose the contour of Fig. 1 and add the integral relation Denoting as above the temporal and spatial mesh steps by τ and h, and applying the midpoint integration rule we obtain the system: Its Gröbner basis form for the elimination order with u x ≻ u ≻ f and θ t ≻ θ x is given by The obtained difference scheme is the standard explicit scheme with forward time and forward space differencing. It is wellknown that schemes of this type are unstable [5,8,9]. Furthermore, by using implicit schemes one can provide the von Neumann stability 1 . However, all such schemes are usually not very satisfactory when one considers non-smooth or discontinuous solutions (shock waves) of Burgers equation.

Lax method
To exploit more flexibility and freedom in our difference elimination approach to generation of finite-difference schemes, we go back to the original differential equation (26) and consider it together with the integral relations providing the elimination. For discretization of the relations we combine the midpoint rule for integration over x with the explicit integration over t and apply the Lax method to the last integration: The Gröbner basis based elimination with One can also use the trapezoidal rule for the spatial integrations. This derives other schemes. Since there are three spatial integrals in (28), by selecting either the midpoint or the trapezoidal rule for these integrals, we obtain eight possible variants of the difference schemes. Our computation with the Gröbner bases reveals seven different schemes. Apart from (29) there are 2(u n+1 j+3 + 2u n+1 j+2 + u n+1 j+1 ) − (u n j+4 + 2u n j+3 + 2u n j+2 + 2u n j+1 + u n j ) 8τ Just the scheme (34) is obtained twice in the course of generating eight schemes.

Two-step Lax-Wendroff method
Our Gröbner basis based technique can also be applied to generate other types of difference schemes. For example, one can generate two-step Lax-Wendroff schemes [41]. Let u and f denote the values of u and f on the intermediate time levels. Then, applying again the midpoint rule for the spatial integrals, gives the following difference system: For the elimination ranking with the Gröbner basis contains the Lax-Wendroff scheme With all possible combinations of the trapezoidal and midpoint rules one obtains 49 different Lax-Wendroff schemes.

Differential approximation
To analyze properties of a difference scheme it can be useful to compute its differential approximation [40] that is often called the modified equation(s) of the difference scheme. There are whole classes of different schemes for which their stability properties can be obtained with the aid of the differential approximation [2]. For all that, in many cases, the computation can be easily done with modern computer algebra software. In our research we use Maple [17]. Consider, for example, the schemes (29)-(35).
Their differential approximation for f = u 2 and with collection of the coefficients at τ , h 2 , h 2 /τ is given by: The schemes (29)-(35) differ in the coefficient ( * ) at h 2 only. These coefficients are as follows Thereby, comparison of differential approximations for schemes (29)- (35) shows that • all the schemes provide the same order of approximation in τ , h; • they have identical linear numerical dissipation (viscosity) [4,9] determined by u xx h 2 /(2 τ ); • the schemes possess similar dispersion properties with distinction in the rational coefficients of the differential polynomial in u multiplied by h 2 .
As to scheme (27), the right-hand side of its differential approximation reads This explicitly shows instability of the scheme which does not yield linear numerical viscosity. We obtained also analogous results on stability and on close properties for the different Lax-Wendroff schemes of type (36) and its variations due to the choice of different numerical integration rules for the spatial integrals.

Godunov method
It is especially difficult to simulate numerically nonsmooth and discontinuous solutions which are among most interesting problems in computational fluid and gas dynamics [1,4,5,8]. Most of the known difference schemes fail to handle these singularities. The most appropriate numerical approach to such problems was developed by Godunov [1,42] and based on solving a local Riemann problem [4,6] as a cornerstone of the Godunov's scheme generation. There are special numerical Riemann solvers, for example [43], designed for these purposes and for application to computational fluid dynamics.
Instead of the use of numerical Riemann solvers, we apply the Gröbner bases technique to generate the Godunov-type scheme for inviscid Burgers equation when ν = 0 in (26). For this purpose we discretize the corresponding system in (28) in the following way Here, the third equation contains in its left-hand side the product of two different solutions for the flux function f of the local Riemann problem [43]. Therefore, we add to the system composed of the original differential equation and discrete forms of the integral relations for partial derivatives u t , u x , u xx the nonlinear difference equation on f and f x containing solutions of the local Riemann problem.
Since the Riemann condition on the flux is now a constituent of the difference system, an elimination of all the partial derivatives of u and f gives the difference scheme consistent with that condition. To do the elimination from the nonlinear system (37) we apply the Gröbner factoring approach [44]: if a Gröbner basis contains a polynomial which factors, then the computation is split into the computation of two or more Gröbner bases corresponding to the factors. In doing so, we choose the elimination ranking and compute two Gröbner bases, for every factor in (37). Then we compose the product of two obtained difference polynomials in u and f that gives us the Godunov-type difference scheme: Below (Section 7) we compare schemes (29), (36) and (38) by numerical simulation of a discontinuous solution.

Falkowich-Karman equation
Consider now the nonlinear two-dimensional Falkowich-Karman equation [45] ϕ xx (K − (γ + 1)ϕ x ) + ϕ yy = 0 (39) describing transonic flow in gas dynamics in its non-stationary form This form can be used to find a stationary solution by the steady-state method. We rewrite equation (40) into the conservation law form Here we use a grid decomposition of the three-dimensional domain in (x, y, t) into elementary volumes. Fig. 3 shows an elementary volume. Again we add the integral relations for partial derivatives with the use of the trapezoidal integration rule for ϕ x , ϕ y and the midpoint rule for ϕ t .
Then we obtain the nonlinear operator equations: Because of nonlinearity in the initial differential equation (40), the difference system obtained is also nonlinear. By this reason the Maple package [29] implementing algorithm GröbnerBasis (Section 4) is inapplicable to (41). Since there is no software available for computing difference Gröbner bases for nonlinear systems, we had to perform hand computations in accordance with j + 2 j + 1 j k k + 1 k + 2 P P P P P P P P P! ! ! ! ! ! ! P P P P P P P P P P ! ! ! ! ! ! P P P P P P P P P P P P P P P P P P! ! ! ! ! ! ! P P P P P P P P P P ! ! ! ! ! ! P P P P P P P P P u u u u u u u ! ! ! ! ! ! n n + 1 Figure 3. Cell for Falkowich-Karman equation.
the above described algorithm and with an assistance of Maple to check some intermediate results.
In these calculations we used the lexicographical ranking such that ϕ x ≻ ϕ y ≻ ϕ t ≻ ϕ and θ x ≻ θ y ≻ θ t . The resulting Gröbner basis has the form: The last element is the finite-difference scheme for equation (40): By construction, this scheme is fully conservative and does not involve switches that is typical for computing transonic flows [47].  In its stationary form scheme (42) is related to equation (39). Here symbols D x and D y are the forward differencing operators and D xx and D yy are the central second-order differencing operators with respect to x and y.
The stencil for scheme (43) is shown in Fig. 4. It should be noted that, unlike the original differential equation (40) which is quadratically nonlinear, both schemes (42) and (43) have the the cubic nonlinearity in the grid function. This is the rigorous algebraic consequence of the difference system (41). In accordance to the wellknown fact [46], that algebraic elimination of variables from a nonlinear system leads generally to increase of its degree of nonlinearity.
As an application of this scheme, in the next section we consider an example of one-dimensional transonic flow with shock-wave taken from [47]. Here H(y) is the Heaviside step function [48] whose derivative is the Dirac delta function: Physically, the solution (46) defined by the initial condition (45) represents a shock wave which moves with constant speed (u l + u r )/2 without changing its shape. In our numerical simulation the values of u l and u l were chosen as 0.8 and 0.2. The pictures below demonstrate the numerical solution of the Riemann problem (44), (45) at time t = 2/3. Solid line shows the exact solution (46), and the numerical results are depicted by green dots. For the ratio τ /h of mesh steps which is called Courant (or Courant-Friedrichs-Levy) number [1] we have chosen the two values 0.9 and 0.1.
All schemes are numerically stable. For schemes (29) and (36) their stability is analytically showed by the differential approximation (Section 6.4). Because of the nonlinearity in f x in the third equation of Godunov scheme (38) we did not compute the differential approximation for this scheme.
As can be expected (cf. [3]), the dispersion effects in schemes (29)-(36) become stronger for the smaller value of the Courant number (Figs. 6 and 8). Qualitatively [49,50], the behavior of Lax scheme (29) in Figs. 5, 6 is typical for the classical first-order schemes when they are applied to problem (44)- (45) whereas the Lax-Wendroff scheme (36) behaves as the secondoder method. The Godunov scheme (38), as a shock capturing one (cf. [49,50]), is much better for numerical description of solution (46) than the schemes (29)-(36), and does not reveal its sensitivity to the value of the Courant number.

Falkowich-Karman equation
Now we consider the application of difference scheme (43) to the one-dimensional stationary transonic flow in a channel with a straight density jump [47]. The exact shock-wave solution of equation (39) at 0 ≤ x ≤ 1 is shown in Figs. 11 and 12 by solid red line. Circles depict the numerical data obtained from difference scheme (43). As an initial approximation, the parabola was chosen satisfying the following boundary conditions: at the left, both the function and its derivative are fixed by the values from the exact solution; at the right, the only function is bound to the exact solution.
As one can see from Fig. 12, scheme (43) possesses a stable and uniform convergence to the exact shock-wave solution. Because, by its construction, the scheme is fully conservative, it does not reveal non-uniqueness of solutions that is typical for the traditional difference schemes [47].
Moreover, the size of the shock transition zone is just one spatial mesh step that is a consequence of preserving at the discrete level of all algebraic properties of the initial PDE (40). This is a result of algebraic difference elimination provided by the Gröbner bases method. Another merit of scheme (43) is that it does not involve switches that is typical for computing transonic flow as we already pointed out in Section 7.
This example shows a principal possibility of constructing difference schemes for transonic flow without switches and with the same stencil for both subsonic and supersonic flow.

Conclusion
In the present paper we have shown that the Gröbner bases method, being a universal algorithmic tool for linear difference algebra, can be effectively applied to the construction of differences schemes for linear PDEs with two independent variables and with rational function coefficients. Owing to the Gröbner bases, this construction is an algorithmic procedure. It consists in elimination of partial derivatives from the system of difference equations composed from a discrete version of the original PDEs (on an orthogonal uniform grid) and numerically approximated integral relations between the unknown functions and their partial derivatives. As this takes place, the difference scheme obtained may depend on the choice of the integration contour and numerical approximations for integral relations.
The method is especially efficient when a PDE or a system of PDEs admits the conservation law form. In this case the difference schemes obtained are fully conservative. The structure of a scheme generated may depend on the choice of integration contour and numerical integration rules. In so doing, it is not clear a priori which integration rule leads to a better scheme.
We also described an efficient algorithm for the construction of Gröbner bases for linear difference ideals. The algorithm is based on the concept of Janet-like reductions. Its first implementation in Maple is already available, and we used this implementation in the generation of all linear difference schemes presented in the paper.
For classical linear PDEs such as the Laplace equation, the Heat equation, the Wave equation and the Advection equation our algorithmic technique leads to the well-known finite difference schemes. For Burgers equation we generated several schemes based on the Lax and Lax-Wendroff methods and computed their numerical dissipation and dispersion by the differential approximation (modified equation) method. By example of Burgers equation we also demonstrate that it is possible to combine the Godunov method with Gröbner bases to derive a shock capturing scheme.
The non-traditional cubic nonlinear difference scheme generated by our difference elimination method for the Falkowich-Karman equation describing transonic flow in gas dynamics possesses a number of attractive properties in comparison with traditional schemes. Among them there are a stable convergence in time to the exact solution with a one-dimensional shock wave and absence of switches. It should be noted, however, that due to its cubic nonlinearity, scheme (43) has a slower convergence in comparison with the traditional schemes specially optimized for numerical simulation of transonic flows in gas dynamics. By this reason one needs additional research for optimizing nonlinear schemes obtained by the difference elimination.
As we already mentioned in the introduction, algorithm GröbnerBasis admits a generalization to polynomial-nonlinear systems of difference equations exactly in the same way as the differential involutive algorithm of paper [32]. In doing so, if every equation in the initial system is linear with respect to the highest ranking difference term and this property of the system is not violated during its completion to involution, then algorithm GröbnerBasis will work correctly and provide the desirable output. Such is indeed the case for system (41). In the most general case of a difference system with polynomial nonlinearity, it can be split into a finite number of subsystems such that every subsystem can be converted into the Gröbner basis form by applying our algorithm. The underlying splitting algorithm is a difference analogue of that described in [31]. The latter algorithm is similar to the splitting algorithm implemented in the library package diffalg in Maple.
The above described approach can be also generalized to PDEs with three and more independent variables. Thus, if PDEs admit the conservation law form, then one can use multidimensional analogues of equations (1) and (2) together with their elementary volume discretization.