Strong Consistency and Thomas Decomposition of Finite Difference Approximations to Systems of Partial Differential Equations
SStrong Consistency and Thomas Decomposition ofFinite Difference Approximations to Systems ofPartial Differential Equations
Vladimir P. Gerdt ∗ Daniel Robertz † Yuri A. Blinkov ‡ Abstract
For a wide class of polynomially nonlinear systems of partial differential equationswe suggest an algorithmic approach that combines differential and difference algebra toanalyze s(trong)-consistency of finite difference approximations. Our approach is appli-cable to regular solution grids. For the grids of this type we give a new definition ofs-consistency for finite difference approximations which generalizes our definition givenearlier for Cartesian grids. The algorithmic verification of s-consistency presented in thepaper is based on the use of both differential and difference Thomas decomposition. First,we apply the differential decomposition to the input system, resulting in a partition of itssolution space. Then, to the output subsystem that contains a solution of interest we ap-ply a difference analogue of the differential Thomas decomposition which allows to checkthe s-consistency. For linear and some quasi-linear differential systems one can also applydifference Gr¨obner bases for the s-consistency analysis. We illustrate our methods andalgorithms by a number of examples, which include Navier-Stokes equations for viscousincompressible flow.
In the given paper we consider systems of partial differential equations (PDE): f = · · · = f p = 0 , F := { f , . . . , f p } ⊂ R . (1)Here f i ( i = 1 , . . . , p ) are elements of the differential polynomial ring R := K{ u } , the ring ofpolynomials in the dependent variables u := { u (1) , . . . , u ( m ) } ( differential indeterminates ) andtheir partial derivatives, which are obtained by applying the power products of the pairwisecommuting derivation operators ∆ := { ∂ , . . . , ∂ n } ( ∂ j ≡ ∂ x j ). We shall assume that thecoefficients of the polynomials are rational functions in a := { a , . . . , a l } , a finite number ofparameters (constants), whose coefficients are rational numbers, i.e., K := Q ( a ). One can alsoextend the last field to Q ( a , x ), where x := { x , . . . , x n } is the set of independent variables.Equations (1) arise in mathematical descriptions of many processes in natural sciences,e.g., in continuous mechanics and physics, whose dynamics evolve in space-time. Apart from ∗ Laboratory of Information Technologies, Joint Institute for Nuclear Research, 6 Joliot-Curie Str., 141980Dubna, Russian Federation and Peoples’ Friendship University of Russia (RUDN), 6 Miklukho-Maklaya Str,Moscow, 117198, Russian Federation, [email protected] † School of Engineering, Computing and Mathematics, University of Plymouth, 2-5 Kirkby Place, DrakeCircus, Plymouth PL4 8AA, Devon, United Kingdom, [email protected] ‡ Faculty of Mathematics and Mechanics, National Research Saratov State University, 83 AstrakhanskayaSt, Saratov, 410012, Russian Federation, [email protected] a r X i v : . [ c s . S C ] S e p ery special cases, exact solutions to the governing PDE system are unknown and only numer-ical solutions can provide valuable information in the study of the process under investigation.For their numerical solution the differential equations in (1) have to be replaced by discretecounterparts. The most widely used methods of discretization and numerical solving are themethod of finite elements, the method of finite volumes and the method of finite differences.The last method is historically the first [52] and is based on replacing differential equationsby difference ones defined on a chosen solution grid. In order to construct a numerical so-lution, the devised finite difference approximation (FDA) to PDE is complemented with anappropriate discretization of initial or/and boundary condition(s) for the PDE. As this takesplace, the quality of the numerical solution to PDE crucially depends on the quality of itsFDA (difference scheme).The main requirement for an FDA is the convergence of a numerical solution to a solutionof PDE in a limit when the grid spacings tend to zero. However, except for a very limited classof problems (see [54], Thm. 1.5.1), convergence cannot be directly established. In practice,given an FDA, its consistency and stability are analyzed as the necessary conditions forconvergence. Consistency implies reduction of the FDA to the original PDE when the gridspacings tend to zero, and stability provides boundedness of the error in the solution undersmall perturbation in the numerical data. It is pertinent to note that in the case of nonlinearFDA (scheme) its theoretical stability analysis is highly conjectural and it is usually studied“experimentally” by division of the grid spacings in halves or by comparison with the exactsolution if it is known.One of the most challenging problems is to construct such difference approximations ofequations (1) which preserve their basic algebraic properties in the discrete setting, e.g.,the continuous identities and theorems of vector calculus, symmetries and conservation laws.Such discretizations, which are called compatible or mimetic [3, 4, 10], and sometimes structurepreserving [11], are more likely to produce highly accurate and stable numerical results, aswas observed in numerous computational experiments (cf. [34]). The most universal approachto determine invariant solutions of initial and boundary value problems for systems (1) andto derive their conservation laws is the Lie symmetry analysis [44]. Certain counterpartsof continuous symmetries for single differential equations were studied for finite differenceschemes in [14, Ch. 4].In [23, 18] for systems (1) and Cartesian (i.e., rectilinear and equisized) solution grids weintroduced the novel concept of strong consistency , or s-consistency , of FDA to PDE, whichstrengthens the concept of consistency. Loosely speaking, s-consistency of an FDA meansnot only approximation of (1) by the FDA, but also approximation of every element in theradical differential ideal, generated by { f , . . . , f p } , by an element in the perfect differenceideal generated by the difference polynomials in the FDA. In the subsequent papers [1, 2], bycomputational experiments with two-dimensional incompressible Navier-Stokes equations, itwas shown that FDA which are s-consistent have much better numerical behavior than FDAwhich are not. To verify s-consistency of linear FDA (to linear PDE) in [23] we used the algo-rithms and software for constructing differential and difference Janet/Gr¨obner bases [8, 24].The generalization to nonlinear PDE given in [18], based on the concept of difference Gr¨obnerbasis, is not algorithmic, because the difference polynomial ring is non-Noetherian [22, 36]and the basis may be infinite. In the conference paper [25] we suggested an algorithm forverification of s-consistency on Cartesian grids that is based on investigating the FDA bya difference analogue of differential Thomas decomposition . In the special case of binomialperfect difference ideals another kind of decomposition was suggested in [16].The notion of differential Thomas decomposition and its algorithmic construction stemmedfrom the Riquier-Janet theory [49, 30]. Wu Wen-tsun was the first who showed [61] that2his theory can be used for algorithmic construction of algebraic Gr¨obner bases. JosephM. Thomas [56, 57] generalized the Riquier-Janet theory to polynomially nonlinear systemsand showed how to decompose algebraic and differential systems into triangular subsystemswith disjoint solution sets. Thomas called these subsystems simple since their structurealleviates their algebraic analysis. The first algorithmization and implementation in Maple ofThomas’ approach for systems of algebraic and ordinary differential equations was achieved byDongming Wang [59, 60]. The further algorithmic development of Thomas decomposition foralgebraic systems and its full algorithmization for PDE systems, incorporating the involutivealgorithm for constructing Janet bases [17], together with an implementation in Maple, wasrealized in [6, 51, 21]. Thomas decomposition provides regular differential chains [29], whichallow to test membership to the radical differential ideal through differential Janet reduction.Related methods are the Rosenfeld-Gr¨obner algorithm [5] and the rif-algorithm [47]. Incontrast to regular differential chains generated by the Rosenfeld-Gr¨obner algorithm, theThomas decomposition and the rif-algorithm provide partitions of the solution sets. Becauseof the last property these decompositions lend themselves to the s-consistency analysis.However, the concept of s-consistency, as it was introduced in [23, 18, 25], is applicableto Cartesian grids only. Its generalization to more general regular grids, whose grid spacingsmay be pairwise different, requires certain modifications and extensions. These modificationsand extensions are presented in the given paper and illustrated by a number of examples thatinclude incompressible Navier-Stokes equations and overdetermined PDE systems. For someexamples we performed not only symbolic but also numeric analysis.The paper is organized as follows. Section 2 contains a description of differential Thomasdecomposition into simple differential systems which, in addition to equations, may includeinequations. The illustrative Example 1 is given and the fundamental property of simplesystems used in the s-consistency analysis is formulated in Proposition 1. In Section 3 weconsider finite difference approximations to the PDE system (1) on a regular grid (6) anddefine the differential and difference ideals generated by PDE and FDA, respectively. Theconcept of difference Gr¨obner basis together with related definitions and the simplest formof difference Buchberger algorithm are outlined in Section 4. Then, in Section 5 we givethe definition of s -consistency of FDA to PDE for the grid (6). In addition, we present thecriterion of s-consistency in terms of difference Gr¨obner bases. As an example of application ofthis criterion, we construct and analyze certain s-consistent FDA to the incompressible Navier-Stokes equations (Section 6). In Section 7 we define simple and quasi-simple difference systemsand describe the algorithm of Thomas decomposition into quasi-simple difference systems .We prove correctness and termination of the algorithm and show how it provides the fullyalgorithmic check of s-consistency. Two examples (Examples 5 and 6) of quasi-linear PDEand different FDA are analyzed with respect to their s-consistency in Section 8. Concludingremarks are given in Section 9. The proof of admissibility of the difference monomial orderingused in Example 6 is postponed to Appendix A.
Let K be the field of (complex) meromorphic functions on some connected open subset D of C n with coordinates x , . . . , x n . A system of polynomial partial differential equations andinequations (or differential system for brevity) for m unknown functions u (1) , u (2) , . . . , u ( m ) of x , . . . , x n is given by f = 0 , f = 0 , . . . f p = 0 , g (cid:54) = 0 , g (cid:54) = 0 , . . . g q (cid:54) = 0 , (2)3here p ∈ Z ≥ , q ∈ Z ≥ , and all f i and g j are elements of the differential polynomial ring K{ u } , endowed with the set ∆ = { ∂ , . . . , ∂ n } of commuting derivations. Most commonly,a solution of (2) is an m -tuple ( φ , . . . , φ m ) of locally analytic functions on D which satisfyevery equation and inequation of (2). Around any point z of the domain each function φ i hasan expansion as convergent power series (cid:88) k ∈ ( Z ≥ ) n c k ( x − z ) k k ! , ( x − z ) k = ( x − z ) k · · · ( x n − z n ) k n , k ! = k ! · · · k n ! , with certain coefficients c k ∈ C .Given a differential system (2) the determination of (even just formal) power series solu-tions around some point z is in general a non-trivial task, because integrability conditions needto be taken into account and the system of simultaneous algebraic equations and inequationsfor the coefficients c k requires splitting into different cases due to nonlinearity. Example 1 (cf. also Ex. 2.1.46 in [51]) . For simplicity we choose z = (0 ,
0) for investigatingformal power series solutions of the overdetermined system of quasilinear PDE f := u x − u = 0 ,f := u y,y − u = 0 , for u = u ( x, y ) , D = C . (3)Each of the two differential equations by itself is equivalent to (cid:88) ( k ,k ) ∈ ( Z ≥ ) c ( k +1 ,k ) x k y k k ! k ! − (cid:88) ( k ,k ) ∈ ( Z ≥ ) c ( k ,k ) x k y k k ! k ! = 0and (cid:88) ( k ,k ) ∈ ( Z ≥ ) c ( k ,k +2) x k y k k ! k ! − (cid:88) ( k ,k ) ∈ ( Z ≥ ) c ( k ,k ) x k y k k ! k ! = 0 , respectively, and hence equivalent to a system of algebraic recurrence equations c , = c , ,c , = 2 c , c , ,c , = 2 c , c , ,c , = 2 c , c , + 2 c , c , , ... c , = c , ,c , = 3 c , c , ,c , = 3 c , c , ,c , = 3 c , c , + 6 c , c , c , , ...respectively. However, the integrability condition0 = ∂ y f − ∂ x f = 3 u u x − u u y,y − u y (mod f = 0 , f = 0)= 3 u − u − u y = ( u − √ u y ) ( u + √ u y )reveals that further conditions on c ( k ,k ) are implied when f = 0, f = 0 is considered as asystem. Taking the above factorization into account, we obtain c , = c , / √ ,c , = √ c , c , ,c , = √ c , c , ,c , = √ c , c , + c , c , ) , ... ∨ c , = − c , / √ ,c , = −√ c , c , ,c , = −√ c , c , ,c , = −√ c , c , + c , c , ) , ...4he method of Thomas decomposition does not require polynomial factorization. If theabove factorization is ignored, the newly-discovered consequence 2 u y − u = 0 translates intoalgebraic conditions on the Taylor coefficients c ( k ,k ) as follows: c , − c , = 0 , c , c , − c , c , = 0 , c , c , − c , c , = 0 , c , c , + c , ) − c , ( c , c , + 3 c , ) = 0 , ... (4)In this example the process of finding integrability conditions is complete because furthercross-derivatives reduce to zero modulo the previous equations. The given system does notimpose any conditions on the Taylor coefficient c , , whose value can be chosen arbitrarily,and the possible values of all other Taylor coefficients are determined by the above algebraicequations. Taking the total order of differentiation into account, a systematic way of solvingthese algebraic equations is to solve each equation for the underlined variable. In order toensure both square-freeness of the first polynomial equation for c , in (4) and solvability ofall subsequent equations, a case distinction whether u ( x, y ) is the zero function or not is alsomade. Therefore, a Thomas decomposition of system (3) is u x − u = 0 , u y − u = 0 ,u (cid:54) = 0 , ∨ u = 0 . If the additional effort in factorizing the integrability condition 2 u y − u = 0 is spent, aThomas decomposition of the same system is also given by u x − u = 0 ,u y − u / √ , ∨ u x − u = 0 ,u y + u / √ . Note that even when no polynomial factorization is performed, a Thomas decomposition ofa PDE system is not uniquely determined in general.Computing a Thomas decomposition of a differential system is a finite process whichconstructs a generating set of all integrability conditions systematically and performs casesplittings, if necessary, so as to obtain a generating set of recurrence relations for c k around ageneric center of expansion. This process is steered by a total order (cid:23) on the set of symbolsrepresenting derivatives of unknown functions:Mon(∆) u := { ∂ k u ( α ) = ∂ k · · · ∂ k n n u ( α ) | ≤ α ≤ m, k ∈ ( Z ≥ ) n } . (5)(We shall mainly be working with the strict total order (cid:31) associated with (cid:23) .) Definition 1.
Let R = K{ u } be the differential polynomial ring and f ∈ R \ K .1. A ranking (cid:31) on R is a total order on Mon(∆) u such that for all 1 ≤ α ≤ m and all k (cid:54) = we have ∂ k u ( α ) (cid:31) u ( α ) , and such that ∂ k u ( α ) (cid:31) ∂ k u ( α ) implies ∂ k + k (cid:48) u ( α ) (cid:31) ∂ k + k (cid:48) u ( α ) for all k (cid:48) ∈ ( Z ≥ ) n . A ranking (cid:31) is said to be orderly if | k | > | k | implies ∂ k u ( α ) (cid:31) ∂ k u ( α ) for any α . 5. The leader ld( f ) of the differential polynomial f with respect to a ranking (cid:31) is thehighest ranked derivative in Mon(∆) u that effectively occurs in f .3. The coefficient of the highest power of ld( f ) in f is the initial of f , denoted by init( f ).It is itself a differential polynomial in derivatives that are ranked lower than ld( f ) withrespect to (cid:31) .4. The discriminant disc( f ) is the discriminant of f as a polynomial in ld( f ).5. The separant of f is the differential polynomial sep( f ) := ∂f /∂ ld( f ). Example 2. If x , y are the independent variables, u the dependent variable, and f = u y u x,y + u ∈ K{ u } , then, with respect to any orderly ranking (cid:31) on K{ u } , we have ld( f ) = u x,y andinit( f ) = u y and sep( f ) = 2 u y u x,y . Note that, generally, the separant of f is the initial ofany proper derivative of f , e.g., ∂ x f = 2 u y u x,y u x,x,y + u x,y + 5 u u x has leader u x,x,y andinitial 2 u y u x,y .The determination of all integrability conditions of a system of polynomially nonlinearPDE is facilitated by a combination of Euclid’s algorithm with case distinctions and comple-tion to involution as performed by Janet’s algorithm. Before recalling the latter ingredientwe outline the former aspect. In what follows we assume that a ranking (cid:31) on R is fixed.Note that any linear combination with coefficients in R of (the left hand sides of) equations f = 0, . . . , f p = 0 and their partial derivatives in a differential system is a consequence of thatsystem, and these consequences form a differential ideal of R . Every differential polynomial f ∈ R \ K is considered as a univariate polynomial in ld( f ) whose coefficients are themselvesunivariate polynomials in their leaders. In this way an algebraic and a differential reduction are defined for all pairs ( f , f ) ∈ ( R \ K ) , producing a differential polynomial f that iseither in K or has a leader that is ranked lower than ld( f ) with respect to (cid:31) .a) If ld( f ) = ld( f ) =: v and d := deg v ( f ) ≥ d := deg v ( f ), then let f = c f − c v d − d f , where c is a suitable power of init( f ) and c ∈ R such that the d -th power of v cancelsin f .b) If ld( f ) = ∂ k ld( f ) =: v for ∂ k ∈ Mon(∆) u , k (cid:54) = , and d := deg v ( f ), then let f = c f − c v d − ∂ k f , where c = sep( f ) and c ∈ R such that the d -th power of v cancels in f .Note that f is an element of the differential ideal containing f and f in any case.If f = 0, f = 0 are two equations in a differential system, then replacing f = 0by f = 0 is supposed to not alter the solution set of the system. This is ensured if thedifferential polynomial c does not vanish on the solution set of the system. Note that c ischosen as a power of init( f ) or sep( f ). If Euclid’s algorithm considers separately the casesobtained by adding the inequation init( f ) (cid:54) = 0 (resp. sep( f ) (cid:54) = 0) or the equation init( f ) = 0(resp. sep( f ) = 0) to the system, the above replacement of f = 0 by f = 0 is justified withthe imposed inequation, and the solution sets corresponding to the branches of computationdefine a partition of the solution set of the original system.Ignoring that the indeterminates represent unknown functions of a PDE system, Euclid’salgorithm deals with a system S of algebraic equations, say in, y , . . . , y r , totally ordered6y the fixed ranking. The solution set Sol( S ) in C r of that algebraic system is investigatedwith respect to a sequence of projections from C r to affine subspaces which corresponds tothe ordering, say, y (cid:31) y (cid:31) . . . (cid:31) y r , of the indeterminates: π : C r −→ C r − : ( a , a , . . . , a r ) (cid:55)−→ ( a , a , a , . . . , a r ) ,π : C r −→ C r − : ( a , a , . . . , a r ) (cid:55)−→ ( a , a , . . . , a r ) , ... ... π r − : C r −→ C : ( a , a , . . . , a r ) (cid:55)−→ a r . Euclid’s algorithm, performing case distinctions with regard to the vanishing of initials init( f )and discriminants disc( f ) of (non-constant) polynomials f , produces a finite collection ofalgebraic systems having the following property. Definition 2.
Let S = { f = 0 , . . . , f p = 0 , g (cid:54) = 0 , . . . , g q (cid:54) = 0 } be an algebraic system,i.e., f i , g j ∈ K [ y , . . . , y r ]. Then S is said to be simple if the following four conditions aresatisfied.1. None of f , . . . , f p , g , . . . , g q is constant.2. The leaders of f , . . . , f p , g , . . . , g q are pairwise distinct.3. For every h ∈ { f , . . . , f p , g , . . . , g q } , if ld( h ) = y k , then the equation init( h ) = 0 hasno solution in π k (Sol( S )).4. For every h ∈ { f , . . . , f p , g , . . . , g q } , if ld( h ) = y k , then the equation disc( h ) = 0 hasno solution in π k (Sol( S )).(Note that in 3. and 4. we have init( h ) , disc( h ) ∈ K [ y k +1 , . . . , y r ].) Definition 3.
An algebraic system S as in Definition 2 is said to be quasi-simple if condi-tions 1.–3. (but not necessarily 4.) are satisfied. (Such systems are also called regular , cf. [59,p. 107], [29], [31], [37].)Our strategy for handling integrability conditions builds on Janet division . Note first thatthe leader of the derivative of an equation f = 0 is the corresponding derivative of ld( f ).Hence, for each α ∈ { , . . . , m } , the monomials ∂ k , k ∈ ( Z ≥ ) n , for which ∂ k u ( α ) is theleader of a consequence of a differential system, form a set of monomials that is closed undermultiplication by ∂ , . . . , ∂ n .Suppose a set of monomials is closed under multiplication by the elements of a certainsubset µ of ∆ = { ∂ , . . . , ∂ n } . If that set of monomials consists of all such multiples of a singlemonomial, then we call the set a cone . Let M be a finite set of monomials. Janet divisionassigns to each m ∈ M a set µ ( m, M ) ⊆ ∆ of multiplicative variables so as to decompose theset of all multiples of M into disjoint cones. Denoting by Mon( µ ) the set of all monomials inthe elements of µ , we have (cid:91) m ∈ M Mon(∆) m ⊇ (cid:93) m ∈ M Mon( µ ( m, M )) m . In case of equality the set M is said to be Janet complete .In our context we call the multiplicative variables admissible derivations .7 efinition 4. Let M be a finite set of monomials in ∂ , . . . , ∂ n . For 1 ≤ j ≤ n we let ∂ j bean admissible derivation for ∂ i · · · ∂ i n n ∈ M if and only if i j = max { k j | ∂ k · · · ∂ k n n ∈ M with k = i , k = i , . . . , k j − = i j − } . Example 3.
Let M = { ∂ ∂ , ∂ ∂ , ∂ ∂ , ∂ ∂ } . These four monomials are assigned thesets of admissible derivations { ∂ , ∂ , ∂ } , { ∂ , ∂ } , { ∂ , ∂ } and { ∂ } , respectively.We extend Janet division as well as the notion of Janet completeness from finite sets ofmonomials to finite sets { f , . . . , f p } of differential polynomials in R \ K by assigning f i theset of admissible derivations µ i := µ ( θ i , { θ , . . . , θ p } ), where θ i ∈ Mon( { ∂ , . . . , ∂ n } ) is suchthat ld( f i ) = θ i u ( α i ) for a certain α i .By restricting the differential reduction process introduced in b) above to reduction stepsfor which ∂ k is a monomial in admissible derivations for f , we obtain the Janet reductionprocess. The remainder of a differential polynomial f modulo { f , . . . , f p } , or modulo T = { ( f , µ ) , . . . , ( f p , µ p ) } , is called the Janet normal form of f modulo T and is denoted byNF( f, T, (cid:31) ). Definition 5.
Let T = { ( f , µ ) , . . . , ( f p , µ p ) } be Janet complete. Then the differentialsystem { f = 0 , . . . , f p = 0 } , or T , is said to be passive ifNF( ∂f i , T, (cid:31) ) = 0 for all ∂ ∈ µ i = { ∂ , . . . , ∂ n } \ µ i , i = 1 , . . . , p . A suitable combination of Euclid’s algorithm with case distinctions and differential reduc-tions of differential polynomials that are obtained by applying non-admissible derivations de-fines a process that returns a Thomas decomposition in finitely many steps [51, Thm. 2.2.57],[6, Sect. 3.4], namely, a finite collection of differential systems, whose solution sets partitionthe solution set of the original differential system, and such that each output system has thefollowing property.
Definition 6.
A differential system S as in (2) is said to be simple if the following threeconditions hold.1. S is simple as an algebraic system (in the finitely many indeterminates occurring in it,ordered by the ranking (cid:31) ; cf. Definition 2).2. { f = 0 , . . . , f p = 0 } is passive (cf. Definition 5).3. The left hand sides g , . . . , g q are Janet reduced (i.e., in Janet normal form) modulothe equations { f = 0 , . . . , f p = 0 } .A simple differential system S allows to decide, by differential reduction, whether or nota given equation f = 0, where f ∈ R , is a consequence of S . Proposition 1 ([51], Prop. 2.2.50) . Let S be a simple differential system, defined over thedifferential polynomial ring R , and let E be the differential ideal of R which is generated by f , . . . , f p . Moreover, let q be the product of the initials and separants of all f , . . . , f p .Then the differential ideal E : q ∞ := { f ∈ R | q r f ∈ E for some r ∈ Z ≥ } is radical. Given f ∈ R , we have f ∈ E : q ∞ if and only if the Janet normal form of f modulo { f , . . . , f p } is zero. Difference approximations to PDE systems
To approximate the differential system (1) by a difference system we shall consider a regularcomputational grid (mesh) as the set of points { ( z + k h , . . . , z n + k n h n ) | k i ∈ Z } , (6)where ( z , . . . , z n ) ∈ R n and 0 < h i ∈ R are fixed. Definition 7.
A vector function ˜ u = { ˜ u (1) , . . . , ˜ u ( m ) } which assigns to each grid node a valueis called grid vector function . We shall denote such function by ˜u k ,...,k n := ˜u ( z + k h , . . . , z n + k n h n ) . From now on we shall consider h , . . . , h n as parameters and denote by h := { h , . . . , h n } and Mon( h ) := (cid:110) n (cid:89) i =1 h µ i i (cid:12)(cid:12)(cid:12) µ ∈ ( Z ≥ ) n (cid:111) the set of mesh steps (grid spacings) and the monoid of monomials generated by the elementsin h , respectively. The total degree of an element m ∈ Mon( h ) will be denoted by deg( m ).We assume that coefficients of the differential polynomials in F (cf. (1)) do not vanish inthe grid points. The coefficients on the grid as rational functions in { a , h } are elements ofthe difference field [38] with differences { σ , . . . , σ n , σ − , . . . , σ − n } acting on a grid function˜ u ( α ) k ,...,k n as the shift operators σ ± j ˜ u ( α ) k ,...,k j ,...,k n = ˜ u ( α ) k ,...,k j ± ,...,k n , α ∈ { , . . . , m } , j ∈ { , . . . , n } . (7)Let Mon(Σ) be the free commutative semigroup generated by Σ = { σ , . . . , σ n } ,Mon(Σ) := { σ i · · · σ i n n | i , . . . , i n ∈ Z ≥ } , (8)˜ K = Q ( a , h ), and ˜ R the ring of difference polynomials over ˜ K . The elements in ˜ R are polyno-mials in the difference indeterminates ˜ u ( α ) ( α = 1 , . . . , m ) and in their shifts σ i · · · σ i n n ˜ u ( α ) , i , . . . , i n ∈ Z , with coefficients in ˜ K . Remark 1.
Since the shift operators σ j admit the (formal) power series expansion σ j = (cid:88) k ≥ h kj k ! ∂ kj , σ − j = (cid:88) k ≥ ( − k h kj k ! ∂ kj , a difference polynomial f (˜ u ) ∈ ˜ R admits the Taylor expansion around a grid point (6).The standard technique to obtain FDA to the PDE system (1) is to replace the derivativesoccurring in (1) by finite differences. In order to use the method of difference Gr¨obner bases(Section 4) or/and difference Thomas decomposition (Section 7) one has to apply appropriatepower products of the right-shift operators (7) to remove negative shifts in indices which maybe introduced from expressions like ∂ j u ( α ) = ˜ u ( α ) k ,...,k j ,...,k n − σ − j ˜ u ( α ) k ,...,k j ,...,k n h j + O ( h j ) ,∂ j u ( α ) = σ j ˜ u ( α ) k ,...,k j ,...,k n − σ − j ˜ u ( α ) k ,...,k j ,...,k n h j + O ( h j ) . In the sequel we shall consider discretization of (1) as a finite set of difference polynomials˜ f = · · · = ˜ f p = 0 , ˜ F := { ˜ f , . . . , ˜ f p } ⊂ ˜ R . (9)9 efinition 8. The differential (resp. difference) ideal generated by a polynomial set F (resp. ˜ F ) , denoted by I := [ F ] (resp. ˜ I := [ ˜ F ]), is the smallest subset of R (resp. ˜ R ) containing F (resp. ˜ F ) and satisfying( ∀ ∂ i ∈ { ∂ , . . . , ∂ n } ) ( ∀ a, b ∈ I ) ( ∀ c ∈ R ) [ a + b ∈ I , a · c ∈ I , ∂ i a ∈ I ]and, respectively,( ∀ σ i ∈ { σ , . . . , σ n } ) ( ∀ ˜ a, ˜ b ∈ ˜ I ) ( ∀ ˜ c ∈ ˜ R ) [ ˜ a + ˜ b ∈ ˜ I , ˜ a · ˜ c ∈ ˜ I , σ i ˜ a ∈ ˜ I ] . Let
I ⊂ R be a differential ideal. Then the set √I := { p ∈ R | p k ∈ I , k ∈ N > } (10)is a differential ideal.If I = √I , then I is called radical or perfect differential ideal . Given F ⊂ R , the radicaldifferential ideal generated by F , denoted by (cid:74) F (cid:75) , is the smallest radical differential ideal of R containing F .In the difference case, the radical (cid:112) ˜ I of ˜ I is defined similarly to Eq. (10). However, thenotion of perfect difference ideal is significantly distinct from that of perfect differential idealin differential algebra [50]. Definition 9.
The perfect difference ideal [38] generated by a set ˜ F ⊂ ˜ R , denoted by (cid:74) ˜ F (cid:75) ,is the smallest difference ideal of ˜ R containing ˜ F and such that for any ˜ f ∈ ˜ R , θ , . . . , θ r ∈ Mon(Σ) and k , . . . , k r ∈ Z ≥ we have( θ ˜ f ) k · · · ( θ r ˜ f ) k r ∈ (cid:74) ˜ F (cid:75) = ⇒ ˜ f ∈ (cid:74) ˜ F (cid:75) . It is clear that [ ˜ F ] ⊆ (cid:112) ˜ F ⊆ (cid:74) ˜ F (cid:75) . In difference algebra perfect ideals are analoguesof radical ideals in commutative [13] and differential algebra [50, 29]. In particular, thedifference Hilberts Nullstellensatz is formulated in terms of perfect difference ideals (cf. [12],Ch. 4, Thm. 4 and [38], Thm. 2.6.5). For this reason we give the following definition. Definition 10.
We shall say that a differential (resp. difference) polynomial f ∈ R (resp.˜ f ∈ ˜ R ) is a differential-algebraic (resp. difference-algebraic ) consequence of (1) (resp. of (9))if f (resp. ˜ f ) is an element of the perfect differential (resp. difference) ideal generated by (1)(resp. (9)).Some recent results on the relation between the difference Hilberts Nullstellensatz andsolvability are presented in [45, 46]. The notion of difference Gr¨obner basis was introduced and studied in [18, 36, 22]. It is adifference analogue of the notion of differential standard basis introduced in [43], where afinite standard basis is called
Gr¨obner basis . In this paper we prefer to use the approach todifference Gr¨obner bases suggested in [18].
Definition 11. A ranking on ˜ R is defined in the same way as in Definition 1 by replacingthe action of ∂ i by the action of σ i and ∆ by Σ.10 efinition 12. A total ordering (cid:65) on the set of difference monomials M := (cid:110) ( θ ˜ u (1) ) i · · · ( θ m ˜ u ( m ) ) i m (cid:12)(cid:12)(cid:12) θ j ∈ Σ , i j ∈ Z ≥ , ≤ j ≤ m (cid:111) is an admissible (difference) monomial ordering if it extends a ranking and satisfies( a ) ( ∀ ˜ t ∈ M \ { } ) [˜ t (cid:65) , ( b ) ( ∀ θ ∈ Σ) ( ∀ ˜ t, ˜ v, ˜ w ∈ M ) [ ˜ v (cid:65) ˜ w ⇐⇒ ˜ t · θ ◦ ˜ v (cid:65) ˜ t · θ ◦ ˜ w ] . For examples of admissible monomial orderings we refer to Appendix A.Given an admissible ordering (cid:65) , every difference polynomial ˜ f has the leading monomial lm( ˜ f ) ∈ M with leading coefficient lc( ˜ f ). In what follows every difference polynomial is tobe normalized (i.e., monic) by division by its leading coefficient. Definition 13.
If for v, w ∈ M the equality w = t · θ ◦ v holds with θ ∈ Mon(Σ) and t ∈ M we shall say that v divides w and write v | w . It is easy to see that this divisibility relationyields a partial order. Definition 14.
Given a difference ideal ˜ I and an admissible monomial ordering (cid:65) , a subset˜ G ⊂ ˜ I is a (difference) Gr¨obner basis for ˜ I if [ ˜ G ] = ˜ I and( ∀ ˜ f ∈ I )( ∃ ˜ g ∈ ˜ G ) [ lm(˜ g ) | lm( ˜ f ) ] . Definition 15.
A polynomial ˜ p ∈ ˜ R is said to be head reducible modulo ˜ q ∈ ˜ R to ˜ r if˜ r = ˜ p − m · θ ◦ ˜ q and m ∈ M , θ ∈ Mon(Σ) are such that lm(˜ p ) = m · θ ◦ lm(˜ q ). In thiscase transformation from ˜ p to ˜ r is an elementary reduction , denoted by ˜ p −→ ˜ q ˜ r . Given a set˜ F ⊂ ˜ R , ˜ p is head reducible modulo ˜ F (denotation: ˜ p −→ ˜ F ) if there is ˜ f ∈ ˜ F such that ˜ p is headreducible modulo ˜ f . A polynomial ˜ p is head reducible to ˜ r modulo ˜ F if there is a finite chainof elementary reductions ˜ p −→ ˜ F ˜ p −→ ˜ F ˜ p −→ ˜ F · · · −→ ˜ F ˜ r . (11)If no monomial in ˜ r from (11) is head reducible modulo ˜ F , then ˜ r is in head normal formmodulo ˜ F and we write ˜ r = HNF(˜ p, ˜ F , (cid:65) ). Similarly, one can define tail reduction and (full)normal form (denotation: NF(˜ p, ˜ F , (cid:65) ) . A polynomial set ˜ F with more than one element is (head) interreduced if ( ∀ ˜ f ∈ ˜ F ) [ ˜ f = (H)NF( ˜ f , ˜ F \ { ˜ f } , (cid:65) ) ] . (12)Admissibility of (cid:65) , as in commutative algebra, provides termination of the chain (11) forany ˜ p and ˜ F . Then (H)NF(˜ p, ˜ F , (cid:65) ) can be computed by the difference version of a multivariatepolynomial division algorithm [7, 13]. If ˜ G is a Gr¨obner basis of [ ˜ G ], then from Definitions 14and 15 it follows ˜ f ∈ [ ˜ G ] ⇐⇒ (H)NF( ˜ f , ˜ G, (cid:65) ) = 0 . Thus, if an ideal has a finite Gr¨obner basis, then its construction solves the ideal membershipproblem in the same way as in commutative [7, 13] and differential [43, 63] algebra. Thealgorithmic characterization of difference Gr¨obner bases and their construction in differencepolynomial rings employ difference S -polynomials.11 efinition 16. Given an admissible ordering and monic difference polynomials ˜ p and ˜ q , apolynomial S (˜ p, ˜ q ) := m · θ ◦ ˜ p − m · θ ◦ ˜ q is called S -polynomial associated to ˜ p and ˜ q if m · θ ◦ lm(˜ p ) = m · θ ◦ lm(˜ q ) with co-prime m · θ and m · θ (for ˜ p = ˜ q we shall say thatthe S -polynomial is associated with ˜ p ). Proposition 2.
Given a difference ideal ˜ I ⊂ ˜ R and an admissible monomial ordering (cid:65) , aset of polynomials ˜ G ⊂ ˜ I is a Gr¨obner basis of ˜ I if and only if (H)NF( S (˜ p, ˜ q ) , ˜ G, (cid:65) ) = 0 (13) for all S -polynomials associated with polynomials in ˜ G .Proof. This follows from Definitions 14, 15 and 16 in line with the standard proof of theanalogous theorem for Gr¨obner bases in commutative algebra [7, 13] and with the proof of asimilar theorem for standard bases in differential algebra [43].
Definition 17.
Given a system (9) of difference equations, the conditions (13) with ˜ p, ˜ q ∈ ˜ F are said to be the (Gr¨obner) passivity conditions for the system (9).Let ˜ I = [ ˜ F ] be a difference ideal generated by a finite set ˜ F ⊂ ˜ R of difference poly-nomials with non-negative shifts. Then for a fixed admissible monomial ordering the algo-rithm DifferenceGr¨obnerBasis given below, if it terminates, returns a Gr¨obner basis ˜ G of˜ I . The subalgorithm Interreduce invoked in line 9 performs mutual (head) interreductionof the elements in ˜ H and returns a set satisfying (12).Algorithm DifferenceGr¨obnerBasis is a difference analogue of the simplest version ofBuchberger’s algorithm (cf. [7, 13, 43]). Its correctness is provided by Theorem 2. Thealgorithm always terminates when the input polynomials are linear. If this is not the case,the algorithm may not terminate. This means that the repeat-until loop (lines 2–8) maybe infinite as in the differential case [43, 63]. One can improve the algorithm by taking intoaccount Buchberger’s criteria to avoid some useless zero reductions in line 5. The differencecriteria are similar to the differential ones [43].
Algorithm 1:
DifferenceGr¨obnerBasis
Input: ˜ F ⊂ ˜ R \ { } , a finite set of non-zero polynomials; (cid:65) , a monomial ordering Output: ˜ G , a (head) interreduced Gr¨obner basis of [ ˜ F ] ˜ G ← ˜ F repeat ˜ H ← ˜ G for S -polynomials ˜ s associated with elements in ˜ H do ˜ g ← (H)NF(˜ s, ˜ H, (cid:65) ) if ˜ g (cid:54) = 0 then ˜ G ← ˜ G ∪ { ˜ g } until ˜ G = ˜ H ˜ G ← Interreduce ( ˜ G ) return ˜ G Consistency
Let the PDE system (1) and its finite difference discretization (9) on the regular grid (6) begiven.
Definition 18.
We shall say that a difference equation ˜ f ( ˜u ) = 0, ˜ f ∈ ˜ R , implies the set ofdifferential equations F := { r k ( u ) = 0 | r k ∈ R , k = 1 , . . . , K, K ∈ N > } and we write ˜ f (cid:66) F for K > f (cid:66) r for K = 1, if Taylor expansion of ˜ f about a gridpoint, after clearing denominators containing the elements in h = ( h , . . . , h n ) by multiplyingby an appropriate q ( h ), yields q ( h ) · ˜ f ( ˜u ) = K (cid:88) k =1 m k ( h ) · r k ( u ) + O ( d + 1) , (14)where ( ∀ k ) [ m k ( h ) ∈ Mon( h ) , deg( m k ( h )) = d ] for some d ∈ Z > and O ( d + 1) denotesterms whose total degree in h , . . . , h n is larger than d . Definition 19.
A difference system (9) is weakly consistent or w-consistent with PDE (1) if( ∀ j ∈ { , . . . , p } ) [ ˜ f j (cid:66) f j ] . (15)It is clear that if one considers a single PDE f ( u ) = 0 and a difference equation ˜ f ( ˜u ) = 0with implication ˜ f (cid:66) f , i.e., q ( h ) · ˜ f ( ˜u ) = m ( h ) · f ( u ) + O (deg( m ) + 1) , m ∈ Mon( h ) , (16)then it is always possible to redefine ˜ f (cid:48) := q ( h ) ˜ f /m ( h ), and there is a limit (cf. Proposition 3) h →
0, i.e. ( ∀ i )[ h i → f (cid:48) −−−→ h → f + O ( h ) , taking Remark 1 into account.The condition (15) means that Eqs. (9) approximate Eqs. (1) and by this reason we callEqs. (9) finite difference approximation (FDA) to Eqs. (1).
Remark 2.
Given ˜ f (˜ u ), computation of f ( u ) is straightforward and has been implementedas routine ContinuousLimit in the Maple package
LDA [24] (Linear Difference Algebra).
Example 4. ([54], Ex. 1.4.2) We consider the one-way wave equation f := u t + a u x = 0 , u = u ( t, x ) , a = const , (17)where a is a constant, t represents time, and x represents the spatial variable. For the classicalLax-Friedrichs discretization the difference form of Eq. (17) for the grid with t n +1 − t n = h , x m +1 − x m = h is given by˜ f := ˜ u n +1 m − (cid:0) ˜ u nm +1 + ˜ u nm − (cid:1) h + a ˜ u nm +1 − ˜ u nm − h = 0 , (18)where ˜ u nm := ˜ u ( nh , mh ) is a smooth grid function. The Taylor expansion of (18) aroundthe point ( t = nh , x = mh ) reads˜ f −−−−−→ h ,h → ( u t + a u x ) h + 12 h u tt − h u xx + 16 a h h u xxx + · · · . So ˜ f (cid:66) f as h , h →
0, and ˜ f /h → f as h → h /h → Definition 20.
An FDA (9) to a PDE system (1) with p > strongly consistent or s-consistent if ( ∀ ˜ f ∈ (cid:74) ˜ F (cid:75) ) ( ∃ F ⊂ (cid:74) F (cid:75) ) [ ˜ f (cid:66) F ] . (19)In [18, Definition 12] we defined s-consistent FDA to PDE for Cartesian grids h = h = · · · = h n = h as the ones satisfying the condition( ∀ ˜ f ∈ (cid:74) ˜ F (cid:75) ) ( ∃ f ∈ (cid:74) F (cid:75) ) [ ˜ f (cid:66) f ]that corresponds to the Taylor expansion˜ f ( ˜u ) = h k f ( u ) + O ( h k +1 )of the form (16) and to the consistency condition ˜ f −−−→ h → f with ˜ f = ˜ f /h k .The s-consistency property (19) implies the existence of a limit h → h , is a differential-algebraic consequence of (1). Proposition 3.
Let FDA (9) to PDE system (1) be defined on a regular solution grid (6) . Ifthe perfect difference ideal generated by the set ˜ F satisfies the condition (19) of s-consistency,then there is a limit h → such that ( ∀ ˜ f ∈ (cid:74) ˜ F (cid:75) ) ( ∃ f ∈ (cid:74) F (cid:75) , µ = O ( h ) , d ∈ Z ≥ ) (cid:34) ˜ fµ d −−−→ h → f (cid:35) . (20) Proof.
Let ˜ f ( ˜u ) ∈ (cid:74) ˜ F (cid:75) . Then the equation ˜ f ( ˜u ) = 0 is a difference-algebraic consequence ofsystem (9). Without loss of generality one may assume that the denominators in ˜ f containingelements in h have been cleared. Now we consider the following limit h → h i := a i µ , < a i ∈ Q , µ ∈ R > , i = 1 , . . . , n , µ → . Then, from the Taylor expansion (14) we obtain˜ f ( ˜u ) = µ d (cid:32) K (cid:88) k =1 m k ( a ) · r k ( u ) + O ( µ ) (cid:33) , where a := ( a , . . . , a n ) ∈ Q n and (cid:80) Kk =1 m k ( a ) · r k ( u ) ∈ (cid:74) F (cid:75) .Let F be a PDE system (1) and ˜ F be a w-consistent FDA (9). In practice, FDA (scheme)can be obtained from PDE by approximation of the partial derivatives occurring in PDE withappropriate finite differences. Another way of discretization is to apply the method suggestedin [20]. In either case, to verify the s-consistency condition (19), one has to reformulate thiscondition to make it algorithmic.The first step in this direction is to use the following statement.14 heorem 1. A difference approximation (9) to a differential system (1) is s-consistent if andonly if a Gr¨obner basis ˜ G ⊂ ˜ R of the difference ideal [ ˜ F ] satisfies ( ∀ ˜ g ∈ ˜ G ) ( ∃ G ⊂ (cid:74) F (cid:75) ) [ ˜ g (cid:66) G ] . (21) Proof.
For a Cartesian grid with h = h = · · · = h n the proof was given in [18, Thm. 3]. Itsextension to the grid (6) is straightforward.If the difference Gr¨obner basis ˜ G is finite and we can construct it in finitely many steps,then we can compute the Taylor expansion (14) of every element ˜ f ( ˜u ) ∈ ˜ G and obtain˜ fµ d −−−→ µ → f = K (cid:88) k =1 m k ( q ) · r k ( u ) . Furthermore, to check the radical ideal membership f ∈ (cid:74) F (cid:75) , one can compute algorith-mically a differential Thomas decomposition (cf. Section 2, [21, 6, 51]) of the PDE (1) (withrespect to any ranking (cid:31) ) into the subsystems S , . . . , S t with disjoint solution sets and usethe relation f ∈ (cid:74) F (cid:75) ⇐⇒ NF( f, S , (cid:31) ) = . . . = NF( f, S t , (cid:31) ) = 0 , (22)where NF( f, S i , (cid:31) ) denotes the Janet normal form of f modulo S i , as defined above (beforeDefinition 5, cf. also [51, Prop. 2.2.50]).Since the difference polynomial ring ˜ R is non-Noetherian [22], the computation of aGr¨obner basis ˜ G is not algorithmic. Hence, Algorithm 1 may not terminate (cf. also [18, 22]).However, instead one can apply the difference triangular decomposition as developed in Sec-tion 7.If the input PDE system is linear, then the condition (21) is algorithmic and can beverified by computing Gr¨obner bases of the ideals generated by F and ˜ F [23]. All relatedcomputations can be done by using the relevant routines of the Maple packages LDA [24]and
Janet [8].
The Navier-Stokes equations for a three-dimensional incompressible flow of constant viscositycan be written as f := ∂ u + ∂ v + ∂ w = 0 ,f := ∂ t u + u∂ u + v∂ u + w∂ u + ∂ p − ∇ u = 0 ,f := ∂ t v + u∂ v + v∂ v + w∂ v + ∂ p − ∇ v = 0 ,f := ∂ t w + u∂ w + v∂ w + w∂ w + ∂ p − ∇ w = 0 . (23)Here u ( x , t ) is the velocity vector u := ( u, v, w ), x := ( x , x , x ) is the vector of Cartesiancoordinates, p ( x , t ) is the pressure, ∂ i := ∂ x i ( i ∈ { , , } ), ∇ is the Laplace operator ∇ := ∂ + ∂ + ∂ and Re is the Reynolds number. Remark 3.
The PDE system (23) consist of the continuity equation or incompressibilitycondition f = 0 that represents conservation of mass and the momentum equations f i = 0 ( i ∈{ , , } ) that represent conservation of momentum. As a consequence of these fundamentalconservation laws, the system is invariant under a permutation of the coordinates providedthe corresponding permutation is applied to the components of u .15ne can rewrite system (23) equivalently as F := ∂ u + ∂ v + ∂ w = 0 ,F := ∂ t u + ∂ ( u ) + ∂ ( uv ) + ∂ ( uw ) + ∂ p − ∇ u = 0 ,F := ∂ t v + ∂ ( uv ) + ∂ ( v ) + ∂ ( vw ) + ∂ p − ∇ v = 0 ,F := ∂ t w + ∂ ( uw ) + ∂ ( vw ) + ∂ ( w ) + ∂ p − ∇ w = 0 , (24)where the nonlinear parts in the momentum equations are in divergence form or conservativeform , and modulo the continuity equation the systems (23) and (24) coincide.For the Navier-Stokes equations in the form (23), one can conveniently use vector notation,which has the advantage of brevity, and rewrite these equations as ∇ · u = 0 , ∂ t u + ( u · ∇ ) u + ∇ p − ∇ u = 0 , (25)where ∇ := ( ∂ , ∂ , ∂ ) is the nabla operator.Let (cid:31) be the ranking that compares first the monomials in the partial derivations ∂ t , ∂ , ∂ , ∂ (cf. Eq. (5)) with respect to the lexicographic ordering and then, in the case of equaldifferential monomials, compares differential indeterminates (dependent variables) as ∂ t (cid:31) ∂ (cid:31) ∂ (cid:31) ∂ and p (cid:31) u (cid:31) v (cid:31) w . (26)The (non-admissible) prolongation ∂ t F = ∇ · ∂ t u = 0 of the continuity equation and itsreduction modulo the vector momentum equation yields the pressure Poisson equation ∇ p + ∇ · ( u · ∇ ) u = 0 , (27)which is the integrability condition (cf. [53], p. 50) to Eqs. (25) and to (24) as well. Clearly,the differential system (25), (27) is passive and simple (cf. Definition 6). We mention thatthe arbitrariness of analytic solutions to the incompressible Navier-Stokes equations can alsobe represented by the differential counting polynomial ∞ (cid:96) + (cid:96) + (cid:96) +4 [35, Example 4.7].Eq. (27) can be expressed in terms of the continuity and momentum equations as F := ∂ F + ∂ F + ∂ F + 1Re (cid:0) ∂ F + ∂ F + ∂ F (cid:1) − ∂ t F . (28)It is significant that both Eqs. (27) and (28) preserve permutational symmetry in line withRemark 3.Now we consider the following class of FDA to (23) defined on the four-dimensional grid (6) D · ˜u = 0 , D t ˜u + ( ˜u · D ) ˜u + D ˜ p − ˜u = 0 , (29)where D t approximates ∂ t , D = ( D , D , D ) approximates ∇ and ˜∆ approximates ∇ . Itis clear that system (29) is w-consistent with Eqs. (25).As an example of such finite difference approximations on the grid (6), one can considerthe following one D t = σ t − τ , D i = σ i − σ − i h i , ˜∆ = (cid:88) j =1 σ i − σ − i h j , (30)where i ∈ { , , } and τ, h i ∈ R > . 16f one considers a difference analogue of Eq. (26) satisfying σ t (cid:31) σ (cid:31) σ (cid:31) σ and ˜ p (cid:31) ˜ u (cid:31) ˜ v (cid:31) ˜ w , (31)then completion of (29) to a passive form by Algorithm 1 is equivalent to enlargement of thissystem with the difference integrability condition( D · D ) ˜ p + D · ( ˜u · D ) ˜u − D · ˜∆ ˜u = 0 . (32)Eq. (32) approximates Eq. (27) and is obtained, in full analogy with the differential case, bythe prolongation D · D t ˜ u = 0 of the discrete continuity equation in (29) and its reductionmodulo the discrete vector momentum equation. Remark 4.
Because of equality D · ˜∆ ˜u = ˜∆ D · ˜u , the last term in Eq. (32) can be omittedif one considers a solution satisfying the discrete continuity equation D · ˜ u = 0. In such acase instead of Eq. (32) one can use( D · D ) ˜ p + D · ( ˜u · D ) ˜u = 0 . (33)The left-hand sides of Eqs. (29) and (33) form a difference Gr¨obner basis of the idealthey generate in Q (Re , h) { ˜u , ˜ p } . Hence, by Theorem 1, FDA (29), (33) is s-consistent withEqs. (23), (27).Similarly, we can approximate Eqs. (24) as follows: ˜ F := D ˜ u + D ˜ v + D ˜ w = 0 , ˜ F := D t ˜ u + D ˜ u + D (˜ u ˜ v ) + D (˜ u ˜ w ) + D ˜ p − ˜∆ ˜ u = 0 , ˜ F := D t ˜ v + D (˜ u ˜ v ) + D (˜ v ) + D (˜ v ˜ w ) + D ˜ p − ˜∆ ˜ v = 0 , ˜ F := D t ˜ w + D (˜ u ˜ w ) + D (˜ v ˜ w ) + D ( ˜ w ) + D ˜ p − ˜∆ ˜ w = 0 , (34)and complete this system to a passive form by performing the head reduction of D · D t ˜ u = 0modulo the momentum equations in (34). As a result, we obtain another discretization ofEq. (28): ˜ F := (cid:0) D · D (cid:1) ˜ p + (cid:0) D · D (cid:1) ˜ u + (cid:0) D · D (cid:1) ˜ v + (cid:0) D · D (cid:1) ˜ w +2 (cid:2)(cid:0) D · D (cid:1) (˜ u ˜ v ) + (cid:0) D · D (cid:1) (˜ u ˜ w ) + (cid:0) D · D (cid:1) (˜ v ˜ w ) (cid:3) − (cid:2)(cid:0) D · ˜∆ (cid:1) ˜ u + (cid:0) D · ˜∆ (cid:1) ˜ v + (cid:0) D · ˜∆ (cid:1) ˜ w (cid:3) = 0 . (35)We emphasize that the right-hand sides of Eqs. (34) and (35) are tale redundant modulo ˜ F .However, we prefer to use this redundant form since it inherits the permutational symmetryof Eqs. (24) and their divergence (conservative) form. Proposition 4.
FDA ˜ F := { ˜ F , ˜ F , ˜ F , ˜ F , ˜ F } to F := { F , F , F , F , F } is s-consistent.Proof. By inspection of the leading terms, it is easy to see that ˜ F is a head reduced differenceGr¨obner basis of [ ˜ F ] = (cid:74) ˜ F (cid:75) , and Theorem 1 implies the s-consistency.To compare different FDA to Eqs. (23), we compare their numerical behavior with theexact non-stationary two-dimensional solution [32] originally derived by Taylor [55] in hisstudy of decaying vortex flow. This solution is widely used as a benchmark for numerical17olving of Navier-Stokes equations (see, for example, [42]) and we have already used it in [1]and [2]. u = − e − t Re cos( x ) sin( y ) ,v = e − t Re sin( x ) cos( y ) ,p = − e − t Re (cos(2 x ) + cos(2 y )) . (36)We consider here four difference approximations to the two-dimensional form of Eqs. (25) and(27) with the grid functions˜ u nj,k := ˜ u ( jh, kh, nτ ) , ˜ p nj,k := ˜ p ( jh, kh, nτ ) , ( j, k, n ) ∈ Z , (37)and the following approximations of partial derivatives D t = σ t − τ , D i = σ i − σ − i h , ˜∆ = σ + σ − σ − + σ − h , (38)where i ∈ { , } and τ, h ∈ R > .FDA1 [2] ˜ F (1)0 := D ˜ u + D ˜ v = 0 , ˜ F (1)1 := D t ˜ u + D ˜ u + D (˜ u ˜ v ) + D ˜ p − ˜∆ ˜ u = 0 , ˜ F (1)2 := D t ˜ v + D (˜ u ˜ v ) + D (˜ v ) + D ˜ p − ˜∆ ˜ v = 0 , ˜ F (1)3 := (cid:0) D · D + D · D (cid:1) ˜ p + (cid:0) D · D (cid:1) ˜ u + (cid:0) D · D (cid:1) ˜ v +2 (cid:0) D · D (cid:1) (˜ u ˜ v ) − (cid:2)(cid:0) D · ˜∆ (cid:1) ˜ u + (cid:0) D · ˜∆ (cid:1) ˜ v (cid:3) = 0 . (39)FDA2 [2] ˜ F (2)0 := D ˜ u + D ˜ v = 0 , ˜ F (2)1 := D t ˜ u + ˜ u D ˜ u + ˜ v D ˜ u + D ˜ p − ˜∆ ˜ u = 0 , ˜ F (2)2 := D t ˜ v + ˜ u D ˜ v + ˜ v D ˜ v + D ˜ p − ˜∆ ˜ v = 0 , ˜ F (2)3 := ˜∆ ˜ p + ( D ˜ u ) + 2 ( D ˜ v )( D ˜ u ) + ( D ˜ v ) = 0 . (40)FDA3 ˜ F (3)0 := D ˜ u + D ˜ v = 0 , ˜ F (3)1 := D t ˜ u + D ˜ u + D (˜ u ˜ v ) + D ˜ p − ˜∆ ˜ u = 0 , ˜ F (3)2 := D t ˜ v + D (˜ u ˜ v ) + D (˜ v ) + D ˜ p − ˜∆ ˜ v = 0 , ˜ F (3)3 := (cid:0) D · D + D · D (cid:1) ˜ p + (cid:0) D · D (cid:1) ˜ u + (cid:0) D · D (cid:1) ˜ v +2 (cid:0) D · D (cid:1) (˜ u ˜ v ) = 0 . (41)FDA4 ˜ F (4)0 := D ˜ u + D ˜ v = 0 , ˜ F (4)1 := D t ˜ u + ˜ u D ˜ u + ˜ v D ˜ u + D ˜ p − ˜∆ ˜ u = 0 , ˜ F (4)2 := D t ˜ v + ˜ u D ˜ v + ˜ v D ˜ v + D ˜ p − ˜∆ ˜ v = 0 , ˜ F (4)3 := (cid:0) D · D + D · D (cid:1) ˜ p + D (cid:0) ˜ u D ˜ u (cid:1) + D (cid:0) ˜ v D ˜ u (cid:1) + D (cid:0) ˜ v D ˜ v (cid:1) − (cid:2)(cid:0) D · ˜∆ (cid:1) ˜ u + (cid:0) D · ˜∆ (cid:1) ˜ v (cid:3) = 0 . (42)18ll these FDA are explicit as difference schemes, w-consistent and they inherit permu-tational symmetry of the Navier-Stokes equations. The first approximation, FDA1, givenby Eqs. (39), was constructed in the paper [2] where its s-consistency was established. Thedifference approximation FDA4, given by Eqs. (42) with the last term omitted, D · ˜∆ ˜ u ,which is reduced to zero modulo ˜ F (4)0 (cf. Remark 4), was derived in [19] by the methodsuggested in [20]. Remark 5.
The equation ˜ F (2)3 = 0 in FDA2 provides a compact finite difference discretizationof the Poisson pressure equation (27). In this case( D ˜ u ) + 2 ( D ˜ v )( D ˜ u ) + ( D ˜ v ) −−−→ h → u x + 2 v x u y + v y (43)whereas ∇ · ( u · ∇ ) u = u x + 2 v x u y + v y + uu xx + vv yy + vu xy + uv xy , (44)and the right-hand sides of Eqs. (43) and (44) are equal modulo the continuity equation u x + v y = 0.However, FDA2 is s-inconsistent, since ˜ F (2)3 (cid:54)∈ (cid:74) ˜ I (cid:75) , where (cid:74) ˜ I (cid:75) ⊂ ˜ R the ideal generated by { ˜ F (2)0 , ˜ F (2)1 , ˜ F (2)2 } by virtue of inequality D · D = σ + σ − σ − + σ − h (cid:54) = ˜∆ , and there are consequences of Eqs. (40) implying the differential equations (cf. Definition 18)which are not consequences of (25). Below we explicitly demonstrate this for the linearizedversion of Eqs. (23) called Stokes equations.The approximation FDA3 given by Eqs. (41) differs from FDA1 in the structure of discretePoisson pressure equation. In contrast to equation ˜ F (3)3 = 0 in Eqs. (41), the equation ˜ F (1)3 = 0in (39) contains extra part D · ˜∆ ˜ u and can be omitted if ˜ u satisfies the discrete continuityequation (see Remark 4).In the difference approximation FDA4, Eqs. (42), the discrete pressure Poisson equation˜ F (4)3 = 0, as opposed to that in FDA2, provides s-consistency of FDA4.We compare these four schemes by using the following absolute/relative error formula e ng = max j,k | g nj,k − g ( x j , y k , t n ) | | g ( x j , y k , t n ) | , (45)where g ∈ { u, v, p } and g ( x, y, t ) belongs to the exact solution (36).The governing differential system (25), (27) is mixed elliptic-parabolic: the momentumequations are parabolic and the pressure Poisson equation is elliptic. It should be particularlyemphasized that in our construction of a numerical solution for the initial-value problemwith initial data taken from Eqs. (36) at t = 0, we use the discrete momentum equationsto determine velocities and the discrete pressure Poisson equation to determine pressure. Inother words, we use the classical pressure-Poisson formulation of the Navier-Stokes equations(cf. [48], Sect. 3.2) to solve the initial-value problem numerically. However, as this takesplace, to construct numerical solution of the above given difference approximations we do notexploit the discrete divergence-free constraint D · ˜ u = 0 and use (cf. [1, 2]) and use it only forverification of the obtained results.We compute numerical solutions in the domain [0 , π ] × [0 , π ] × [0 ,
6] with the Reynoldsnumber Re = 100. Figures 1–3 contain the computed error for three different choices of h (error in u and v coincides). We let τ = 0 . × h .19he results of our computational experiments shown in Figures 1 and 2. Except FDA3,which is unstable (Fig. 2, top) by the lack of mass conservation (violation of the incompressibil-ity condition), the other approximations clearly demonstrate the second order of convergencewith respect to h , during which FDA1 far exceeds the others in accuracy. In addition, basedon the numerical velocities obtained for the last FDA, the continuity equation is accurateto 10 − (for h = 0 . || D · ˜u || F := (cid:16) (cid:88) j,k | D ˜ u nj,k + D ˜ v nj,k | (cid:17) . (46)The superiority of FDA1 over the others FDAs is due to incorporation of s-consistency,conservativity (divergence form) of the nonlinear terms in momentum equations and presenceof the last term in the pressure Poisson equation (32). In its turn, FDA2 provides loweraccuracy in comparison with FDA4 because of its s-inconsistency. On the other hand, FDA2the more stable than FDA3, since the former unlike the latter was constructed with applicationof the incompressibility condition (see Remark 5). If one correlates FDA1 with FDA4, thenthe last one does not have a conservation law (divergence) form and by this reason its accuracyis not so good.Nearly all known finite difference approaches to solving incompressible Navier-Stokes equa-tions in terms of (’primitive’) variables { ˜ u , ˜ p } started from the (fractional step) projectionmethod based on presentation of the vector momentum equation in (25) as a Helmholtz-Hodge decomposition [9, 27, 48] and on the use of the continuity equation (incompressibilitycondition) for correction of the velocity vector ˜ u on every time step. Our numerical experi-ments with the scheme FDA1 show, contrastingly, that it is sufficient to attain the fulfilmentof the continuity equation by initial and/or boundary conditions.To illustrate the fact that s-inconsistency has an adverse effect on the solution space of thediscretized equations, we consider the incompressible Stokes equations flow given in vectornotation by ∂ t u + ∇ p − u = 0 , ∇ · u = 0 . (47)The linear PDE system (47) approximates (25) when Re is small (cf. [40], Ch. 22 · ρ g ( ρ is the liquid density, g is the acceleration due to gravity)in the right-hand side of the first equation, Eqs. (47) have numerous applications in the de-scription of fluid displacement processes that take place in porous media (see [15] and thereferences therein) and that are related to both chemical and physical phenomena.For Eqs. (47) the pressure Poisson equation (27) becomes the pressure Laplace equation ∆ p = 0 , (48)and let the PDE system (47)–(48) be discretized as follows D t ˜u + D ˜ p − ˜u = 0 , D · ˜u = 0 , ˜∆ ˜ p = 0 , (49)where the difference approximations of partial derivatives given by Eqs. (30).The passive form of Eqs. (49) reads D t ˜u + D ˜ p − ˜u = 0 , D · ˜u = 0 , ( D · D ) ˜ p = 0 , ˜∆ ˜ p = 0 , ˜ q = 0 . t FDA1 (u, v), h = 0.100 p, h = 0.100(u, v), h = 0.050 p, h = 0.050(u, v), h = 0.025 p, h = 0.0250 1 2 3 4 5 6 t FDA2 (u, v), h = 0.100 p, h = 0.100(u, v), h = 0.050 p, h = 0.050(u, v), h = 0.025 p, h = 0.025
Figure 1: Taylor decaying problem: error with different h in the computed solution withFDA1 scheme and second-order standard discretizations FDA221 t FDA3 (u, v), h = 0.100 p, h = 0.100(u, v), h = 0.050 p, h = 0.050(u, v), h = 0.025 p, h = 0.0250 1 2 3 4 5 6 t FDA4 (u, v), h = 0.100 p, h = 0.100(u, v), h = 0.050 p, h = 0.050(u, v), h = 0.025 p, h = 0.025
Figure 2: Taylor decaying problem: error with different h in the computed solution withFDA3 and FDA4 scheme 22 t || D u || h = 0.025 FDA1FDA2FDA3FDA4
Figure 3: Taylor decaying problem: computed value of error (46) for all four schemesHere the difference polynomial ˜ q is the reduced S -polynomial of the two preceding equationsand, in accordance with Definition 18, it implies0˜ q := S (cid:16) ( D · D ) ˜ p, ˜∆ ˜ p (cid:17) (cid:66) ∂ i ∂ j p , i, j ∈ { , , } . (50)It is clear that FDA (49) is w-consistent with Eqs. (47) and (48). However, none of thedifferential equations occurring in (50) is a consequence of Eqs. (47). This can be explicitlyverified with the Maple packages LDA [24] and
Janet [8] by computing the related normalforms (cf. (22)). Therefore, FDA (49) is s-inconsistent with the PDE system (47), (48) inaccordance with Theorem 1.The difference polynomials in Eqs. (49) generate the perfect difference ideal whose element˜ q , as indicated in Eq. (50), in the continuous limit implies additional equations for the pressure.These equations together with the Laplace equation (48) restrict the pressure component tothe exact solution p := c + c x + c x + c x + c x + c x x + c x + c x x + c x + c x x + c x x x , where c + c + c = 0 , where c i (0 ≤ i ≤
10) are arbitrary functions of t satisfying the above constraint. It isclear that, in comparison with the s-consistent approximations, the s-inconsistent one (49)significantly decreases the domain of solutions to the governing differential system that canbe successfully constructed by numerical methods.23 Difference Thomas decomposition and s-consistency check
Let ˜ S be a system of polynomial partial difference equations and inequations ˜ f = 0 , . . . , ˜ f s = 0 , ˜ f s +1 (cid:54) = 0 , . . . , ˜ f s + t (cid:54) = 0 ( s, t ∈ Z ≥ ) . (51)Here ˜ f , . . . , ˜ f s + t are elements of the difference polynomial ring ˜ R in ˜ u (1) , . . . , ˜ u ( m ) withcommuting automorphisms Σ = { σ , . . . , σ n } . A ranking (cid:31) on ˜ R is fixed, so that leaders,initials and discriminants of non-constant difference polynomials are defined as in Definition 1.In this section we develop a difference analogue of the Thomas decomposition method fordifferential systems (cf. Section 2). The resulting main Algorithm 4 is based, in particular,on Algorithm 2 for auto-reduction of a difference system, which precedes the assignment ofadmissible automorphisms by Janet division, and on Algorithm 3 computing Janet normalforms modulo a Janet complete difference system, in order to check passivity.Algorithm 2 performs reductions on a finite system of difference equations (given by ˜ L ),if possible, removing zero remainders from the system. If a reduction occurs that results in anon-zero remainder, the original polynomial is replaced by this remainder and the algorithmstops. Since in that case further splitting of the system may be necessary to ensure non-va-nishing of initials, this situation is indicated by a flag returned to the main Algorithm 4.The following notation is used in what follows. For a difference system ˜ S as above let ˜ S = (resp. ˜ S (cid:54) = ) be the set { ˜ f , . . . , ˜ f s } (resp. { ˜ f s +1 , . . . , ˜ f s + t } ). Let ˜ E be a difference ideal of ˜ R and let ∅ (cid:54) = ˜ Q ⊆ ˜ R be multiplicatively closed and closed under σ , . . . , σ n . Then define˜ E : ˜ Q := { ˜ f ∈ ˜ R | ˜ q ˜ f ∈ ˜ E for some ˜ q ∈ ˜ Q } . Moreover, for U ⊆ Mon(Σ) ˜ u and v ∈ Mon(Σ) ˜ u we define U : v := { θ ∈ Mon(Σ) | θ v ∈ U } . Algorithm 2:
Auto-reduce for difference algebra
Input: ˜ L ⊂ ˜ R \ ˜ K finite and a ranking (cid:31) on ˜ R such that ˜ L = ˜ S = for some finitedifference system ˜ S which is quasi-simple as an algebraic system (in the finitely manyindeterminates ˜ u ( α ) J occurring in it, totally ordered by (cid:31) ) Output: a ∈ { true , false } and ˜ L (cid:48) ⊂ ˜ R \ ˜ K finite such that[ ˜ L (cid:48) ] : ˜ Q = [ ˜ L ] : ˜ Q , where ˜ Q is the smallest multiplicatively closed subset of ˜ R containing all init( θ ˜ f ),where ˜ f ∈ ˜ L and θ ∈ ld( ˜ L \ { ˜ f } ) : ld( ˜ f ), and which is closed under σ , . . . , σ n , and, incase a = true , there exist no ˜ f , ˜ f ∈ ˜ L (cid:48) , ˜ f (cid:54) = ˜ f , such that we have v := ld( ˜ f ) = θ ld( ˜ f ) for some θ ∈ Mon(Σ) and deg v ( ˜ f ) ≥ deg v ( θ ˜ f ) ˜ L (cid:48) ← ˜ L while ∃ ˜ f , ˜ f ∈ ˜ L (cid:48) , ˜ f (cid:54) = ˜ f and θ ∈ Mon(Σ) such that we have v := ld( ˜ f ) = θ ld( ˜ f ) and deg v ( ˜ f ) ≥ deg v ( θ ˜ f ) do ˜ L (cid:48) ← ˜ L (cid:48) \ { ˜ f } ; v ← ld( ˜ f ) ˜ r ← init( θ ˜ f ) ˜ f − init( ˜ f ) v d θ ˜ f , d := deg v ( ˜ f ) − deg v ( θ ˜ f ) if ˜ r (cid:54) = 0 then return ( false , ˜ L (cid:48) ∪ { ˜ r } ) return ( true , ˜ L (cid:48) ) 24ince leaders are dealt with in decreasing order with respect to (cid:31) , and no ranking admitsinfinitely decreasing chains (cf. [33, Ch. 0, Sect. 17, Lemma 15]), Algorithm 2 terminates . Its correctness follows from the definition of ˜ E : ˜ Q .Before presenting Algorithm 3 which computes the Janet normal form of a differencepolynomial, we adapt the discussion of Janet division preceding Definition 4 to the differencecase.Janet division associates (with respect to a total ordering of Σ) to each ˜ f i = 0 withld( ˜ f i ) = θ i ˜ u ( α ) the set µ i := µ ( θ i , ˜ G α ) ⊆ Σ (resp. µ i := Σ \ µ i ) of admissible (resp. non-admissible ) automorphisms , where˜ G α := { θ ∈ Mon(Σ) | θ ˜ u ( α ) ∈ { ld( ˜ f ) , . . . , ld( ˜ f s ) } } . We call { ˜ f = 0 , . . . , ˜ f s = 0 } or T := { ( ˜ f , µ ) , . . . , ( ˜ f s , µ s ) } Janet complete if each ˜ G α equalsits Janet completion, α = 1, . . . , m .Let ˜ r ∈ ˜ R . If some v ∈ Mon(Σ)˜ u occurs in ˜ r for which there exists ( ˜ f , µ ) ∈ T such that v = θ ld( ˜ f ) for some θ ∈ Mon( µ ) and deg v (˜ r ) ≥ deg v ( θ ˜ f ), then ˜ r is Janet reducible modulo T . In this case, ( ˜ f , µ ) is called a Janet divisor of ˜ r . If ˜ r is not Janet reducible modulo T ,then ˜ r is also said to be Janet reduced modulo T . Iterated pseudo-reductions of ˜ r modulo T yield its Janet normal form
NF(˜ r, T, (cid:31) ), which is the Janet reduced difference polynomial ˜ r (cid:48) returned by Algorithm 3. Definition 21.
Let T = { ( ˜ f , µ ) , . . . , ( ˜ f s , µ s ) } be Janet complete. The difference system { ˜ f = 0 , . . . , ˜ f s = 0 } or T is said to be passive , if the following Janet passivity conditionshold: NF( σ ˜ f i , T, (cid:31) ) = 0 for all σ ∈ µ i = Σ \ µ i , i = 1 , . . . , s . (52)Note that Eqs. (52), in general, form a proper subset of the Gr¨obner passivity condi-tions (13) (cf. [17]). Definition 22.
Let (cid:31) be a ranking on ˜ R , and fix a total ordering on Σ with respect towhich Janet division is defined. A difference system ˜ S as in (51) is said to be simple (resp., quasi-simple ) if the following three conditions are satisfied.1. ˜ S is simple (resp., quasi-simple) as an algebraic system (in the finitely many indeter-minates which occur in the equations and inequations of ˜ S , totally ordered by (cid:31) ; cf.Definitions 2 and 3).2. The difference system { ˜ f = 0 , . . . , ˜ f s = 0 } is passive.3. The left hand sides ˜ f s +1 , . . . , ˜ f s + t are Janet reduced modulo the passive differencesystem { ˜ f = 0 , . . . , ˜ f s = 0 } . Theorem 2.
Let ˜ S be a quasi-simple difference system over ˜ R as in (51). Let ˜ E be thedifference ideal of ˜ R generated by ˜ f , . . . , ˜ f s and let ˜ Q be the smallest subset of ˜ R which ismultiplicatively closed, closed under σ , . . . , σ n and contains the initials ˜ q i := init( ˜ f i ) for all i = 1 , . . . , s . Then a difference polynomial ˜ f ∈ ˜ R is an element of ˜ E : ˜ Q = { ˜ f ∈ ˜ R | ( θ (˜ q )) r . . . ( θ s (˜ q s )) r s ˜ f ∈ ˜ E for some θ , . . . , θ s ∈ Mon(Σ) , r , . . . , r s ∈ Z ≥ } if and only if the Janet normal form of ˜ f modulo { ˜ f , . . . , ˜ f s } is zero. lgorithm 3: Janet-reduce for difference algebra
Input: ˜ r ∈ ˜ R , T = { ( ˜ f , µ ) , ( ˜ f , µ ) , . . . , ( ˜ f s , µ s ) } , and a ranking (cid:31) on ˜ R , where T isJanet complete (with respect to (cid:31) ) Output: (˜ r (cid:48) , ˜ b ) ∈ ˜ R × ˜ R such that (1) if ˜ r ∈ ˜ K or T = ∅ , then ˜ r (cid:48) = ˜ r , ˜ b = 1, (2)otherwise ˜ r (cid:48) is Janet reduced modulo T and˜ r (cid:48) + [ ˜ f , . . . , ˜ f s ] = ˜ b · ˜ r + [ ˜ f , . . . , ˜ f s ] , where ˜ b is in the multiplicatively closed set generated by s (cid:91) i =1 { θ init( ˜ f i ) | θ ∈ Mon(Σ) , ld(˜ r ) (cid:31) θ ld( ˜ f i ) } ∪ { } ˜ r (cid:48) ← ˜ r ; ˜ b ← if ˜ r (cid:48) (cid:54)∈ ˜ K then v ← ld(˜ r (cid:48) ) while ˜ r (cid:48) (cid:54)∈ ˜ K and there exist ( ˜ f , µ ) ∈ T and θ ∈ Mon( µ ) such that v = θ ld( ˜ f ) and deg v (˜ r (cid:48) ) ≥ deg v ( θ ˜ f ) do ˜ r (cid:48) ← init( θ ˜ f ) ˜ r (cid:48) − init(˜ r (cid:48) ) v d θ ˜ f , d := deg v (˜ r (cid:48) ) − deg v ( θ ˜ f ) ˜ b ← init( θ ˜ f ) · ˜ b for each coefficient ˜ c of ˜ r (cid:48) (as a polynomial in v ) do (˜ r (cid:48)(cid:48) , ˜ b (cid:48) ) ← Janet-reduce (˜ c , T , (cid:31) ) replace the coefficient ˜ b (cid:48) · ˜ c in ˜ b (cid:48) · ˜ r (cid:48) with ˜ r (cid:48)(cid:48) and replace ˜ r (cid:48) with this result ˜ b ← ˜ b (cid:48) · ˜ b return (˜ r (cid:48) , ˜ b ) 26 roof. By definition of ˜ E : ˜ Q , every element ˜ f ∈ ˜ R for which Algorithm 3 yields Janet normalform zero is an element of ˜ E : ˜ Q .Let ˜ f ∈ ˜ E : ˜ Q , ˜ f (cid:54) = 0. Then there exist ˜ q ∈ ˜ Q and k , . . . , k s ∈ Z ≥ and ˜ c i,j ∈ ˜ R \ { } , α i,j ∈ Mon(Σ), j = 1, . . . , k i , i = 1, . . . , s , such that˜ q ˜ f = s (cid:88) i =1 k i (cid:88) j =1 ˜ c i,j α i,j ( ˜ f i ) . (53)Among all pairs ( i, j ) for which α i,j involves a non-admissible automorphism for ˜ f i = 0 letthe pair ( i (cid:63) , j (cid:63) ) be such that α i (cid:63) ,j (cid:63) (ld( ˜ f i (cid:63) )) is maximal with respect to the ranking (cid:31) . Let σ be a non-admissible automorphism for ˜ f i (cid:63) = 0 which divides the monomial α i (cid:63) ,j (cid:63) . Since { ˜ f = 0 , . . . , ˜ f s = 0 } is a passive difference system, there exist ˜ b ∈ ˜ Q and l , . . . , l s ∈ Z ≥ and˜ d i,j ∈ ˜ R \ { } as well as β i,j ∈ Mon(Σ), j = 1, . . . , l i , i = 1, . . . , s , such that˜ b · ( σ ˜ f i (cid:63) ) = s (cid:88) i =1 l i (cid:88) j =1 ˜ d i,j β i,j ( ˜ f i ) , where each β i,j involves only admissible automorphisms for ˜ f i = 0. Let γ i (cid:63) ,j (cid:63) := α i (cid:63) ,j (cid:63) /σ andmultiply (53) by γ i (cid:63) ,j (cid:63) (˜ b ) to obtain γ i (cid:63) ,j (cid:63) (˜ b ) · ˜ q ˜ f = s (cid:88) i =1 k i (cid:88) j =1 ˜ c i,j · γ i (cid:63) ,j (cid:63) (˜ b ) · α i,j ( ˜ f i ) . In this equation we replace γ i (cid:63) ,j (cid:63) (˜ b ) · α i (cid:63) ,j (cid:63) ( ˜ f i (cid:63) ) = γ i (cid:63) ,j (cid:63) (˜ b · σ ( ˜ f i (cid:63) ))by γ i (cid:63) ,j (cid:63) s (cid:88) i =1 l i (cid:88) j =1 ˜ d i,j β i,j ( ˜ f i ) . Since γ i (cid:63) ,j (cid:63) β i (cid:63) ,j (cid:63) involves fewer non-admissible automorphisms for ˜ f i = 0 than α i (cid:63) ,j (cid:63) , iterationof this substitution process will rewrite equation (53) in such a way that every α i,j (ld( ˜ f i ))involving non-admissible automorphisms for ˜ f i = 0 will be less than α i (cid:63) ,j (cid:63) (ld( ˜ f i (cid:63) )) with respectto (cid:31) . A further iteration of this substitution process will therefore produce an equation as(53) with no α i,j involving any non-admissible automorphisms for ˜ f i = 0.This shows that for every ˜ f ∈ ( ˜ E : ˜ Q ) \ { } there exists a Janet divisor of ld( ˜ f ) in thepassive set defined by ˜ f = 0, . . . , ˜ f s = 0. Corollary 1.
In the situation of Theorem 2 let ˜ S be simple. Then the difference ideal ˜ E : ˜ Q is radical.Proof. Let ˜ f ∈ ˜ R and r ∈ N be such that ˜ f r ∈ ˜ E : ˜ Q . We will show that ˜ f ∈ ˜ E : ˜ Q . Since˜ f r ∈ ˜ E : ˜ Q and the difference system ˜ S is (quasi-) simple, there exist ˜ q ∈ ˜ Q , k , . . . , k s ∈ Z ≥ ,˜ c i,j ∈ ˜ R , α i,j ∈ Mon(Σ), j = 1, . . . , k i , i = 1, . . . , s , such that˜ q ˜ f r = s (cid:88) i =1 k i (cid:88) j =1 ˜ c i,j α i,j ( ˜ f i ) , (54)27here each α i,j only involves admissible automorphisms for ˜ f i = 0 and where ˜ q is a productof powers of init( α i,j ( ˜ f i )), i = 1, . . . , s , j = 1, . . . , k i .Let V ⊂ Mon(Σ)˜ u be minimal such that the (non-difference) polynomial algebra ˜ K [ V ] ⊂ ˜ R contains all indeterminates occurring in (54). Note that V is finite and recall that ˜ S is simpleas an algebraic system (cf. Definition 2). Now define the algebraic system (over ˜ K [ V ] with V totally ordered by (cid:31) )˜ S (cid:48) = { α i,j ( ˜ f i ) = 0 | i = 1 , . . . , s, j = 1 , . . . , k i } ∪ { ˜ f s +1 (cid:54) = 0 , . . . , ˜ f s + t (cid:54) = 0 } . Then ˜ S (cid:48) is simple. In fact, the leaders of all equations and inequations in ˜ S (cid:48) are pairwisedistinct, because the cones of monomials in σ , . . . , σ n defined by applying admissible auto-morphisms to the leaders of ˜ f , . . . , ˜ f s are disjoint (cf. the discussion before Definition 4),and vanishing of init( α i,j ( ˜ f i )) = α i,j (init( ˜ f i )) or disc( α i,j ( ˜ f i )) = α i,j (disc( ˜ f i )) on the solutionset of the algebraic system ˜ S (cid:48) is prevented by the simplicity of ˜ S .Let I V be the (algebraic) ideal of ˜ K [ V ] that is generated by the equations in ˜ S (cid:48) , andlet ˜ q (cid:48) be the product of the initials of all equations in ˜ S (cid:48) . Then equation (54) shows that˜ f r ∈ I V : (˜ q (cid:48) ) ∞ . Since the algebraic system ˜ S (cid:48) is simple, [51, Prop. 2.2.7] shows that the ideal I V : (˜ q (cid:48) ) ∞ is radical. Hence, ˜ f ∈ I V : (˜ q (cid:48) ) ∞ ⊂ ˜ E : ˜ Q .Let Ω ⊆ R n be open and connected and fix z ∈ Ω. Denoting the grid in (6) by Γ z , h , wedefine F Ω , z , h := { ˜ u : Γ z , h ∩ Ω → C | ˜ u is the restriction to Γ z , h ∩ Ω ofsome locally analytic function u on Ω } , and for a system ˜ S of partial difference equations and inequations as in (51) we define thesolution set Sol Ω , z , h ( ˜ S ) := { ˜ u ∈ F Ω , z , h | ˜ f i (˜ u ) = 0 , ˜ f s + j (˜ u ) (cid:54) = 0 for all i = 1 , . . . , s, j = 1 , . . . , t } . Definition 23.
Let ˜ S be a finite difference system over ˜ R and (cid:31) a ranking on ˜ R . A differencedecomposition of ˜ S (with respect to (cid:31) ) is a finite collection of quasi-simple difference systems˜ S , . . . , ˜ S r over ˜ R such that Sol Ω , z , h ( ˜ S ) = Sol Ω , z , h ( ˜ S ) (cid:93) . . . (cid:93) Sol Ω , z , h ( ˜ S r ).Given a finite difference system ˜ S over ˜ R , Algorithm 4, presented below, constructs adifference decomposition of ˜ S in finitely many steps. In step 11 Decompose refers to analgorithm which computes a smallest superset of ˜ G = { ˜ f , . . . , ˜ f s } in ˜ R that is Janet completeas defined on page 24 (cf., for example, [51, Algorithm 2.1.6]). Theorem 3.
Algorithm 4 terminates and is correct.Proof.
Algorithm 4 maintains a set Q of difference systems that still have to be dealt with.Given that termination of all subalgorithms has been proved, termination of Algorithm 4 isequivalent to the condition that Q = ∅ holds after finitely many steps.Apart from step 1, new systems are inserted into Q in steps 18 and 20. We consider thesystems that are at some point an element of Q as the vertices of a tree. The root of thistree is the input system ˜ S . The systems which are inserted into Q in steps 18 and 20 are thevertices of the tree whose ancestor is the system ˜ L that was extracted from Q in step 3 whichin the following steps produced these new systems. Since the for loop beginning in step 5terminates, the degree of each vertex in the tree is finite. We claim that every branch of thetree is finite, i.e., that the tree has finite height, hence, that the tree has only finitely manyvertices. 28 lgorithm 4: DifferenceDecomposition
Input:
A finite difference system ˜ S over ˜ R , a ranking (cid:31) on ˜ R , and a total ordering onΣ (used by Decompose ) Output:
A difference decomposition of ˜ S Q ← { ˜ S } ; T ← ∅ repeat choose ˜ L ∈ Q and remove ˜ L from Q compute a decomposition { A , . . . , A r } of ˜ L , considered as an algebraic system,into quasi-simple systems (cf. Definition 3) for i = 1 , . . . , r do if A i = ∅ then // no equation and no inequation return {∅} else ( a, ˜ G ) ← Auto-reduce ( A = i , (cid:31) ) // Algorithm 2 if a = true then J ← Decompose ( ˜ G ) P ← { NF( σ ˜ f , J, (cid:31) ) | ( ˜ f , µ ) ∈ J, σ ∈ µ } // Algorithm 3 if P ⊆ { } then // J is passive replace each ˜ g (cid:54) = 0 in A i with NF(˜ g, J, (cid:31) ) (cid:54) = 0 if (cid:54)∈ A (cid:54) = i then insert { ˜ f = 0 | ( ˜ f , µ ) ∈ J } ∪ { ˜ g (cid:54) = 0 | ˜ g ∈ A (cid:54) = i } into T else if P ∩ ˜ K ⊆ { } then insert { ˜ f = 0 | ( ˜ f , µ ) ∈ J } ∪ { ˜ f = 0 | ˜ f ∈ P \ { }} ∪ { ˜ g (cid:54) = 0 | ˜ g ∈ A (cid:54) = i } into Q else insert { ˜ f = 0 | ˜ f ∈ G }∪{ ˜ g (cid:54) = 0 | ˜ g ∈ A (cid:54) = i } into Q until Q = ∅ return T
29n case of step 20 the new system contains an equation which resulted from a non-trivialdifference reduction in step 9. When this new system will be extracted from Q in a laterround, a decomposition into quasi-simple algebraic systems will be computed in step 4. Thismay produce new branches of the tree, but along any of these branches, after finitely manysteps the condition a = true in step 10 will hold, because the order of the shifts in leaders ofthe arising equations is bounded by the maximum order of shifts in leaders of the ancestorsystem ˜ L .In case of step 18 we are going to show that after finitely many steps a difference equationis obtained whose leader has not shown up as a leader of an equation in any preceding systemin the current branch of the tree. First of all, the passivity check (step 12) yielded an equation˜ f = 0, ˜ f ∈ P \ ˜ K , which is Janet reduced modulo J . Hence, either ld( ˜ f ) is not containedin the multiple-closed set generated by ld( ˜ G ), or there exists ( ˜ f (cid:48) , µ (cid:48) ) ∈ J such that ld( ˜ f (cid:48) )is a Janet divisor of ld( ˜ f ), but the degree of ˜ f in ld( ˜ f ) is smaller than the degree of ˜ f (cid:48) inld( ˜ f (cid:48) ). In the first case the above claim holds. The second case cannot repeat indefinitely:First of all, if ld( ˜ f ) = ld( ˜ f (cid:48) ), then in a later round, either a pseudo-reduction of ˜ f (cid:48) modulo˜ f will be performed if the initial of ˜ f does not vanish, or init( ˜ f ) = 0 has been added as anew equation (with lower ranked leader). Since this leads to a sequence in Mon(Σ) whichstrictly decreases, infinite chains are excluded in this situation. If case ld( ˜ f ) (cid:54) = ld( ˜ f (cid:48) ) occursrepeatedly, then a sequence (( θ i ˜ u ( α ) ) e i ) i =1 , , ,... of leaders of newly inserted equations arises,where θ i ∈ Mon(Σ), α ∈ { , . . . , m } , e i ∈ Z ≥ , such that e i +1 < e i holds (and where also θ i | θ i +1 ). Any such sequence is finite. Hence, the first case arises after finitely many steps.Therefore, termination follows from Dickson’s Lemma.In order to prove correctness, we note that a difference system is only inserted into T ifstep 12 confirmed passivity. Such a system is quasi-simple as an algebraic system because (upto auto-reduction in step 9 and Janet completion in step 11) it was returned as one system A i in step 4. Condition (3) in Definition 22 is ensured by step 14. Hence, all differencesystems in T are quasi-simple. Splitting of a system only arises in step 4 by adding anequation init( ˜ f ) = 0 and the corresponding inequation init( ˜ f ) (cid:54) = 0, respectively, to the twonew systems replacing the given one. Since no solutions are lost or gained, this leads to apartition as required by Definition 23.Given a simple differential system and its w-consistent discretization on the regular grid (6),Algorithm 5 allows to verify strong consistency of the latter. Correctness of the algorithm follows from Definition 19 (extended to inequations), Def-inition 20 and passivity of the subsystems returned by Algorithm 4. Their solution spacespartition the solution space of the input FDA. Thereby, any subsystem ˜ L i in the output with b i = true is s-consistent with L i , where ˜ L i −−−−→ | h |→ L i and w-consistent if b i = false . If b i = true for all i , then ˜ S is s-consistent with S . Termination follows from that of the subalgorithms.
In this section we consider two systems of quasi-linear PDEs for unknown functions of twoindependent variables. 30 lgorithm 5:
S-ConsistencyCheck
Input:
A simple differential system S over R , a differential ranking (cid:31) on R , adifference ranking > on ˜ R , a total ordering on Σ (used by Decompose ) and adifference system ˜ S consisting of equations that are w-consistent with S Output: ˜ L = { ( ˜ L , b ) , . . . , ( ˜ L r , b r ) } , where ˜ L i is s-consistent (resp. w-consistent) with L i ←−−−− | h |→ ˜ L i if b i = true (resp. false ) ˜ L = { ˜ L , . . . , ˜ L k } ← DifferenceDecomposition ( ˜
S, > ) for i = 1 , . . . , k do if ∃ ˜ f ∈ ˜ L (cid:54) = i such that ˜ f (cid:66) F with F ∩ (cid:74) S = (cid:75) (cid:54) = ∅ then // Definition 18 ˜ L ← ˜ L \ { ˜ L i } else b i ← true for ˜ f ∈ ˜ L = do compute F ⊂ R such that ˜ f (cid:66) F // Remark 2 if ∃ f ∈ F such that NF( f, S = , (cid:31) ) (cid:54) = 0 then b i ← false ; break return { ( ˜ L i , b i ) | ˜ L i ∈ ˜ L } Example 5.
Let us consider the overdetermined PDE system (cid:40) u x − u = 0 ,u y + u = 0 , u = u ( x, y ) . (55)Since the cross-derivative ∂ y ( u x − u ) − ∂ x ( u y + u ) reduces to zero modulo the given equations,the differential system (55) is simple (with respect to any ranking). The exact general solutionof (55) can easily be found with Maple: u ( x, y ) = 1 y − x + C , (56)where C is an arbitrary constant (the corresponding counting polynomial [35] being ∞ ). Forthe numerical comparison of the following finite difference approximations we consider thedomain [0 , × [0 ,
10] with Cartesian grid defined by h = h = h = 1 / C = 12. The error will be computed as( e g ) j,k := | g j,k − g ( x j , y k ) | | g ( x j , y k ) | , maximum error := max j,k ( e g ) j,k , where g is the exact solution (56).We investigate the system of difference equations ˜ u i +1 ,j − ˜ u i,j h − ˜ u i,j = 0 , ˜ u i,j +1 − ˜ u i,j h + ˜ u i,j = 0 , (57)which is obtained as discretization of (55) by replacing u x and u y by the corresponding forwarddifferences, with step sizes h and h , respectively. For simplicity we ignore case distinctions31 y Figure 4: Computed error for discretization (57), maximum error = 0.046222672401600905and pursue the generic case only. We fix an orderly ranking on the difference polynomial ring˜ R = ˜ K{ ˜ u } with automorphisms σ , σ .Denote by ˜ f and ˜ f the left hand sides in (57). Then (57) is simple as an algebraicsystem, but the passivity check (cf. Definition 17) reveals the consequence h σ ˜ f − h σ ˜ f − (cid:0) h ˜ u i +1 ,j + h h ˜ u i,j + h ˜ u i,j − (cid:1) h ˜ f − (cid:0) h ˜ u i,j +1 − h h ˜ u i,j + h ˜ u i,j + 1 (cid:1) h ˜ f = h h ( h + h ) ˜ u i,j . Note that h h ( h + h ) is non-zero. By adding ˜ f := ˜ u i,j to system (57) we obtain thequasi-simple difference system { ˜ f = 0 , ˜ f = 0 , ˜ f = 0 } . The continuous limit of ˜ f for h → h → u , which is not in the radical differential idealcorresponding to (55). Hence, (57) is not s-consistent with (55).Figure 4 shows the error computed for FDA (57) relative to the exact solution (56).Hereafter, in our numerical computation we chose the grid spacings as h = h = h = 0 . x = 10 , y = 0), which is closest tothe pole in (56). Because of s-inconsistency of this discretization we did not compute theassociated modified PDE (cf. [62]). In general terms a modified PDE for a given FDA is onethat a numerical solution to FDA satisfies to a higher accuracy than the initial PDE (see, forexample, the textbooks [41], Sections 5.5–5.6 and [58], Section 7.7). The method of modifiedequation provides a useful tool for evaluating such important properties of finite differenceschemes as order of accuracy, consistency, stability, dissipation and dispersion.Next we consider the discretization obtained by replacing u x by the forward difference asbefore and u y by the backward difference: ˜ u i +1 ,j − ˜ u i,j h − ˜ u i,j = 0 , ˜ u i,j +1 − ˜ u i,j h + ˜ u i,j +1 = 0 , (58)32 y Figure 5: Computed error for discretization (58), maximum error = 0.038621319610305495again with step sizes h and h , respectively. Denote by ˜ f (cid:48) and ˜ f (cid:48) the left hand sides in (58).The passivity check reveals the consequence (with underlined leader) h (cid:16) σ ˜ f (cid:48) − h ˜ u i +1 ,j +1 σ ˜ f (cid:48) (cid:17) − (cid:0) h h ˜ u i,j +1 + h ˜ u i,j +1 + 1 (cid:1) h h σ ˜ f (cid:48) − (cid:16) h h ˜ u i,j +1 − h h ( h − h ) ˜ u i,j +1 + h h ˜ u i,j + h − h h + h (cid:17) h ˜ f (cid:48) = − h ( h − h ) (cid:16) (2 h ˜ u i,j + 1) ˜ u i,j +1 − h ˜ u i,j − ˜ u i,j (cid:17) . The continuous limit of this difference polynomial ˜ f (cid:48) is given by (cid:0) u y + u (cid:1) h h − (cid:0) u y + u (cid:1) h h . Now a pseudo-reduction of ˜ f (cid:48) modulo ˜ f (cid:48) yields h h ( h − h ) (2 h ˜ u i,j + 1) ˜ f (cid:48) + (cid:0) h (2 h ˜ u i,j + 1) ˜ u i,j +1 + h ˜ u i,j + 3 h ˜ u i,j + 1 (cid:1) ˜ f (cid:48) = h h ( h − h ) ˜ u i,j . For h (cid:54) = h we define ˜ f (cid:48) := ˜ u i,j and obtain the quasi-simple difference system { ˜ f (cid:48) = 0 , ˜ f (cid:48) =0 , ˜ f (cid:48) = 0 } , and we conclude that (58) is not s-consistent with (55). However, if h = h , then(58) is a simple difference system and it is s-consistent with (55).Now we perform the Taylor expansion of the left-hand sides in Eqs. (58) with h = h = h and explicitly write the first-order terms f := u x − u + u x,x h + O ( h ) = 0 ,f := u y + u + (cid:16) uu y + u y,y (cid:17) h + O ( h ) = 0 . (59)From Eqs. (59) we obtain the modified PDE for scheme (58) f − (cid:18)
12 ( f ) x + uf (cid:19) h = u x − u + u h + O ( h ) = 0 ,f − (cid:18)
12 ( f ) y + uf (cid:19) h = u y + u − u h + O ( h ) = 0 , (60)33 y Figure 6: Computed error for discretization (61), maximum error = 0.002446111801807538which shows that scheme (58) has first order accuracy. Furthermore, Eqs. (60) allow anobvious modification (cf. [52], p. 80) of FDA (58) to one with second order accuracy given by ˜ u i +1 ,j − ˜ u i,j h − ˜ u i,j − h ˜ u i,j = 0 , ˜ u i,j +1 − ˜ u i,j h + ˜ u i,j +1 + h ˜ u i,j = 0 . (61)The corresponding decrease of numerical error for FDA (61) in comparison with FDA (58)(Fig. 5) is shown in Fig. 6.Now we consider the discretization obtained by replacing both u x and u y by centraldifferences: ˜ u i +2 ,j − ˜ u i,j h − ˜ u i +1 ,j = 0 , ˜ u i,j +2 − ˜ u i,j h + ˜ u i,j +1 = 0 , (62)with step sizes h and h , respectively. Denote by ˜ f (cid:48)(cid:48) and ˜ f (cid:48)(cid:48) the left hand sides in (62). Thepassivity check yields the consequence (with underlined leaders) (cid:16) h σ ˜ f (cid:48)(cid:48) − h σ ˜ f (cid:48)(cid:48) − h h (cid:0) ˜ u i +2 ,j +1 + 2 h ˜ u i +1 ,j +1 + ˜ u i,j +1 (cid:1) σ ˜ f (cid:48)(cid:48) − h h (cid:16) ˜ u i +1 ,j +2 − h ˜ u i +1 ,j +1 + ˜ u i +1 ,j (cid:17) σ ˜ f (cid:48)(cid:48) + h ˜ f (cid:48)(cid:48) − h ˜ f (cid:48)(cid:48) (cid:17) / h h ˜ u i +1 ,j +1 (cid:16) ( h + h ) ˜ u i +1 ,j +1 − ˜ u i +1 ,j + ˜ u i,j +1 (cid:17) . The continuous limit of this difference polynomial ˜ f (cid:48)(cid:48) is given by u (cid:0) u y + u (cid:1) h h − u (cid:0) u x − u (cid:1) h h . f (cid:48)(cid:48) defined by σ ˜ f (cid:48)(cid:48) − h h (cid:104) ( h + h ) (cid:16) ˜ u i +2 ,j +1 + 4 h ˜ u i +1 ,j +1 + 4 h ˜ u i,j +1 ˜ u i +1 ,j +1 + ˜ u i,j +1 (cid:17) + ˜ u i +1 ,j +1 − ˜ u i +2 ,j (cid:105) (˜ u i +2 ,j +1 + 2 h ˜ u i +1 ,j +1 + ˜ u i,j +1 ) σ ˜ f (cid:48)(cid:48) + 2 h h (cid:0) h ˜ u i +1 ,j +1 + ˜ u i,j +1 (cid:1) ˜ f (cid:48)(cid:48) − h (cid:0) h ˜ u i +1 ,j +1 + ˜ u i,j +1 (cid:1) ˜ f (cid:48)(cid:48) = h h (cid:16) h ˜ u i +1 ,j +1 + ˜ u i,j +1 (cid:17) (cid:16) h ( h ˜ u i +1 ,j + h ˜ u i,j +1 ) ˜ u i +1 ,j +1 +˜ u i +1 ,j +1 + ( h + h ) ˜ u i,j +1 − h ˜ u i +1 ,j − ˜ u i,j (cid:17) , as well as the difference polynomial ˜ f (cid:48)(cid:48) defined by σ ˜ f (cid:48)(cid:48) − h h (cid:104) ( h + h ) (cid:16) ˜ u i +1 ,j +2 + 4 h ˜ u i +1 ,j +1 − h ˜ u i +1 ,j ˜ u i +1 ,j +1 + ˜ u i +1 ,j (cid:17) − ˜ u i +1 ,j +1 + ˜ u i,j +2 (cid:105) (˜ u i +1 ,j +2 − h ˜ u i +1 ,j +1 + ˜ u i +1 ,j ) σ ˜ f (cid:48)(cid:48) − h h (cid:0) h ˜ u i +1 ,j +1 − ˜ u i +1 ,j (cid:1) ˜ f (cid:48)(cid:48) − h (cid:0) h ˜ u i +1 ,j +1 − ˜ u i +1 ,j (cid:1) ˜ f (cid:48)(cid:48) = h h (cid:16) h ˜ u i +1 ,j +1 − ˜ u i +1 ,j (cid:17) (cid:16) h ( h ˜ u i +1 ,j + h ˜ u i,j +1 ) ˜ u i +1 ,j +1 +˜ u i +1 ,j +1 − ( h + h ) ˜ u i +1 ,j + 2 h ˜ u i,j +1 − ˜ u i,j (cid:17) . The continuous limit of the difference polynomial ˜ f (cid:48)(cid:48) + ˜ f (cid:48)(cid:48) is given by2 u u y ( u y + 2 u ) h h + 4 u ( u x + u y ) h h − u u x ( u y − u ) h h , whose Janet normal form modulo (55) is2 h h ( h − h ) ( h + h ) u . (63)Since the differential polynomial u is not in the radical differential ideal corresponding to(55), we conclude that (62) is not s-consistent with (55) unless h = h .Let h = h = h in (62). If one handles this case along the same lines as above, onemay encounter an enormous growth of expressions. We demonstrate here how to benefit fromapplying inverse shifts to difference polynomials when possible, i.e., when no indeterminateswith negative shifts are introduced by this process. Note that the perfect closure of thedifference ideal ˜ I generated by (62) contains the reflexive closure of ˜ I , i.e., all ˜ f ∈ ˜ R suchthat σ ˜ f ∈ ˜ I for some σ ∈ Mon(Σ).Denote that left hand sides of (62), for h = h = h , again by ˜ f (cid:48)(cid:48) and ˜ f (cid:48)(cid:48) . Similarly to theprevious discussion, the passivity check yields a difference polynomial˜ f (cid:48)(cid:48) := ˜ u i +1 ,j +1 (cid:16) h ˜ u i +1 ,j +1 − ˜ u i +1 ,j + ˜ u i,j +1 (cid:17) . The difference polynomial˜ f (cid:48)(cid:48) + 2 h ˜ u i +1 ,j +1 ( σ ˜ f (cid:48)(cid:48) − σ ˜ f (cid:48)(cid:48) ) = ˜ u i +1 ,j +1 (˜ u i +2 ,j +1 − ˜ u i +1 ,j +2 − h ˜ u i +1 ,j +1 )can be shifted back by one step in each of the two grid directions, producing˜ f (cid:48)(cid:48) := ˜ u i,j (˜ u i +1 ,j − ˜ u i,j +1 − h ˜ u i,j ) . u (cid:0) u x − u y − u (cid:1) h . Passivity checks yield˜ f (cid:48)(cid:48) := σ ˜ f (cid:48)(cid:48) − h ˜ u i +1 ,j ˜ f (cid:48)(cid:48) = − ˜ u i +1 ,j (˜ u i +1 ,j +1 − ˜ u i,j )and ˜ f (cid:48)(cid:48) := σ ˜ f (cid:48)(cid:48) + 2 h ˜ u i,j +1 ˜ f (cid:48)(cid:48) = ˜ u i,j +1 (˜ u i +1 ,j +1 − ˜ u i,j ) , whose continuous limits are given by ± u ( u x + u y ) h . We obtain the decomposition into simple difference systems ˜ f (cid:48)(cid:48) = ˜ u i,j +1 (˜ u i +1 ,j +1 − ˜ u i,j ) = 0 , ˜ f (cid:48)(cid:48) = ˜ u i,j (˜ u i +1 ,j − ˜ u i,j +1 − h ˜ u i,j ) = 0 , ˜ u i,j (cid:54) = 0 , ∨ ˜ u i,j = 0 , the first one confirming s-consistency of (62) with (55) provided h = h = h .Now we present our numerical experiments with difference equations in FDA (62) provided h = h = h . Again we perform the Taylor expansion of their left-hand sides up to terms oforder h g : = u x − u + ( u x,x + u x,y − uu x − uu y ) h − (cid:0) ( u x,x + 2 u x,y + u y,y ) u + ( u x + u y ) (cid:1) h + (cid:0) u x,x,x + u x,x,y + u x,y,y (cid:1) h + O ( h ) = 0 ,g : = u y + u + ( u y,y + u x,y + 2 uu x + 2 uu y ) h + (cid:0) ( u x,x + 2 u x,y + u y,y ) u + ( u x + u y ) (cid:1) h + (cid:0) u x,x,y + u x,y,y + u y,y,y (cid:1) h + O ( h ) = 0 . (64)Based on Eqs. (64), the modified PDE for scheme (62) with h = h = h reads g − ( g ) x h + (cid:18)
13 ( g ) x,x −
13 ( g ) x u − g u x − g u h (cid:19) h = u x − u + u h + O ( h ) = 0 ,g − ( g ) y h + (cid:18)
13 ( g ) y,y + 13 ( g ) y u + 13 g u y − g u h (cid:19) h = u y + u − u h + O ( h ) = 0 . (65)Thus, the scheme (62) has second order accuracy, and the last can be increased to fourthorder as follows: ˜ u i +2 ,j − ˜ u i,j h − ˜ u i +1 ,j − h ˜ u i +1 ,j = 0 , ˜ u i,j +2 − ˜ u i,j h + ˜ u i,j +1 + h ˜ u i,j +1 = 0 . (66)The numerical behavior of schemes (62) and (66) in the above described initial value problemfor (55) with the initial data defined by the exact solution (56) is shown in Fig. 7 and Fig. 8,respectively. One can see that the experimental numerical accuracy is scaled in accordancewith the theoretical accuracy h for (62) and h for (66).36 y Figure 7: Computed error for discretization (62), maximum error = 0.0026833687620488877 x y Figure 8: Computed error for discretization (66), maximum error = 1.0482200407964845e-0537 xample 6.
Let us consider another quasi-linear system of PDEs u x − u v = 0 ,u y + u v = 0 , u = u ( x, y ) , v = v ( x, y ) . (67)We define the differential polynomial ring R = K{ u, v } with commuting derivations ∂ x and ∂ y , and we use the elimination ranking (cid:31) on R satisfying . . . (cid:31) v yy (cid:31) v xy (cid:31) v xx (cid:31) v y (cid:31) v x (cid:31) v (cid:31) . . . (cid:31) u yy (cid:31) u xy (cid:31) u xx (cid:31) u y (cid:31) u x (cid:31) u The passivity check involves a pseudo-reduction of the second equation in (67), multiplied by u , modulo the derivative of the first equation with respect to y . Hence, the case distinction u = 0 ∨ u (cid:54) = 0 is made. If u = 0, then (67) reduces to (cid:40) u = 0 ,v y = 0 , (68)which is a simple differential system. If u (cid:54) = 0, then pseudo-reduction yields u ( v y + u v ) + ∂ y ( u x − u v ) = ( u − u y ) v + u x,y . Another pseudo-reduction of the last differential polynomial modulo the first equation in (67)gives the simple differential system (with underlined leaders) u (cid:54) = 0 ,u x − u v = 0 ,u u x,y + ( u − u y ) u x = 0 . (69)We investigate the system of difference equations ˜ u i +1 ,j − ˜ u i,j h − ˜ u i,j ˜ v i,j = 0 , ˜ v i,j +1 − ˜ v i,j h + ˜ u i,j ˜ v i,j = 0 , (70)which is obtained as discretization of (67) by replacing ∂ x u and ∂ y v by the correspondingforward differences, with step sizes h and h , respectively. Denote by ˜ f and ˜ f the lefthand sides in (70). We use the difference ranking (cid:31) (cf. Definition 11) on ˜ R = ˜ K{ ˜ u, ˜ v } thatcorresponds to the differential one, namely . . . (cid:31) ˜ v i +2 ,j (cid:31) ˜ v i,j +1 (cid:31) ˜ v i +1 ,j (cid:31) ˜ v i,j (cid:31) . . . (cid:31) ˜ u i +2 ,j (cid:31) ˜ u i,j +1 (cid:31) ˜ u i +1 ,j (cid:31) ˜ u i,j . The passivity check yields the consequence h ˜ u i,j +1 ˜ f + σ ˜ f = ( h ˜ u i,j −
1) ˜ u i,j +1 ˜ v i,j + ˜ u i +1 ,j +1 − ˜ u i,j +1 h . The continuous limit of this difference polynomial ˜ f is given by u x − u v . u i,j +1 does not vanish. If ˜ u i,j +1 = 0, then ˜ u i,j = 0,and we obtain the simple difference system (cid:40) ˜ u i,j = 0 , ˜ v i,j +1 − ˜ v i,j = 0 . Otherwise, we continue with the generic case by applying pseudo-reduction to ˜ f modulo ˜ f ,which yields the remainder˜ u i,j ˜ f + ( h ˜ u i,j −
1) ˜ u i,j +1 ˜ f = (˜ u i,j ˜ u i +1 ,j +1 + ( h ˜ u i,j −
1) ˜ u i,j +1 ˜ u i +1 ,j − h ˜ u i,j ˜ u i,j +1 ) /h . The continuous limit of this difference polynomial ˜ f is given by( u u x,y + ( u − u y ) u x ) h . We obtain the simple difference system h ˜ f = ˜ u i +1 ,j − ˜ u i,j − h ˜ u i,j ˜ v i,j = 0 ,h ˜ f = ˜ u i,j ˜ u i +1 ,j +1 + ( h ˜ u i,j −
1) ˜ u i,j +1 ˜ u i +1 ,j − h ˜ u i,j ˜ u i,j +1 = 0 , ˜ u i,j (cid:54) = 0 , which is s-consistent with (69).Alternatively, the same difference system (70) can be checked for s-consistency with (69)by using difference Gr¨obner bases. If we choose the monomial ordering (cid:23) = ≥ TOPdegrevlex (cf.Appendix A), then the leading monomial of the equations in (70) are the underlined ones in˜ f = (˜ u i +1 ,j − ˜ u i,j ) /h − ˜ u i,j ˜ v i,j , ˜ f = (˜ v i,j +1 − ˜ v i,j ) /h + ˜ u i,j ˜ v i,j . Reduction of the S-polynomial of ˜ f and ˜ f modulo (70) yields˜ v i,j +1 h ˜ f − ˜ u i +1 ,j h ˜ f + h h ˜ u i,j ˜ v i,j ( ˜ f + ˜ f ) − h ˜ v i,j ( ˜ f − ˜ f ) = 0 . Hence, by Proposition 2, ˜ f and ˜ f form a difference Gr¨obner basis (cf. Definition 14) of thedifference ideal defined by (70), confirming s-consistency of the difference approximation (70)with the PDE system (67) again (cf. Theorem 1).Next we consider the discretization obtained by replacing ∂ x u by the forward differenceas before and ∂ y v by the backward difference: ˜ u i +1 ,j − ˜ u i,j h − ˜ u i,j ˜ v i,j = 0 , ˜ v i,j +1 − ˜ v i,j h + ˜ u i,j +1 ˜ v i,j +1 = 0 , (71)again with step sizes h and h , respectively. Denote by ˜ f (cid:48) and ˜ f (cid:48) the left hand sides in (71).We use the degree-reverse lexicographical ranking (cid:31) with ˜ v (cid:31) ˜ u on the difference polynomialring ˜ R = ˜ K{ ˜ u, ˜ v } , namely . . . (cid:31) ˜ v i +2 ,j (cid:31) ˜ u i +2 ,j (cid:31) ˜ v i,j +1 (cid:31) ˜ u i,j +1 (cid:31) ˜ v i +1 ,j (cid:31) ˜ u i +1 ,j (cid:31) ˜ v i,j (cid:31) ˜ u i,j . f (cid:48) and ˜ f (cid:48) are ˜ u i +1 ,j and ˜ v i,j +1 , respectively. Since these involve differentindeterminates ˜ u and ˜ v , passivity is ensured, but a case distinction regarding the initial of ˜ f (cid:48) leads to a splitting. Thus a difference decomposition of (71) is ( h ˜ u i,j +1 + 1) ˜ v i,j +1 − ˜ v i,j = 0 , ˜ u i +1 ,j − ˜ u i,j − h ˜ u i,j ˜ v i,j = 0 ,h ˜ u i,j +1 + 1 (cid:54) = 0 , ∨ ˜ v i,j = 0 ,h ˜ u i,j + 1 = 0 . Since the continuous limits of the equations in the first simple system are in the radicaldifferential ideal corresponding to (67) due to w-consistency of (71), we conclude that (71) iss-consistent with (67). For h → (cid:23) = ≥ TOPdegrevlex with ˜ u (cid:31) ˜ v (cf. Appendix A). Then the leading monomials of ˜ f (cid:48) and ˜ f (cid:48) are ˜ u i +1 ,j and ˜ u i,j +1 ˜ v i,j +1 , respectively. The passivity check yields σ ˜ f (cid:48) − h ˜ v i +1 ,j +1 ( σ ˜ f (cid:48) + ˜ f (cid:48) )= ˜ u i,j +1 ˜ v i +1 ,j +1 − h ˜ v i,j +1 ˜ v i +1 ,j +1 − ˜ v i +1 ,j +1 − h ˜ v i,j ˜ v i +1 ,j +1 + ˜ v i +1 ,j h . The continuous limit of this difference polynomial ˜ f (cid:48) is given by v y + u v . A further reductionyields ˜ v i,j +1 ˜ f (cid:48) − ˜ v i +1 ,j +1 ˜ f (cid:48) = − ( h ˜ v i,j +1 ˜ v i +1 ,j +1 − h ˜ v i,j ˜ v i,j +1 ˜ v i +1 ,j +1 − ˜ v i,j ˜ v i +1 ,j +1 + ˜ v i,j +1 ˜ v i +1 ,j ) /h . The continuous limit of this difference polynomial ˜ f (cid:48) is given by (cid:0) v v x,y − v y ( v x + v ) (cid:1) h . Note that the coefficient of h is the linear combination( v ∂ x − v x − v ) ( v y + u v ) − v ( u x − u v )of the original equations in (67). The reduction of the final S-polynomial is h ˜ v i,j +1 ˜ f (cid:48) + h ˜ u i,j +1 ˜ f (cid:48) − ( h ˜ v i,j ˜ v i,j +1 + ˜ v i,j ) ˜ f (cid:48) +( h ˜ v i,j − h ˜ v i,j +1 + 1) ˜ f (cid:48) + ˜ v i +1 ,j ˜ f (cid:48) = 0 . Hence, we again conclude that (71) is s-consistent with (67).
We extended the notion of s(trong)-consistency for FDA (9), introduced and studied in [24, 18,25] for the Cartesian grids ( h = h = · · · = h n ), to the regular ones, where the grid spacings h i may be pairwise different. This notion for a finite difference discretization (9) of PDE (1)satisfying the condition (19) in Definition 5.1 means that any element ˜ f in the differenceideal [ ˜ F ], as well as any in its perfect closure (cid:74) ˜ F (cid:75) , after appropriate normalization (cf. (20)),40pproximates an element f ∈ (cid:74) F (cid:75) in the radical differential ideal. Thereby, the algebraicproperties of discrete (finite difference) equations, characterized by the perfect differenceideal they generate, mimic the algebraic properties of differential equations characterized bythe radical differential ideal generated by these equations.By using the method of difference Gr¨obner bases we derived a new s-consistent and con-servative FDA (34)–(35) to the three-dimensional incompressible Navier-Stokes equations.This discretization allows to solve numerically the last equations in the pressure-Poissonformulation when the pressure is determined from the Poisson pressure equation and thevelocities from the momentum equations. Our numerical experiments with the two-dimen-sional analogue (39) of the new scheme have clearly demonstrated its superiority over theother two-dimensional schemes (40)–(42). In particular, the scheme reveals, at the discretelevel, a surprisingly high accuracy preservation of the mass conservation law (continuity equa-tion). This law is satisfied by the initial condition, but is not employed in the subsequentconstruction of the numerical solution.In general, the techniques of difference Gr¨obner bases cannot be applied to the s-consistencyanalysis of FDA to nonlinear PDE systems, since termination of Algorithm 1 is not guaran-teed. Instead, the fully algorithmic triangular difference Thomas decomposition, designed lastyear in our conference paper [25] and described in Section 7, can be applied. Before its appli-cation, we suggest first to apply the differential Thomas decomposition to the input PDE (1).Each output subsystem is simple (cf. Definition 6), which, e.g., clarifies the arbitrariness ofpower series solutions to the system [35] and thus facilitates formulating well-posed initialvalue problems in the sense of Hadamard [28] (cf. [21, Example 14] for the case of the three-dimensional Navier-Stokes equations). Disjointness of the decomposition, i.e., partition of thesolution space by the output subsystems, allows to confine the investigation to the uniquesimple system admitting a solution of interest. After a discretization of the input simpledifferential system providing its w-consistency (cf. Definition 19) we apply Algorithm 4 to theFDA. For different simple systems different ways to discretize may be chosen. Algorithm 4 isthe main one, it provides the difference Thomas decomposition into quasi-simple subsystems(cf. Definition 22). It is based on two subalgorithms: Algorithm 2 performing difference auto-reduction and Algorithm 3 computing Janet normal forms of difference polynomials. Finally,Algorithm 5 performs the s-consistency analysis for the input FDA with the simple differen-tial system. Since the difference Thomas decomposition partitions the solution space of theFDA, s-consistency holds if and only if every difference equation in each output subsystemapproximates an element in the radical differential ideal generated by the elements in theinput simple differential system. In the recent paper [39] it is argued that if a differential(or difference) decomposition algorithm terminates on every input, then one can provide acomputable upper bound for the size of its output in terms of the input, i.e., an upper boundfor number of output subsystems, their order and degree. Because of the termination of bothdecomposition algorithms, the upper bound estimation approach of paper [39] is applicableto differential and difference Thomas decompositions.For illustration we applied both methods, the one based on difference Gr¨obner bases andthe one based on difference Thomas decomposition, to the s-consistency analysis of finitedifference discretizations of two first-order quasi-linear PDE systems (Section 7) with twoindependent variables. The first PDE system (55) is overdetermined and has a consequenceof the conservation law form ∂ x u + ∂ y u = 0 . (72)If one approximates the partial derivatives in Eqs. (55) by forward differences, then its differ-ence S -polynomial in the continuous limit yields the equation u = 0, which does not follow41rom Eqs. (55). Therefore, the FDA (57) is s-inconsistent. Another FDA (58) combiningthe forward and backward differences for derivatives yields an S -polynomial that shows thats-consistency is equivalent to h = h . The third discretization (62) based on approximationof both partial derivatives by central differences also has a passivity condition, whose con-tinuous limit (63) allows to conclude that h = h is a necessary condition for s-consistency.In case h = h = h subsequent passivity checks produce rather large expressions. How-ever, by applying backward shifts to intermediate difference polynomials when possible, thes-consistency analysis is drastically simplified, yielding differential polynomials as continuouslimits which occur in the left-hand sides of Eqs. (55) or their sum (72) or their difference.Hence, for Cartesian grid, i.e., for equisized grid spacings, both difference approximations (58)and (62) are s-consistent. This example shows that s-consistency may place constraints on thegrid spacings. Furthermore, for both s-consistent FDA (58) and (62) we constructed modifiedequations and applied them to analyze the actual accuracy of those FDA and to increase theiraccuracy. Additionally, we used the exact solution (56) to (55) for the numerical constructionof this solution, verifying experimentally the theoretically predicted accuracy.The second quasi-linear PDE system (67) of Section 8 has two dependent variables. First,we apply to this system the differential Thomas decomposition which splits Eqs. (67) into twosimple systems (68) and (69). It is easy to see that any (w-consistent) FDA to system (68)is s-consistent due to the lack of passivity conditions. As to FDA for (69) one can replacepartial derivatives by the corresponding forward differences. This produces a simple differencesystem providing an s-consistent approximation to (69). Moreover, we also established thecompatibility of the differential and the difference Thomas decomposition by starting with (67)and discretizing its equations by forward differences. The difference Thomas decompositionagain produces the same discrete version of (69). Alternatively, we applied the method ofdifference Gr¨obner bases to verify our results. Next we considered the discretization of (67) byusing the forward difference to approximate u x and the backward difference to approximate v y . We established the s-consistency by means of the Gr¨obner basis method using one of themonomial orderings described in Appendix A.Concerning implementation of (nonlinear) difference Gr¨obner basis construction, the onlyone is realized in [36, 22]. There the problem of computation in a difference polynomial ringis reduced to computation in the ring of commutative polynomials whose set of variables isextended with their shifts obtained by the action of the elements in (8). In this case, under theassumption of an admissible monomial ordering compatible with the order function definedin [22, Def. 4.1] and for difference ideals that admit finite Gr¨obner bases one can use thealgorithm designed in [22, Alg. 4.1] to compute such a basis in a finite number of steps.This algorithm, implemented in Maple [22], may cause a quite considerable growth of thenumber of variables involved in the computation, and thus restricts applicability to rathersmall problems. The difference Thomas decomposition has not yet been implemented. Allcomputations with difference polynomials presented in the paper were done “by hand” usingMaple for simplification of intermediate expressions. Acknowledgements
The work of V.P. Gerdt (Sections 1,3-5,9; Section 6, except Fig.1-5;Algorithm 5 of Section 7) is supported by the Russian Science Foundation under grant No.20-11-202574.
A Monomial ordering for difference Gr¨obner basis
Using the identification ˜ u ( r ) a ,a = σ a ˜ u ( r ) for a ∈ ( Z ≥ ) , 1 ≤ r ≤ m , we define a total ordering (cid:65) on the set of monomials in the infinitely many indeterminates σ a ˜ u ( r ) , where a ∈ ( Z ≥ ) ,42 ∈ { , } , as follows: d (cid:89) i =1 σ a i ˜ u ( r i ) (cid:65) e (cid:89) i =1 σ b i ˜ u ( s i ) : ⇐⇒ d (cid:88) i =1 | a i | > e (cid:88) i =1 | b i | or (cid:16) d (cid:88) i =1 | a i | = e (cid:88) i =1 | b i | and( ( a j , r j ) , . . . , ( a j d , r j d ) ) > lex ( ( b k , s k ) , . . . , ( b k e , s k e ) ) (cid:17) , where ( a j , r j ) ≥ ( a j , r j ) ≥ . . . ≥ ( a j d , r j d )and ( b k , s k ) ≥ ( b k , s k ) ≥ . . . ≥ ( b k e , s k e )are the tuples ( a i , r i ) and ( b i , s i ) arranged in decreasing order with respect to the ordering ≥ used for breaking ties, and where > lex is the lexicographic ordering comparing tuple entrieswith respect to ≥ . The ordering ( a j i , r j i ) ≥ ( b k i , s k i ) is assumed to respect addition of amulti-index c to a j i and b k i . In Section 8, where we let ˜ u = ˜ u (1) and ˜ v = ˜ u (2) , we choose ≥ = ≥ TOPdegrevlex , which is defined by(( a , a ) , r ) ≥ TOPdegrevlex (( a (cid:48) , a (cid:48) ) , r (cid:48) ) : ⇐⇒ a + a > a (cid:48) + a (cid:48) or (cid:16) a + a = a (cid:48) + a (cid:48) and (cid:0) a < a (cid:48) or a = a (cid:48) and r ≤ r (cid:48) (cid:1)(cid:17) . For elimination purposes one may choose ≥ = ≥ POTdegrevlex , which is defined by(( a , a ) , r ) ≥ POTdegrevlex (( a (cid:48) , a (cid:48) ) , r (cid:48) ) : ⇐⇒ r > r (cid:48) or (cid:16) r = r (cid:48) and (cid:0) a + a > a (cid:48) + a (cid:48) or a + a = a (cid:48) + a (cid:48) and a ≤ a (cid:48) (cid:1)(cid:17) . It is clear that we have ˜ t (cid:65) t (cid:54) = 1. Suppose that the differencemonomials ˜ v and ˜ w satisfy ˜ v (cid:65) ˜ w and let ˜ t be another difference monomial and θ ∈ Σ. Theneither the sum of shifts in ˜ v is greater than the sum of shifts in ˜ w , in which case the samestatement holds for ˜ t · θ ◦ ˜ v compared to ˜ t · θ ◦ ˜ w , or the sums of shifts are equal and, in the abovenotation, either the lexicographic ordering identifies an index i such that ( a j i , r j i ) > ( b k i , s k i )with respect to the ordering used for breaking ties, or ( ( b k , s k ) , . . . , ( b k e , s k e ) ) is a properprefix of ( ( a j , r j ) , . . . , ( a j d , r j d ) ). In the latter situations, application of θ is respected by ≥ , whereas multiplication by ˜ t leads to insertion of the pair corresponding to ˜ t at appropriatepositions in the above tuples, which is respected by the lexicographic ordering. Hence, weconclude ˜ t · θ ◦ ˜ v (cid:65) ˜ t · θ ◦ ˜ w in any case. Therefore, according to Definition 12, (cid:65) is an admissibledifference monomial ordering. 43 eferences [1] P. Amodio, Yu. A. Blinkov, V. P. Gerdt and R. La Scala. On consistency of finite differ-ence approximations to the Navier–Stokes equations . In: Computer Algebra in ScientificComputing / CASC 2013, Lecture Notes in Computer Science, vol. 8136, Springer, Cham(2013), pp. 46–60.[2] P. Amodio, Yu. A. Blinkov, V. P. Gerdt and R. La Scala.
Algebraic construction andnumerical behavior of a new s-consistent difference scheme for the 2D Navier–Stokesequations . Appl. Math. and Comput., (2017), 408–421.[3] D. N. Arnold, P. B. Bochev, R. B. Lehoucq, R. A. Nicolaides and M. Shashkov (eds.).Compatible Spatial Discretizations. Springer, 2006.[4] L. Beir˜ao da Veiga, K. Lipnikov and G. Manzini.
The Mimetic Finite Difference Methodfor Elliptic Problems . Modelling, Simulation & Applications, vol. 11. Springer, Cham(2014).[5] F. Boulier, D. Lazard, F. Ollivier and M. Petitot.
Computing representations for radicalsof finitely generated differential ideals . Appl. Algebra Eng. Commun. Comput., (2009),no. 1, 73–121.[6] T. B¨achler, V. Gerdt, M. Lange-Hegermann and D. Robertz. Algorithmic Thomas de-composition of algebraic and differential systems . J. Symb. Comput., (2012), no. 10,1233–1266.[7] T. Becker and V. Weispfenning. Gr¨obner Bases: A Computational Approach to Commu-tative Algebra . Graduate Texts in Mathematics, vol. 141. Springer, New York (1993).[8] Yu. A. Blinkov, C. F. Cid, V. P. Gerdt, W. Plesken and D. Robertz.
The MAPLEPackage Janet: II. Linear Partial Differential Equations . In: V. G. Ganzha, E. W. Mayrand E. V. Vorozhtsov (eds.), CASC 2003. Proc. 6th Int. Workshop on Computer Algebrain Scientific Computing, pp. 41–54. TU M¨unchen, 2003. Package
Janet is freely availableat http:///algebra.data.rwth-aachen.de/software/Janet [9] D. L. Brown, R. Cortez and M. L. Minion.
Accurate Projection Methods for the In-compressible NavierStokes Equations . Journal of Computational Physics, (2001),464–499.[10] J. E. Castillo and G. F. Miranda.
Mimetic Discretization Methods . CRC Press, BocaRaton (2013).[11] S. H. Christiansen, H. Z. Munthe-Kaas and B. Owren.
Topics in structure-preservingdiscretization . Acta Numerica, (2011), 1–119.[12] R. M. Cohn. Difference algebra . Interscience Publishers John Wiley & Sons, New York-London-Sydney (1965).[13] D. Cox, J. Little and D. O’Shea.
Ideals, Varieties and Algorithms. An Introductionto Computational Algebraic Geometry and Commutative Algebra . 3rd Edition. Springer,New York (2007).[14] V. Dorodnitsyn.
Applications of Lie groups to difference equations . Differential and In-tegral Equations and Their Applications, vol. 8, CRC Press, Boca Raton, FL (2011).4415] R. G. Erdmann.
Image-Based Numerical Simulation of Stokes Flow in Porous Media .Electronic PhD Dissertation. The University of Arizona, 2006. https://repository.arizona.edu/handle/10150/195724 [16] Xiao-Shan Gao, Zhang Huang and Chun-Ming Yuan,
Binomial difference ideals , J. Symb.Comput., (2017), 665–706.[17] V. P. Gerdt. Involutive Algorithms for Computing Gr¨obner Bases . Computational Com-mutative and Non-Commutative Algebraic Geometry. IOS Press, Amsterdam (2005), pp.199–225. arXiv:math.AC/0501111.[18] V. P. Gerdt.
Consistency Analysis of Finite Difference Approximations to PDESystems . Mathematical Modelling in Computational Physics / MMCP 2011, Lec-ture Notes in Computer Science, vol. 7125, pp. 28–42. Springer, Berlin (2012).arXiv:math.AP/1107.4269[19] V. P. Gerdt and Yu. A. Blinkov.
Involution and Difference Schemes for the Navier-StokesEquations . Computer Algebra in Scientific Computing / CASC 2009, Lecture Notes inComputer Science, vol. 5743, Springer-Verlag, Berlin (2009), pp. 46–60.[20] V. P. Gerdt, Yu. A. Blinkov and V. V. Mozzhilkin, Gr¨obner Bases and Genera-tion of Difference Schemes for Partial Differential Equations. SIGMA , 051 (2006)arXiv:math.RA/0605334[21] V. P. Gerdt, M. Lange-Hegermann and D. Robertz. The Maple package TDDS for com-puting Thomas decompositions of systems of nonlinear PDEs . Comput. Phys. Commun., (2019), 202–215. arXiv:physics.comp-ph/1801.09942[22] V. Gerdt and R. La Scala.
Noetherian quotients of the algebra of partial difference poly-nomials and Gr¨obner bases of symmetric ideals.
J. Algebra, (2015), 1233–1261.[23] V. P. Gerdt and D. Robertz.
Consistency of Finite Difference Approximations for LinearPDE Systems and its Algorithmic Verification . In: S. Watt (ed.). Proceedings of ISSAC2010, pp. 53–59. Association for Computing Machinery (2010).[24] V. P. Gerdt and D. Robertz.
Computation of difference Gr¨obner bases . Comput. Sc. J.Moldova, , 2(59) (2012), 203–226. Package LDA is freely available on the web page http:///algebra.data.rwth-aachen.de/software/Janet [25] V. P. Gerdt and D. Robertz.
Algorithmic Approach to Strong Consistency Analysis ofFinite Difference Approximations to PDE Systems . In: R. Bradford (ed.). Proceedings ofthe 2019 International Symposium on Symbolic and Algebraic Computation, 15-18 July2019, Beihang University, Beijing, China, pages 163–170 (2019).[26] G. H. Golub and C. F. Van Loan.
Matrix Computations . Johns Hopkins Studies in theMathematical Sciences, 4th Edition. Johns Hopkins University Press, Baltimore (2013).[27] J.-L. Guermond, P. Minev and J. Shen.
An Overview of Projection methods for incom-pressible flows . Computer Methods in Applied Mechanics and Engineering, (2006),6011–6045.[28] J. Hadamard.
Sur les problmes aux drives partielles et leur signification physique . Prince-ton University Bulletin, (1902), 9–52.4529] E. Hubert. Notes on Triangular Sets and Triangulation-Decomposition Algorithms. II:Differential Systems . In: F. Winkler and U. Langer (eds.) SNSC 2001, Lecture Notes inComputer Science, vol. 2630, pp. 40–87. Springer, Berlin (2001).[30] M. Janet.
Le¸cons sur les syst`emes d’´equations aux d´eriv´ees partielles . Cahiers Scien-tifiques, IV. Gauthier-Villars, Paris (1929).[31] M. Kalkbrener.
A generalized Euclidean algorithm for computing triangular representa-tions of algebraic varieties . J. Symb. Comput., (1993), no. 2, 143–167.[32] J. Kim and P. Moin. Application of a fractional-step method to incompressible Navier-Stokes equations . J. Comput. Phys., (1985), no. 2, 308–323.[33] E. R. Kolchin. Differential algebra and algebraic groups . Pure and Applied Mathematics,vol. 54. Academic Press, New York-London (1973).[34] B. Koren, R. Abgrall, P. Bochev, J. Frank and B. Perot (eds.).
Physics-compatible nu-merical methods . J. Comput. Phys., (2014), Part B, 1039–1526.[35] M. Lange-Hegermann.
The differential counting polynomial , Found. Comput. Math., (2018), no. 2, 291–308.[36] R. La Scala. Grbner bases and gradings for partial difference ideals , Math. Comput., (2015) 959–985.[37] F. Lemaire, M. Moreono Maza, W. Pan and Y. Xie. When does (cid:104) T (cid:105) equal sat( T )? J. Symb. Comput., (2011), no. 12, 1291–1305.[38] A. Levin. Difference Algebra . Algebra and Applications, vol. 8. Springer (2008).[39] W. Li, A. Ovchinnikov, G. Pogudin and T. Scanlon.
Algorithms yield upper bounds indifferential algebra . arXiv:math.AC/2005.01608[40] L. M. Milne-Thompson.
Theoretical Hydrodynamics . 5th Edition, Macmillan EducationLTD, Houndmills (1968).[41] P. Moin.
Fundamentals of Engineering Numerical Analysis . 2nd Edition. Cambridge Uni-versity Press, New York (2010).[42] Y. Ning, K. N. Premnath and D. V. Patil.
Numerical study of the properties of the centralmoment lattice Boltzmann method . Internat. J. Numer. Methods Fluids, (2015), no. 2,59–90.[43] F. Ollivier. Standard Bases of Differential Ideals . In: S. Sakata (ed.) AAECC-8. LectureNotes in Computer Science, vol. 508, pp. 304–321. Springer, London (1990).[44] P. J. Olver.
Applications of Lie groups to differential equations . Graduate Texts in Math-ematics, vol. 107, 2nd edition, Springer-Verlag, New York (1993).[45] A. Ovchinnikov, G. Pogudin and T. Scanlon.
Effective difference elimination and Null-stellensatz . arXiv:math.AG/1712.01412[46] G. Pogudin, T. Scanlon and M. Wibmer.
Solving difference equations in sequences: Uni-versality and Undecidability . arXiv:math.AG/1909.032394647] G. J. Reid, A. D. Wittkopf and A. Boulton.
Reduction of systems of nonlinear partialdifferential equations to simplified involutive forms . European J. Appl. Math., (1996),no. 6, 635–666.[48] D. Rempfer. On Boundary Conditions for Incompressible Navier-Stokes Problems . Appl.Mech. Rev., (2006), 107–125.[49] Ch. Riquier. Les syst`emes d’´equations aux d´eriv´ees partielles . Gauthiers-Villars, Paris(1910).[50] J. F. Ritt.
Differential Algebra . American Mathematical Society Colloquium Publica-tions, vol. 33, American Mathematical Society, New York (1950).[51] D. Robertz.
Formal Algorithmic Elimination for PDEs . Lecture Notes in Mathematics,vol. 2121. Springer, Cham (2014).[52] A. A. Samarskii.
Theory of Difference Schemes . Marcel Dekker, New York (2001).[53] W. M. Seiler.
Involution: The Formal Theory of Differential Equations and its Appli-cations in Computer Algebra . Algorithms and Computation in Mathematics, vol. 24.Springer (2010).[54] J. C. Strikwerda.
Finite Difference Schemes and Partial Differential Equations , 2nd Edi-tion. SIAM, Philadelphia (2004).[55] G. I. Taylor. On the decay of vortices in a viscous fluid.
Philosophical Magazine , vol. 46,671–674 (1923).[56] J. M. Thomas.
Differential Systems . AMS Colloquium Publications XXI (1937).[57] J. M. Thomas.
Systems and Roots . William Byrd Press, Richmond, VA (1962).[58] J. W. Thomas.
Numerical Partial Differential Equations: Finite Difference Methods .Springer-Verlag, New York (1995).[59] Dongming Wang.
Elimination Methods . Springer, Vienna (2001).[60] Dongming Wang.
Elimination practice. Software Tools and Applications . Imperial CollegePress, London (2004).[61] Wu Wen-tsun.