Algorithmic approach to strong consistency analysis of finite difference approximations to PDE systems
aa r X i v : . [ c s . S C ] A p r Algorithmic approach to strong consistency analysis of finitedifference approximations to PDE systems
Vladimir P. Gerdt
Joint Institute for Nuclear Research, Dubna, Russiaand Peoples’ Friendship University of Russia (RUDN)Moscow, [email protected]
Daniel Robertz
School of Computing, Electronics and MathematicsUniversity of PlymouthPlymouth, United [email protected]
ABSTRACT
For a wide class of polynomially nonlinear systems of partial differ-ential equations we suggest an algorithmic approach to the s(trong)-consistency analysis of their finite difference approximations onCartesian grids. First we apply the differential Thomas decomposi-tion to the input system, resulting in a partition of the solution set.We consider the output simple subsystem that contains a solutionof interest. Then, for this subsystem, we suggest an algorithm forverification of s-consistency for its finite difference approximation.For this purpose we develop a difference analogue of the differen-tial Thomas decomposition, both of which jointly allow to verifythe s-consistency of the approximation. As an application of ourapproach, we show how to produce s-consistent difference approx-imations to the incompressible Navier-Stokes equations includingthe pressure Poisson equation.
CCS CONCEPTS • Computing methodologies → Algebraic algorithms ; •
Math-ematics of computing → Partial differential equations ; Nonlinearequations ; Discretization ; KEYWORDS
Partial differential equations, Finite difference approximations, Con-sistency, Thomas decomposition
ACM Reference format:
Vladimir P. Gerdt and Daniel Robertz. 2019. Algorithmic approach to strongconsistency analysis of finite difference approximations to PDE systems. In
Proceedings of International Symposium on Symbolic and Algebraic Compu-tation, Beijing, China, July 15-18, 2019 (ISSAC ’19),
Except very special cases, partial differential equations (PDE) ad-mit numerical integration only. Historically first and one of themost-used numerical methods is finite difference method [1] based
Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full cita-tion on the first page. Copyrights for components of this work owned by others thanACM must be honored. Abstracting with credit is permitted. To copy otherwise, or re-publish, to post on servers or to redistribute to lists, requires prior specific permissionand/or a fee. Request permissions from [email protected].
ISSAC ’19, July 15-18, 2019, Beijing, China © 2019 Association for Computing Machinery.ACM ISBN ...$15.00https://doi.org/ on approximation of PDE by difference equations defined on a cho-sen solution grid. To construct a numerical solution, the obtained fi-nite difference approximation (FDA) to PDE is augmented with anappropriate discretization of initial or/and boundary condition(s)providing uniqueness of solution. As this takes place, the qualityof numericalsolution to PDEisdetermined by the quality of itsFDA.Any reasonable discretization must provide the convergence ofa numerical solution to a solution of PDE in the limit when the gridspacings tend to zero. However, except for a very limited class ofproblems, convergence cannot be directly established. In practice,for a given FDA, its consistency and stability are analyzed as thenecessary conditions for convergence. Consistency implies reduc-tion of the FDA to the original PDE when the grid spacings tend tozero and stability provides boundedness of the error in the solutionunder small perturbation in the numerical data.One of the most challenging problems is to construct FDA which,on the one hand, approximates the PDE and, on the other hand,mimics basic algebraic properties and preserves the algebraic struc-ture [2] of the PDE. Such mimetic or algebraic structure preservingFDA are more likely to produce highly accurate and stable numer-ical results (cf. [3]). In [4, 5], for polynomially nonlinear PDE sys-tems and regular solution grids, we introduced the novel conceptof strong consistency , or s-consistency , which strengthens the con-cept of consistency and means that any element of the perfect dif-ference ideal generated by the polynomials in FDA approximatesan element in the radical differential ideal generated by the poly-nomials in PDE. In the subsequent papers [6, 7], by computationalexperiments with two-dimensional incompressible Navier-Stokesequations, it was shown that s-consistent FDA have much betternumerical behavior than FDA which are not s-consistent.For linear PDE one can algorithmically verify [4] s-consistencyof their FDA. In the nonlinear case such verification [5] requiredcomputation of a difference Gröbner basis for FDA. Since differ-ence polynomial rings [8] are non-Noetherian, the difference Gröb-ner basis algorithms [5, 9] do not terminate in general. In compar-ison to differential algebra, fewer computational results have beenobtained in difference algebra. A decomposition technique was de-veloped only for binomial perfect difference ideals [10]. More gen-erally, in the present paper, a difference analogue of the differentialThomas decomposition [11–14] is obtained (see Section 6), whichprovides an algorithmic tool for s-consistency analysis of FDA tosimple PDE subsystems on Cartesian grids (see Section 7). In par-ticular, given an FDA to the momentum and continuity equationsin the Navier-Stokes PDE system for incompressible flow, our ap-proach derives an s-consistent approximation containing the pres-sure Poisson equation (see Section 9).
SSAC ’19, July 15-18, 2019, Beijing, China Vladimir P. Gerdt and Daniel Robertz
Completion to involution is the cornerstone of the differential Thomasdecomposition [11–14]. The underlying completion algorithm [15]is based on the theory of Janet division and Janet bases [13, 15,16] which stemmed from the Riquier-Janet theory [17, 18] of or-thonomic PDE. Joseph M. Thomas [14] generalized the Riquier-Janet theory to non-orthonomic polynomially nonlinear PDE andshowed how to decompose them into the triangular subsystemswith disjoint solution sets. Janet bases are Gröbner ones with addi-tional structure, and Wu Wen-tsun was the first who showed [19]that the Riquier-Janet theory can be used for algorithmic construc-tion of algebraic Gröbner bases. We dedicate this paper to com-memoration of his Centennial Birthday.
In the given paper we consider PDE systems of the form f = · · · = f s = 0 , F := { f , . . . , f s } ⊂ R , s ∈ Z ≥ , (1)where R := K{ u } is the ring of polynomials in the dependent vari-ables u := { u (1) , . . . , u ( m ) } and their partial derivatives obtainedfrom the operator power products in { ∂ , . . . , ∂ n } ( ∂ j = ∂ x j ). Weshall assume that coefficients of the polynomials are rational func-tions in a := { a , . . . , a l } , finitely many parameters (constants),over Q , i.e. K := Q ( a ). One can also extend the last field to Q ( a , x ),where x := { x , . . . , x n } is the set of independent variables. In thiscase we shall assume that coefficients of the differential polynomi-als in F do not vanish in the grid points defined in (2) below.To approximate (1) by a difference system we define a Cartesiancomputational grid (mesh) with spacing 0 < h ∈ R and fixed x by { ( x + k h , . . . , x n + k n h ) | k , . . . , k n ∈ Z } , (2)If the actual solution to (1) is u ( x ), then its approximation at thegrid nodes will be denoted by ˜u k ,..., k n ≈ u ( x + k h , . . . , x n + k n h ).Let ˜ K := Q ( a , h ) and ˜ R be the difference polynomial ring over˜ K , where ˜ K is a difference field of constants [8] with differencesΣ := { σ , . . . , σ n } acting on a grid function ˜ u ( α ) k ,..., k n as the shiftoperators σ ± i ˜ u ( α ) k ,..., k i ,..., k n = ˜ u ( α ) k ,..., k i ± ,..., k n , α ∈ { , . . . , m } . (3)The elements in ˜ R are polynomials in the dependent variables ˜ u ( α ) ( α = 1 , ..., m ) defined on the grid and in their shifts σ i ... σ i n n ˜ u ( α ) ( i j ∈ Z ). However, to provide termination of the decomposition al-gorithm of Sect. 6, we shall consider difference polynomials withnon-negative shifts only. We denote by Mon(Σ) the set of monomi-als in σ , ..., σ n . The coefficients of the polynomials are in ˜ K .The standard method to obtain FDA of such type to the differen-tial system (1) is replacement of the partial derivatives occurringin (1) by finite differences and application of appropriate powerproduct of the forward-shift operators in (3) to eliminate negativeshifts in indices which may come out of expressions like ∂ j u ( α ) ( x ) = u ( α ) k ,..., k j +1 ,..., k n − u ( α ) k ,..., k j − ,..., k n h + O ( h ) . Furthermore, the difference system˜ f = · · · = ˜ f s = 0 , ˜ F := { ˜ f , . . . , ˜ f s } ⊂ ˜ R , (4)is called an FDA to PDE (1) if it is consistent in accordance to:
Definition 2.1.
Given a PDE system (1), a difference system (4)is weakly consistent or w-consistent with (1) if( ∀ j ∈ { , . . . , s } ) [ ˜ f j −−−−→ h → f j ] . This is a universally adopted notion of consistency for a finitedifference discretization of PDE system (1) (cf. [20], Ch.7) and meansthat Eq. (4) reduces to Eq. (1) when the mesh step h tends to zero. Definition 2.2. [4] We say that a difference equation ˜ f ( ˜u ) = 0,˜ f ∈ ˜ R , implies the differential equation f ( u ) = 0, f ∈ R , and write˜ f ⊲ f , if the Taylor expansion of ˜ f about the grid point x , afterclearing denominators containing h , yields˜ f ( ˜u ) = h d f ( u ) + O ( h d +1 ) , d ∈ Z ≥ , (5)and O ( h d +1 ) denotes terms whose degree in h is at least d + 1. Remark 2.3.
Given ˜ f ( ˜ u ), computation of f ( u ) is straightforwardand has been implemented as routine ContinuousLimit in theMaple package
LDA [9, 21] (Linear Difference Algebra).
Definition 2.4. [5] FDA (4) to PDE system (1) is strongly consis-tent or s-consistent if( ∀ ˜ f ∈ n ˜ F o ) ( ∃ f ∈ n F o ) [ ˜ f ⊲ f ] . (6)Here n ˜ F o and n F o denote the perfect difference ideal generated by˜ F in ˜ R and the radical differential ideal generated by F in R . Remark 2.5.
It is clear that if condition (5) holds, then˜ f ( ˜u ) h d −−−−→ h → f ( u ) , (7)that is, ˜ f ( ˜u ) / h d approximates f ( u ). Accordingly, condition (6) meansthat, after clearing denominators, each element of n ˜ F o approxi-mates an element of n F o in the sense of (7). Lemma 2.6.
Let I = [ F ] be a differential ideal of R and ˜ I = [ ˜ F ] a difference ideal of ˜ R such that ( ∀ ˜ f ∈ ˜ I ) ( ∃ f ∈ I ) [ ˜ f ⊲ f ] . Then for the perfect closure n ˜ I o of ˜ I in ˜ R the condition (6) holds. Proof.
Let ˜ G be a (possibly infinite) reduced Gröbner basis of˜ I for an admissible monomial ordering ≻ (cf. [5]). Then˜ f = X ˜ д ∈ ˜ G ⊆ ˜ G X µ a ˜ д , µ σ µ ˜ д , a ˜ д , µ ∈ ˜ R , lm( a ˜ д , µ σ µ ˜ д ) (cid:22) lm( ˜ f ) . Here ˜ f ∈ ˜ I , ˜ G is a finite subset of ˜ G , lm denotes the leading mono-mial of its argument, and we use the multi-index notation µ := ( µ , . . . , µ n ) ∈ Z n ≥ , σ µ := σ µ · · · σ µ n n . In the continuous limit ˜ f implies the differential polynomial f := X д ∈ G X ν b д , ν ∂ ν д , b д , ν ∈ R , where ˜ G ⊲ G . Therefore, ˜ f ⊲ f ∈ [ F ] ⊆ n F o .Let now ˜ p ∈ n ˜ F o \ [ ˜ F ] and θ , . . . , θ r ∈ Mon(Σ) be such that˜ q := ( θ ˜ p ) k · · · ( θ r ˜ p ) k r ∈ [ ˜ F ] , k , . . . , k r ∈ Z ≥ . (8)From Eq. (8) it follows that ˜ q ⊲ q = p k + ··· + k r where ˜ p ⊲ p . Hence, p ∈ n F o . The perfect ideal n ˜ F o can be constructed [8] from [ ˜ F ] lgorithmic approach to strong consistency analysis for FDA to PDE systems ISSAC ’19, July 15-18, 2019, Beijing, China by the procedure in the form called shuffling and based on enlarge-ment of the generator set ˜ F with all polynomials ˜ p occurring in [ ˜ F ]in the form of Eq. (8) and on repetition of such enlargement. It isclear that each such enlargement of the intermediate ideals yieldsin the continuous limit a subset of n F o . (cid:3) The criterion of s-consistency is given by the following theorem.
Theorem 2.7. [5]
A difference approximation (4) to a differentialsystem (1) is s-consistent if and only if a reduced Gröbner basis ˜ G ⊂ ˜ R of the difference ideal [ ˜ F ] ⊂ ˜ R generated by ˜ F satisfies ( ∀ ˜ д ∈ ˜ G ) ( ∃ д ∈ n F o ) [ ˜ д ⊲ д ] . We recall the concept of Janet division. For details we refer to, e.g.,[13, Subsect. 2.1.1], [15], [16, Ch. 3].Let K be a field and R := K [ y , . . . , y n ] the commutative poly-nomial algebra over K with indeterminates y , . . . , y n . We denoteby Mon( R ) the set of monomials in y , . . . , y n and for a subset µ ⊆ { y , ..., y n } we define Mon( µ ) to be the subset of Mon( R ) con-sisting of the monomials involving only indeterminates from µ .If a term ordering on R is fixed and I is an ideal of R , then theset of leading monomials of non-zero polynomials in I are knownto form a set with the following property: Definition 3.1.
A set M ⊆ Mon( R ) is said to be Mon( R ) -multiple-closed if we have rm ∈ M for all m ∈ M and all r ∈ Mon( R ).The smallest Mon( R )-multiple-closed set in Mon( R ) containinga given set G ⊆ Mon( R ) is denoted by h G i . It is well known thatevery Mon( R )-multiple-closed set in Mon( R ) is finitely generatedin that sense and that it has a unique minimal generating set.We adopt Janet’s approach [18] of partitioning a Mon( R )-multiple-closed set M into finitely many subsets of the form Mon( µ ) m , where m ∈ M and µ = µ ( m , M ) ⊆ Mon( R ) (referred to as Janet division). Definition 3.2.
Let G ⊂ Mon( R ) be finite and m = y i · · · y i n n ∈ G .Then y k is said to be a multiplicative variable for m if and only if i k = max { j k | y j · · · y j n n ∈ G with j = i , . . . , j k − = i k − } . This yields a partition { y , . . . , y n } = µ ( m , G ) ⊎ µ ( m , G ), where theelements of µ ( m , G ) (resp. µ ( m , G )) are the multiplicative (resp. non-multiplicative) variables for m . The set G is Janet complete if h G i := [ m ∈ G Mon( R ) m = ] m ∈ G Mon( µ ( m , G )) m . Proposition 3.3.
For every
Mon( R ) -multiple-closed set M thereexists a finite Janet complete set J ⊂ Mon( R ) such that M = h J i . If G ⊂ Mon( R ) is finite, we call the minimal Janet complete set J ⊃ G such that h J i = h G i the Janet completion of G . It is obtainedalgorithmically by adding certain multiples of elements of G to G (which also proves Proposition 3.3), cf., e.g., [13, Algorithm 2.1.6]. Fundamental for both the differential Thomas decomposition (re-called in Sect. 5) as well as its difference analogue to be introducedin Sect. 6 is the Thomas decomposition of an algebraic system Sp = 0 , . . . , p s = 0 , p s +1 = 0 , . . . , p s + t = 0 ( s , t ∈ Z ≥ ) (9) where p , . . . , p s + t ∈ R := K [ z , . . . , z n ]. Here K is a field of char-acteristic zero with algebraic closure K , and R is the commutativepolynomial algebra over K with indeterminates z , . . . , z n . The so-lution set of the algebraic system S in (9) is defined to beSol K ( S ) := { a ∈ K n | p i ( a ) = 0 , p s + j ( a ) = 0 , i = 1 , ..., s , j = 1 , ..., t } . Assuming the indeterminates are ordered as in z ≻ z ≻ . . . ≻ z n ,a sequence of projections from K n is defined correspondingly by π i : K n → K n − i : ( a , ..., a n ) ( a i +1 , ..., a n ) , i = 1 , , ..., n − . For each p ∈ R \ K , this ordering defines the greatest indeterminateld( p ) occurring in p , referred to as leader , the coefficient init( p ) ofthe highest power of ld( p ) in p , called initial , and the discriminant disc( p ) := ( − d ( d − / res( p , ∂ p /∂ ld( p ) , ld( p )) / init( p ), where d isthe degree of p in ld( p ) and where res denotes the resultant. Definition 4.1.
An algebraic system S as in (9) is said to be simple if the following four conditions are satisfied.(1) None of p , . . . , p s , p s +1 , . . . , p s + t is constant.(2) The leaders of p , ..., p s , p s +1 , ..., p s + t are pairwise distinct.(3) For every r ∈ { p , . . . , p s , p s +1 , . . . , p s + t } , if ld( r ) = z k ,then the equation init( r ) = 0 has no solution in π k (Sol K ( S )).(4) For every r ∈ { p , . . . , p s , p s +1 , . . . , p s + t } , if ld( r ) = z k ,then the equation disc( r ) = 0 has no solution in π k (Sol K ( S )).(In (3) and (4), we have init( r ), disc( r ) ∈ K [ z k +1 , . . . , z n ].) Definition 4.2.
An algebraic system S as in (9) is said to be quasi-simple if conditions (1)–(3) (but not necessarily (4)) are satisfied.A Thomas decomposition of an algebraic system S as in (9) isa finite collection of simple algebraic systems S , ..., S r such thatSol K ( S ) = Sol K ( S ) ⊎ ... ⊎ Sol K ( S r ). It can be computed by an al-gorithm combining Euclidean pseudo-reduction and case distinc-tions. For details we refer to [12], [13, Subsect. 2.2.1], [22, Sect. 3.3]. A systemof polynomial partial differential equations and inequations f = 0 , . . . , f s = 0 , f s +1 = 0 , . . . , f s + t = 0 ( s , t ∈ Z ≥ ) (10)is given by elements f , . . . , f s + t of the differential polynomial ring R in u (1) , . . . , u ( m ) with commuting derivations ∆ := { ∂ , . . . , ∂ n } .For α ∈ { , . . . , m } , J ∈ ( Z ≥ ) n we identify u ( α ) J and ∂ J · · · ∂ J n n u ( α ) .Let Ω ⊆ R n be open and connected. The solution set of S on Ω isSol Ω ( S ) := { a = ( a , . . . , a m ) | a k : Ω → C analytic , k = 1 , . . . , m , f i ( a ) = 0 , f s + j ( a ) = 0 , i = 1 , . . . , s , j = 1 , . . . , t } . Definition 5.1. A ranking ≻ on R is a total ordering on the setMon(∆) u := { u ( α ) J | ≤ α ≤ m , J ∈ ( Z ≥ ) n } such that for all j ∈ { , . . . , n } , α , β ∈ { , . . . , m } , J , K ∈ ( Z ≥ ) n we have ∂ j u ( α ) ≻ u ( α ) and, if u ( α ) J ≻ u ( β ) K , then ∂ j u ( α ) J ≻ ∂ j u ( β ) K .A ranking ≻ is orderly if for all α , β ∈ { , . . . , m } , J , K ∈ ( Z ≥ ) n , J + · · · + J n > K + · · · + K n implies u ( α ) J ≻ u ( β ) K . Example 5.2.
Rankings ≻ TOP , lex and ≻ POT , lex on R are given by u ( α ) J ≻ TOP , lex u ( β ) K : ⇔ J ≻ lex K or ( J = K and α < β ) SSAC ’19, July 15-18, 2019, Beijing, China Vladimir P. Gerdt and Daniel Robertz and u ( α ) J ≻ POT , lex u ( β ) K : ⇔ α < β or ( α = β and J ≻ lex K ) , respectively, where ≻ lex compares multi-indices lexicographically.If a ranking ≻ on R is fixed, then for each f ∈ R \ K the leader , initial and discriminant of f are defined as in Section 4. Moreover,sep( f ) := ∂ f /∂ ld( f ) is called the separant of f .Janet division associates (with respect to a total ordering of ∆)to each f i = 0 with ld( f i ) = θ i u ( α ) the set µ i := µ ( θ i , G α ) ⊆ ∆ (resp. µ i := ∆ \ µ i ) of admissible (resp. non-admissible ) derivations , where G α := { θ ∈ Mon(∆) | θu α ∈ { ld( f ) , . . . , ld( f s ) } } . We call { f = 0 , ..., f s = 0 } or T := { ( f , µ ) , ..., ( f s , µ s ) } Janet com-plete if each G α equals its 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 alsosaid to be Janet reduced modulo T . Iterated pseudo-reductions of r modulo T yield its Janet normal form
NF( r , T , ≻ ), a Janet reduceddifferential polynomial, as explained in [13, Algorithm 2.2.45]. Definition 5.3.
Let T = { ( f , µ ) , . . . , ( f s , µ s ) } be Janet complete.Then { f = 0 , . . . , f s = 0 } or T is said to be passive , ifNF( ∂ f i , T , ≻ ) = 0 for all ∂ ∈ µ i = ∆ \ µ i , i = 1 , . . . , s . Definition 5.4.
Let a ranking ≻ on R and a total ordering on ∆be fixed. A differential system S as in (10) is said to be simple if thefollowing three conditions hold.(1) S is simple as an algebraic system (in the finitely manyindeterminates occurring in it, ordered by the ranking ≻ ).(2) { f = 0 , . . . , f s = 0 } is passive.(3) The left hand sides f s +1 , . . . , f s + t are Janet reduced modulothe passive differential system { f = 0 , . . . , f s = 0 } . Proposition 5.5 ([13], Prop. 2.2.50).
Let S be a simple differen-tial system, defined over R , as in (10). Let E be the differential idealof R which is generated by f , . . . , f s and let q be the product of theinitials and separants of all f , . . . , f s . 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 Janetnormal form of f modulo f , . . . , f s is zero.Definition 5.6. A Thomas decomposition of a differential system S as in (10) (with respect to ≻ ) is a finite collection of simple differ-ential systems S , ..., S r such that Sol Ω ( S ) = Sol Ω ( S ) ⊎ ... ⊎ Sol Ω ( S r ).For any differential system S as in (10) and any ranking ≻ on R a Thomas decomposition of S can be computed algorithmically.For more details we refer to, e.g., [12], [13, Subsection 2.2.2], [11]. A system ˜ S of polynomial partial difference equations and inequa-tions ˜ f = 0 , . . . , ˜ f s = 0 , ˜ f s +1 = 0 , . . . , ˜ f s + t = 0 ( s , t ∈ Z ≥ ) (11)is given by elements ˜ f , . . . , ˜ f s + t of the difference polynomial ring˜ R in ˜ u (1) , . . . , ˜ u ( m ) with commuting automorphisms Σ = { σ , . . . , σ n } . For α ∈ { , . . . , m } , J ∈ ( Z ≥ ) n we identify ˜ u ( α ) J and σ J · · · σ J n n ˜ u ( α ) .We denote by ˜ S = (resp. ˜ S = ) the set { ˜ f , ..., ˜ f s } (resp. { ˜ f s +1 , ..., ˜ f s + t } ).A ranking on ˜ R is defined in the same way as in Definition 5.1by replacing the action of ∂ i by the action of σ i and ∆ by Σ.For a subset L of ˜ R we denote by [ L ] the difference ideal of ˜ R generated by L . Let E be a difference ideal of ˜ R and ∅ 6 = Q ⊆ ˜ R bemultiplicatively 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 } . The first algorithm to be introduced performs an auto-reductionof a finite set of difference polynomials.
Algorithm 1:
Auto-reduce for difference algebra
Input: L ⊂ ˜ R \ ˜ K finite and a ranking ≻ on ˜ R such that L = ˜ S = for some finite difference system ˜ S which isquasi-simple as an algebraic system (in the finitely manyindeterminates ˜ u ( α ) J which occur in it, totally ordered by ≻ ) Output: a ∈ { true , false } and L ′ ⊂ ˜ R \ ˜ K finite such that[ L ′ ] : 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, in case a = true , there exist no ˜ f , ˜ f ∈ L ′ , ˜ f = ˜ f , suchthat we have v := ld( ˜ f ) = θ ld( ˜ f ) for some θ ∈ Mon(Σ) anddeg v ( ˜ f ) ≥ deg v ( θ ˜ f ) L ′ ← L while ∃ ˜ f , ˜ f ∈ L ′ , ˜ f = ˜ f and θ ∈ Mon(Σ) such that we have v := ld( ˜ f ) = θ ld( ˜ f ) and deg v ( ˜ f ) ≥ deg v ( θ ˜ f ) do L ′ ← L ′ \ { ˜ f } ; v ← ld( ˜ f ) ˜ r ← init( θ ˜ f ) ˜ f − init( ˜ f ) v d θ ˜ f , d :=deg v ( ˜ f ) − deg v ( θ ˜ f ) if ˜ r = 0 then return ( false , L ′ ∪ { ˜ r } ) return ( true , L ′ )Since leaders are dealt with in decreasing order with respect to ≻ , and no ranking admits infinitely decreasing chains, Algorithm 1 terminates . Its correctness follows from the definition of E : Q .Janet division associates (with respect to a total ordering of Σ)to each ˜ f i = 0 with ld( ˜ f i ) = θ i ˜ u ( α ) the set µ i := µ ( θ i , ˜ G α ) ⊆ Σ (resp. µ i := Σ \ µ i ) of admissible ( 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 com-plete if each ˜ G α equals its 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 isalso said to be Janet reduced modulo T . Iterated pseudo-reductions lgorithmic approach to strong consistency analysis for FDA to PDE systems ISSAC ’19, July 15-18, 2019, Beijing, China of ˜ r modulo T yield its Janet normal form
NF(˜ r , T , ≻ ), which is theJanet reduced difference polynomial ˜ r ′ returned by Algorithm 2. Algorithm 2:
Janet-reduce for difference algebra
Input: ˜ r ∈ ˜ R , T = { ( ˜ f , µ ) , ( ˜ f , µ ) , . . . , ( ˜ f s , µ s ) } , and aranking ≻ on ˜ R , where T is Janet complete (with respect to ≻ ) Output: (˜ r ′ , b ) ∈ ˜ R × ˜ R such that (1) if ˜ r ∈ ˜ K or T = ∅ , then˜ r ′ = ˜ r , b = 1, (2) otherwise ˜ r ′ is Janet-reduced modulo T and˜ r ′ + [ ˜ f , . . . , ˜ f s ] = b · ˜ r + [ ˜ f , . . . , ˜ f s ] , where b is in the multiplicatively closed set generated by s [ i =1 { θ init( ˜ f i ) | θ ∈ Mon(Σ) , ld(˜ r ) ≻ θ ld( ˜ f i ) } ∪ { } ˜ r ′ ← ˜ r ; b ← if ˜ r ′ ˜ K then v ← ld(˜ r ′ ) while ˜ r ′ ˜ K and there exist ( ˜ f , µ ) ∈ T and θ ∈ Mon( µ ) such that v = θ ld( ˜ f ) and deg v (˜ r ′ ) ≥ deg v ( θ ˜ f ) do ˜ r ′ ← init( θ ˜ f ) ˜ r ′ − init(˜ r ′ ) v d θ ˜ f , d :=deg v (˜ r ′ ) − deg v ( θ ˜ f ) b ← init( θ ˜ f ) · b for each coefficient ˜ c of ˜ r ′ (as a polynomial in v ) do (˜ r ′′ , b ′ ) ← Janet-reduce (˜ c , T , ≻ ) replace the coefficient b ′ · ˜ c in b ′ · ˜ r ′ with ˜ r ′′ andreplace ˜ r ′ with this result b ← b ′ · b return (˜ r ′ , b )Algorithm 2 terminates because each coefficient ˜ c of ˜ r ′ is eitherconstant or has a leader which is smaller than ld(˜ r ′ ) with respectto ≻ , and a ranking ≻ does not allow infinitely decreasing chains. Correctness of the algorithm is clear.
Definition 6.1.
Let T = { ( ˜ f , µ ) , . . . , ( ˜ f s , µ s ) } be Janet complete.Then { ˜ f = 0 , . . . , ˜ f s = 0 } or T is said to be passive , ifNF( σ ˜ f i , T , ≻ ) = 0 for all σ ∈ µ i = Σ \ µ i , i = 1 , . . . , s . Definition 6.2.
Let a ranking ≻ on ˜ R and a total ordering on Σbe fixed. A difference system ˜ S as in (11) is said to be simple (resp., quasi-simple ) if the following three conditions hold.(1) ˜ S is simple (resp., quasi-simple) as an algebraic system (inthe finitely many occurring indeterminates, ordered by ≻ ).(2) { ˜ f = 0 , . . . , ˜ f s = 0 } is passive.(3) The left hand sides ˜ f s +1 , . . . , ˜ f s + t are Janet reduced modulothe passive difference system { ˜ f = 0 , . . . , ˜ f s = 0 } . Proposition 6.3.
Let ˜ S be a quasi-simple difference system over ˜ R as in (11). Let E be the difference ideal of ˜ R generated by ˜ f , . . . , ˜ f s and let Q be the smallest subset of ˜ R which is multiplicatively 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. Proof.
By definition of E : Q , every element ˜ f ∈ ˜ R for whichAlgorithm 2 yields Janet normal form zero is an element of E : Q .Let ˜ f ∈ E : Q , ˜ f = 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 X i =1 k i X j =1 c i , j α i , j ( ˜ f i ) . (12)Among all pairs ( i , j ) for which α i , j involves a non-admissible auto-morphism for ˜ f i = 0 let the pair ( i ⋆ , j ⋆ ) be such that α i ⋆ , j ⋆ (ld( ˜ f i ⋆ ))is maximal with respect to the ranking ≻ . Let σ be a non-admissibleautomorphism for ˜ f i ⋆ = 0 which divides the monomial α i ⋆ , j ⋆ .Since { ˜ f = 0 , ..., ˜ f s = 0 } is passive, there exist b ∈ Q , l , . . . , l s ∈ Z ≥ and d i , j ∈ ˜ R \ { } and β i , j ∈ Mon(Σ), j = 1, . . . , l i , i = 1, . . . , s , such that b · ( σ ˜ f i ⋆ ) = s X i =1 l i X j =1 d i , j β i , j ( ˜ f i ) , where each β i , j involves only admissible automorphisms for ˜ f i = 0.Let γ i ⋆ , j ⋆ := α i ⋆ , j ⋆ / σ and multiply (12) by γ i ⋆ , j ⋆ ( b ) to obtain γ i ⋆ , j ⋆ ( b ) · q ˜ f = s X i =1 k i X j =1 c i , j · γ i ⋆ , j ⋆ ( b ) · α i , j ( ˜ f i ) . In this equation we replace γ i ⋆ , j ⋆ ( b ) · α i ⋆ , j ⋆ ( ˜ f i ⋆ ) = γ i ⋆ , j ⋆ ( b · σ ( ˜ f i ⋆ ))by γ i ⋆ , j ⋆ s X i =1 l i X j =1 d i , j β i , j ( ˜ f i ) ! . Since γ i ⋆ , j ⋆ β i ⋆ , j ⋆ involves fewer non-admissible automorphismsfor ˜ f i = 0 than α i ⋆ , j ⋆ , iteration of this substitution process willrewrite equation (12) in such a way that every α i , j (ld( ˜ f i )) involv-ing non-admissible automorphisms for ˜ f i = 0 will be less than α i ⋆ , j ⋆ (ld( ˜ f i ⋆ )) with respect to ≻ . A further iteration of this substi-tution process will therefore produce an equation as (12) with no α i , j involving any non-admissible automorphisms for ˜ f i = 0.This shows that for every ˜ f ∈ ( E : Q ) \ { } there exists a Janetdivisor of ld( ˜ f ) in the passive set defined by ˜ f = 0, . . . , ˜ f s = 0. (cid:3) Let Ω ⊆ R n be open and connected and fix x ∈ Ω. Denoting thegrid in (2) by Γ x , h , we define F Ω , x , h := { ˜ u : Γ x , h ∩ Ω → C | ˜ u is the restriction to Γ x , h ∩ Ω ofsome locally analytic function u on Ω } , and for a system ˜ S of partial difference equations and inequationsas in (11) we define the solution setSol Ω , x , h ( ˜ S ) := { ˜ u ∈ F Ω , x , h | ˜ f i ( ˜ u ) = 0 , ˜ f s + j ( ˜ u ) = 0 for all i = 1 , . . . , s , j = 1 , . . . , t } . Definition 6.4.
Let ˜ S be a finite difference system over ˜ R and ≻ a ranking on ˜ R . A difference decomposition of ˜ S is a finite collec-tion of quasi-simple difference systems ˜ S , . . . , ˜ S r over ˜ R such thatSol Ω , x , h ( ˜ S ) = Sol Ω , x , h ( ˜ S ) ⊎ . . . ⊎ Sol Ω , x , h ( ˜ S r ). SSAC ’19, July 15-18, 2019, Beijing, China Vladimir P. Gerdt and Daniel Robertz
In the following algorithm,
Decompose in step 11 refers to analgorithm which computes a smallest superset of G = { ˜ f , . . . , ˜ f s } in ˜ R that is Janet complete as defined on page 4 (see also Section 3). Algorithm 3:
DifferenceDecomposition
Input:
A finite difference system ˜ S over ˜ R , a ranking ≻ 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 asan algebraic system, into quasi-simple systems (cf.Sect. 4) for i = 1 , . . . , r do if A i = ∅ then // no equation and no inequation return {∅} else ( a , G ) ← Auto-reduce ( A = i , ≻ ) // Alg. 1 if a = true then J ← Decompose ( G ) P ← { NF( σ ˜ f , J , ≻ ) | ( ˜ f , µ ) ∈ J , σ ∈ µ } // Alg. 2 if P ⊆ { } then // J is passive replace each ˜ д =0 in A i with NF( ˜ д , J , ≻ ) = 0 if A = i then insert { ˜ f = 0 | ( ˜ f , µ ) ∈ J } ∪ { ˜ д = 0 | ˜ д ∈ A = i } into T else if P ∩ ˜ K ⊆ { } then insert { ˜ f = 0 | ( ˜ f , µ ) ∈ J } ∪ { ˜ f = 0 | ˜ f ∈ P \ { }} ∪ { ˜ д = 0 | ˜ д ∈ A = i } into Q else insert { ˜ f = 0 | ˜ f ∈ G }∪{ ˜ д = 0 | ˜ д ∈ A = i } into Q until Q = ∅ Theorem 6.5.
Algorithm 3 terminates and is correct.
Proof.
Algorithm 3 maintains a set Q of difference systemsthat still have to be dealt with. Given that termination of all sub-algorithms has been proved, termination of Algorithm 3 is equiva-lent to the condition that Q = ∅ holds after finitely many steps.Apart from step 1, new systems are inserted into Q in steps 18and 20. We consider the systems that are at some point an elementof Q as the vertices of a tree. The root of this tree is the input sys-tem ˜ S . The systems which are inserted into Q in steps 18 and 20are the vertices of the tree whose ancestor is the system L that wasextracted from Q in step 3 which in the following steps producedthese new systems. Since the for loop beginning in step 5 termi-nates, the degree of each vertex in the tree is finite. We claim thatevery branch of the tree is finite, i.e., that the tree has finite height,hence, that the tree has only finitely many vertices. In case of step 20 the new system contains an equation which re-sulted from a non-trivial difference reduction in step 9. When thisnew system will be extracted from Q in a later round, a decom-position into quasi-simple algebraic systems will be computed instep 4. This may produce new branches of the tree, but along anyof these branches, after finitely many steps the condition a = truein step 10 will hold, because the order of the shifts in leaders ofthe arising equations is bounded by the maximum order of shiftsin leaders of the ancestor system L .In case of step 18 we are going to show that after finitely manysteps a difference equation is obtained whose leader has not shownup as a leader of an equation in any preceding system in the cur-rent branch of the tree. First of all, the passivity check (step 12)yielded an equation ˜ f = 0, ˜ f ∈ P \ ˜ K , which is Janet reducedmodulo J . Hence, either ld( ˜ f ) is not contained in the multiple-closed set generated by ld( G ), or there exists ( ˜ f ′ , µ ′ ) ∈ J such thatld( ˜ f ′ ) is a Janet divisor of ld( ˜ f ), but the degree of ˜ f in ld( ˜ f ) issmaller than the degree of ˜ f ′ in ld( ˜ f ′ ). In the first case the aboveclaim holds. The second case cannot repeat indefinitely: First ofall, if ld( ˜ f ) = ld( ˜ f ′ ), then in a later round, either a pseudo-reduc-tion of ˜ f ′ modulo ˜ f will be performed if the initial of ˜ f does notvanish, or init( ˜ f ) = 0 has been added as a new equation (withlower ranked leader). Since this leads to a sequence in Mon(Σ)which strictly decreases, infinite chains are excluded in this sit-uation. If case ld( ˜ f ) = ld( ˜ f ′ ) occurs repeatedly, 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, termina-tion follows from Dickson’s Lemma.In order to prove correctness, we note that a difference systemis only inserted into T if step 12 confirmed passivity. Such a sys-tem is quasi-simple as an algebraic system because (up to auto-reduction in step 9 and Janet completion in step 11) it was returnedas one system A i in step 4. Condition (3) in Definition 6.2 is ensuredby step 14. Hence, all difference systems in T are quasi-simple.Splittings of systems only arise in step 4 by adding an equationinit( ˜ f ) = 0 and the corresponding inequation init( ˜ f ) = 0, respec-tively, to the two new systems replacing the given one. Since nosolutions are lost or gained, this leads to a partition as required byDefinition 6.4. (cid:3) Recall that ˜ S = (resp. ˜ S = ) denotes the set of left hand sides of equa-tions (resp. inequations) in a difference system ˜ S . We shall use thesame notation for differential systems.Clearly, if one approximates the partial derivatives occurringin a simple differential system S by appropriate finite differences,then one obtains a w-consistent approximation ˜ S to S (cf. Sect. 2and 9).The following algorithm verifies s-consistency of such FDA. Correctness of the algorithm follows from Definition 2.1 (extendedto inequations) and from passivity of the output subsystems of Al-gorithm 3. Their solution spaces partition the solution space of theinput FDA. Thereby, any subsystem ˜ L i in the output with b i = true lgorithmic approach to strong consistency analysis for FDA to PDE systems ISSAC ’19, July 15-18, 2019, Beijing, China Algorithm 4:
S-ConsistencyCheck
Input:
A simple differential system S over R , a differentialranking ≻ on R , a difference ranking > on ˜ R , a totalordering on Σ (used by Decompose ) and a difference 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 = i s.t. ˜ f ⊲ f ∈ n S = o then // Def. 2.2 ˜ L ← ˜ L \ { ˜ L i } else b i ← true for ˜ f ∈ ˜ L = do compute f ∈ R such that ˜ f ⊲ f // Rem. 2.3 if NF( f , S = , ≻ ) = 0 then // Alg. 2 b i ← false ; break return { ( ˜ L i , b i ) | ˜ L i ∈ ˜ L } 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.
Example 8.1.
We consider the system of nonlinear PDEs ( ∂ u ∂ x − u = 0 ∂ u ∂ y + u = 0 , u = u ( x , y ) , (13)which is a simple differential system, as it is easily checked that thecross-derivative ∂ y ( u x − u ) − ∂ x ( u y + u ) reduces to zero modulo(13). We investigate the discretized system which is obtained by re-placing ∂ x and ∂ y by the forward differences D +1 , D +2 , respectively: ( D +1 ˜ u − ˜ u = 0 ( A ) D +2 ˜ u + ˜ u = 0 ( B ) (14)This system of nonlinear difference equations is simple as an alge-braic system, but the passivity check reveals the consequence σ A − σ B + ( h ˜ u i +1 , j + h ˜ u i , j + h ˜ u i , j − A +( h ˜ u i , j +1 − h ˜ u i , j + h ˜ u i , j + 1) B = − h ˜ u i , j . The continuous limit of ˜ u i , j for h → u , which is not in the radical differential ideal corresponding to(13). Hence, FDA (14) is not s-consistent with system (13).Now we consider the discretization obtained by replacing ∂ x and ∂ y by D +1 as before and the backward difference D − , respectively: ( D +1 ˜ u − ˜ u = 0 ( C ) D − ˜ u + ˜ u = 0 ( E ) (15) In order to avoid negative shifts, we replace equation ( E ) by σ ( E ).Then this system of nonlinear difference equations is simple be-cause it is algebraically simple and the passivity check yields σ E − ( h ˜ u i , j +1 + h ˜ u i +1 , j +1 + h ˜ u i , j +1 + 1) σ C + C − E − h ( h ˜ u i , j +1 + ˜ u i , j +1 + ˜ u i , j ) E = 0 . We conclude that FDA (15) is s-consistent with system (13).
Example 9.1.
The Navier-Stokes equations for a three-dimensio-nal incompressible viscous flow in vector notation are ∂ u ∂ t + ( u · ∇ ) u + ∇ p − u = 0 , ∇ · u = 0 , (16)where x = ( x , x , x ), u ( x , t ) is the velocity vector u = ( u , v , w ), p ( x , t ) is the pressure and Re is the Reynolds number. For the rank-ing ≻ TOP , lex (Example 5.2) such that ∂ t ≻ ∂ ≻ ∂ ≻ ∂ and p ≻ u ≻ v ≻ w , (17)the (non-admissible) prolongation ∇ · ∂ t u = 0 of the right (continu-ity) equation in (16) and its reduction modulo the left (momentum)equation yields the pressure Poisson equation∆ p + ∇ · ( u · ∇ ) u = 0 , (18)which is the integrability condition (cf. [16], p.50) to (16). Clearly,the differential system (16) and (18) satisfies the simplicity condi-tions (1)–(4) in Definition 4.1. Now we consider the following classof FDA to (16) defined on the four-dimensional grid (2) D t ˜u + ( ˜u · D ) ˜u + D ˜ p − ˜u = 0 , D · ˜u = 0 , (19)where D t approximates ∂ t , D = ( D , D , D ) approximates ∇ and˜∆ approximates ∆. It is clear that system (19) is w-consistent with(16). If one considers the difference analogue of ranking (17) satis-fying σ t ≻ σ ≻ σ ≻ σ and ˜ p ≻ ˜ u ≻ ˜ v ≻ ˜ w , (20)then completion of (19) to a passive form by Algorithm 3 is equiva-lent to enlargement of this system with the integrability condition( D · D ) p + D · ( ˜u · D ) ˜u = 0 . (21)Eq. (21) approximates Eq. (18) and can be obtained, in the full anal-ogy with the differential case, by the prolongation D · D t ˜ u = 0 ofthe discrete continuity equation in system (19) and its reductionby the discrete momentum equation.The left-hand sides of Eqs. (16) and (18) form a difference Gröb-ner basis of the ideal generated by Eqs. (19) in Q (Re , h) { ˜u , ˜ p } . Hence,by Theorem 2.7, FDA (19)–(21) to Eqs. (16), (18) is s-consistent. Remark 9.2.
Formulae (19) and (21) give s-consistent FDA of theNavier-Stokes and pressure Poisson equations in the two-dimen-sional case as well. Examples of such FDA were studied in [6]. Onemore s-consistent two-dimensional FDA was derived in [7]. In itsapproximation of Eq. (18) the redundant to zero term − ˜∆ ( ∇· ˜u )Re wasincluded in the left-hand side of (21). This inclusion improves thenumerical behavior of FDA (cf. [23], Sect.3.2). SSAC ’19, July 15-18, 2019, Beijing, China Vladimir P. Gerdt and Daniel Robertz
Example 9.3.
For the two-dimensional system (16), (18) with gridvelocities ( u , v ) and pressure p we consider the discretization ˜ e (1) := D ˜ u + D ˜ v = 0 , ˜ e (2) := D t ˜ u + ˜ uD ˜ u + ˜ vD ˜ u + D ˜ p − ˜∆ ˜ u = 0 , ˜ e (3) := D t ˜ v + ˜ uD ˜ v + ˜ vD ˜ v + D ˜ p − ˜∆ ˜ v = 0 , ˜ e (4) := ˜∆ ˜ p + ( D ˜ u ) + 2( D ˜ u ) ( D ˜ v ) + ( D ˜ v ) = 0 , (22)where D t = σ t − h , D i = σ i − σ − i h , ˜∆ = σ + σ − σ − + σ − h and i ∈ { , } . Then FDA (22) is w-consistent with (16) and (18).However, it is s-inconsistent since ˜ e (4) n ˜ I o where n ˜ I o ⊂ ˜ R is the perfect closure (see Lemma 2.6) of the ideal generated by { ˜ e (1) , ˜ e (2) , ˜ e (3) } . It follows, as modulo n ˜ I o the equality holds D · ( ˜u · D ) ˜u = ( D ˜ u ) + 2( D ˜ u ) ( D ˜ v ) + ( D ˜ v ) , whereas the difference operator D · D in (21) is not equal to ˜∆: D · D = σ + σ − σ − + σ − h = ˜∆ .
10 CONCLUSIONS
In this paper, for the first time, we devised a universal algorith-mic approach to check s(trong)-consistency of a system of finitedifference equations that approximates a polynomially nonlinearPDE system on a Cartesian solution grid. In our earlier paper [4]we studied this problem for linear PDE systems and showed howto check their s-consistency by using differential and differenceGröbner bases of ideals generated by the polynomials in PDE andFDA. As this takes place, all related computations can be done, forexample, with the Maple packages
LDA [21] and
Janet [24].Extension of the Gröbner basis method to the nonlinear case isnot algorithmic due to the non-Noetherity of differential and differ-ence polynomial rings. On the other hand, the differential Thomasdecomposition (Def. 5.6) and its difference analogue (Def. 6.4) arefully algorithmic (cf. [11–13] and Alg. 3). These decompositionsare essentials of the s-consistency check (Alg. 4). The differentialThomas decomposition is built into Maple 2018 and its implemen-tation for previous versions of Maple is freely available on the web.Algorithm 3 has not been implemented yet.If we are looking for s-consistent FDA to a simple PDE systemand for a (w-consistent) FDA Algorithm 4 returns false , as it takesplace in Example 9.3, then we have to try another FDA and checkthe s-consistency again. In doing so, if we know a minimal gen-erating set for the radical differential ideal generated by the inputsimple differential system, then its FDA should be tried as an inputfor Algorithm 3. Such is indeed the case with the Navier-Stokesequations (Ex. 9.1), for which Algorithm 3 returns s-consistent dis-cretization (19), (21) if it is applied to Eqs. (16) and ranking (20).However, the choice of FDA to the minimal generating set forthe simple differential system as an input for Algorithm 3 not al-ways yields s-consistent FDA, as demonstrated by Example 8.1. Inaddition, designing an algorithm for construction of a minimal gen-erating set for an ideal is an open problem for commutative poly-nomial rings and is probably unsolvable in the differential case.In applications of finite difference methods to PDE systems whichhave integrability conditions, it is important not only to preserve these conditions at the discrete level, but to ensure also that FDAis s-consistent with the PDE system. FDA (19), (21) to the Navier-Stokes equations (16) satisfies this requirement and for this reasonit is appropriate for numerical solution of initial or/and boundary-value problems for (16) in the velocity-pressure formulation.
11 ACKNOWLEDGMENTS
The authors are grateful to the referees for their valuable remarks.The contribution of the first author (V.P.G.) was partially supportedby the Russian Foundation for Basic Research (grant No. 18-51-18005) and by the RUDN University Program (5-100).
REFERENCES [1] A. A. Samarskii.
Theory of Difference Schemes . Marcel Dekker, New York, 2001.[2] S. H. Christiansen, H. Z. Munthe-Kaas and B. Owren.
Topics in structure-preserving discretization . Acta Numerica, 11, 1-119, 2011.[3] B. Koren, R. Abgral, P. Bochev, J. Frank and B. Perot (eds.)
Physics - compatiblenumerical methods . J. Comput. Phys., 257, Part B, 1039–1526, 2014.[4] V. P. Gerdt and D. Robertz.
Consistency of Finite Difference Approximations forLinear PDE Systems and its Algorithmic Verification . in: S. Watt (ed.). Proceedingsof ISSAC 2010, pp. 53–59. Association for Computing Machinery, 2010.[5] V. P. Gerdt.
Consistency Analysis of Finite Difference Approximations to PDE Sys-tems . Mathematical Modelling in Computational Physics / MMCP 2011, LNCS7125, pp. 28–42. Springer, Berlin, 2012. arXiv:math.AP/1107.4269[6] P. Amodio, Yu.A. Blinkov, V. P. Gerdt and R. La Scala.
On consistency of finitedifference approximations to the Navier–Stokes Equations . Computer Algebra inScientific Computing / CASC 2013, LNCS 8136, Springer, Cham, 2013, pp. 46–60.[7] P. Amodio, Yu. A. Blinkov, V. P. Gerdt and R. La Scala.
Algebraic constructionand numerical behavior of a new s-consistent difference scheme for the 2D Navier–Stokes equations . Appl. Math. and Comput., 314, 408–421, 2017.[8] A. Levin.
Difference Algebra . Algebra and Applications, 8, Springer, 2008.[9] V. Gerdt and R. La Scala.
Noetherian quotients of the algebra of partial differencepolynomials and Gröbner bases of symmetric ideals. J. Algebra, 423, 2015, 1233–1261. [10] Xiao-Shan Gao, Zhang Huang and Chun-Ming Yuan,
Binomial difference ideals ,J. Symb. Comput., 80, 2017, 665–706[11] V. P. Gerdt, M. Lange-Hegermann and D. Robertz.
The Maple package TDDS forcomputing Thomas decompositions of systems of nonlinear PDEs . Comput. Phys.Commun., 234, 202–215, 2019. arXiv:physics.comp-ph/1801.09942[12] T. Bächler, V. Gerdt, M. Lange-Hegermann and D. Robertz.
Algorithmic Thomasdecomposition of algebraic and differential systems . J. Symb. Comput., 47(10),2012, 1233–1266. [13] D. Robertz.
Formal Algorithmic Elimination for PDEs , volume 2121 of
LectureNotes in Mathematics . Springer, Cham, 2014.[14] J. M. Thomas.
Differential Systems . AMS Colloquium Publications XXI, 1937;
Sys-tems and Roots . The Wylliam Byrd Press, Rychmond, Virginia, 1962.[15] V. P. Gerdt.
Involutive Algorithms for Computing Gröbner Bases . ComputationalCommutative and Non-Commutative Algebraic Geometry. IOS Press, Amster-dam, 2005, pp. 199–225. arXiv:math.AC/0501111.[16] W. M. Seiler.
Involution: The Formal Theory of Differential Equations and its Ap-plications in Computer Algebra . Algorithms and Computation in Mathematics,24, Springer, 2010.[17] Ch. Riquier.
Les systèmes d’équations aux dérivées partielles . Gauthiers-Villars,Paris, 1910.[18] M. Janet.
Leçons sur les systèmes d’équations aux dérivées partielles . Cahiers Sci-entifiques, IV. Gauthier-Villars, Paris, 1929.[19] Wu Wen-tsun.
On the Construction of Groebner Basis of a Polynomial Ideal Basedon Ruquier-Janet Theory , 5–22, 1990.[20] J. C. Strikwerda. Finite Difference Schemes and Partial Differential Equations , 2ndEdition. SIAM, Philadelphia, 2004.[21] V. P. Gerdt and D. Robertz.
Computation of difference Gröbner bases . Comput. Sc.J. Moldova, 20, 2(59), 203–226, 2012. Package
LDA is freely available on the webpage http://134.130.169.213/Janet/ [22] Dongming Wang.
Elimination Methods . Springer, Wien, 2000.[23] D. Rempfer.
On Boundary Conditions for Incompressible Navier-Stokes Problems .Appl. Mech. Rev., 59, 107–125, 2006.[24] Yu. A. Blinkov, C. F. Cid, V. P. Gerdt, W. Plesken, D. Robertz.
The MAPLE Pack-age Janet: II. Linear Partial Differential Equations . In: V. G. Ganzha, E. W. Mayr,E. V. Vorozhtsov (eds.) CASC 2003. Proc. 6th Int. Workshop on Computer Alge-bra in Scientific Computing, pp. 41–54. TU München (2003). Package
Janet isfreely available on the web pageisfreely available on the web page