Algorithmic approach to strong consistency analysis of finite difference approximations to PDE systems

TL;DR

This paper presents an algorithmic method for analyzing the strong consistency of finite difference schemes for polynomially nonlinear PDE systems, using differential Thomas decomposition and its difference analogue.

Contribution

It introduces a novel algorithm combining differential and difference Thomas decompositions to verify s-consistency of finite difference approximations for nonlinear PDEs.

Findings

Verified s-consistency for Navier-Stokes equations

Developed difference Thomas decomposition algorithm

Provided a systematic approach for PDE approximation analysis

Abstract

For a wide class of polynomially nonlinear systems of partial differential equations we suggest an algorithmic approach to the s(trong)-consistency analysis of their finite difference approximations on Cartesian grids. First we apply the differential Thomas decomposition to the input system, resulting in a partition of the solution set. We consider the output simple subsystem that contains a solution of interest. Then, for this subsystem, we suggest an algorithm for verification of s-consistency for its finite difference approximation. For this purpose we develop a difference analogue of the differential Thomas decomposition, both of which jointly allow to verify the s-consistency of the approximation. As an application of our approach, we show how to produce s-consistent difference approximations to the incompressible Navier-Stokes equations including the pressure Poisson equation.

Equations

$$f_1=\cdots =f_s=0,F:=\{f_1,\dots ,f_s\}\subset \mathcal{R},s\in {\mathbb{Z}}_{\ge 1},$$

$$\{(x_1+k_{1}h,\dots ,x_n+k_{n}h)\mid k_1,\dots ,k_n\in \mathbb{Z}\},$$

$${\sigma }_i^{\pm 1}{\tilde{u}}_{k_1,\dots ,k_i,\dots ,k_n}^{(\alpha )}={\tilde{u}}_{k_1,\dots ,k_i\pm 1,\dots ,k_n}^{(\alpha )},\alpha \in \{1,\dots ,m\}.$$

$${\partial }_j{u}^{(\alpha )}(\mathbf{x})=\frac{u_{k_1,\dots ,k_j+1,\dots ,k_n}^{(\alpha )}-u_{k_1,\dots ,k_j-1,\dots ,k_n}^{(\alpha )}}{2h}+\mathcal{O}(h^2).$$

$${\tilde{f}}_1=\cdots ={\tilde{f}}_s=0,\tilde{F}:=\{{\tilde{f}}_1,\dots ,{\tilde{f}}_s\}\subset \tilde{\mathcal{R}},$$

$$(\forall j\in \{1,\dots ,s\})[{\tilde{f}}_j\underset{h\to 0}{\to }{f}_j].$$

$$\tilde{f}(\tilde{\mathbf{u}})=h^{d}f(\mathbf{u})+\mathcal{O}(h^{d+1}),d\in {\mathbb{Z}}_{\ge 0},$$

$$(\forall \tilde{f}\in [[\tilde{F}]])(\exists f\in [[F]])[\tilde{f}⊳f].$$

$$\frac{\tilde{f}(\tilde{\mathbf{u}})}{h^d}\underset{h\to 0}{\to }f(\mathbf{u}),$$

$$(\forall \tilde{f}\in \tilde{\mathcal{I}})(\exists f\in \mathcal{I})[\tilde{f}⊳f].$$

$$\tilde{f}=\sum _{\tilde{g}\in {\tilde{G}}_1\subseteq \tilde{G}}\sum _{\mu }{a}_{\tilde{g},\mu }{\sigma }^{\mu }\tilde{g},a_{\tilde{g},\mu }\in \tilde{\mathcal{R}},\mathrm{lm}(a_{\tilde{g},\mu }{\sigma }^{\mu }\tilde{g})⪯\mathrm{lm}(\tilde{f}).$$

$$\mu :=({\mu }_1,\dots ,{\mu }_n)\in {\mathbb{Z}}_{\ge 0}^n,{\sigma }^{\mu }:={\sigma }_1^{{\mu }_1}\cdots {\sigma }_n^{{\mu }_n}.$$

$$f:=\sum _{g\in G_1}\sum _{\nu }{b}_{g,\nu }{\partial }^{\nu }g,b_{g,\nu }\in \mathcal{R},$$

$$\tilde{q}:=({\theta }_1\tilde{p}{)}^{k_1}\cdots ({\theta }_r\tilde{p}{)}^{k_r}\in [\tilde{F}],k_1,\dots ,k_r\in {\mathbb{Z}}_{\ge 0}.$$

$$(\forall \tilde{g}\in \tilde{G})(\exists g\in [[F]])[\tilde{g}⊳g].$$

$$i_k=\max \{j_k\mid y_1^{j_1}\cdots y_n^{j_n}\in G\text{mbox}with{j}_1=i_1,\dots ,j_{k-1}=i_{k-1}\}.$$

$$⟨G⟩:=\bigcup _{m\in G}\mathrm{Mon}(R)m=\underset{m\in G}{⨄}\mathrm{Mon}(\mu (m,G))m.$$

$$p_1=0,\dots ,p_s=0,p_{s+1}\ne 0,\dots ,p_{s+t}\ne 0(s,t\in {\mathbb{Z}}_{\ge 0})$$

$${\mathrm{Sol}}_{\overline{K}}(S):=\{a\in {\overline{K}}^n\mid p_i(a)=0,p_{s+j}(a)\ne 0,i=1,...,s,j=1,...,t\}.$$

$${\pi }_i:{\overline{K}}^n\to {\overline{K}}^{n-i}:(a_1,...,a_n)\mapsto (a_{i+1},...,a_n),i=1,2,...,n-1.$$

$$f_1=0,\dots ,f_s=0,f_{s+1}\ne 0,\dots ,f_{s+t}\ne 0(s,t\in {\mathbb{Z}}_{\ge 0})$$

$$\begin{array}[]{c}\operatorname{Sol}_{\Omega}(S):=\{\,a=(a_{1},\ldots,a_{m})\mid a_{k}\colon\Omega\to\mathbb{C}\mbox{ analytic},\,k=1,\ldots,m,\\[1.99997pt] f_{i}(a)=0,\,f_{s+j}(a)\neq 0,\,i=1,\ldots,s,\,j=1,\ldots,t\,\}\,.\end{array}$$

$$\mathrm{Mon}(\Delta )u:=\{u_J^{(\alpha )}\mid 1\le \alpha \le m,J\in ({\mathbb{Z}}_{\ge 0}{)}^n\}$$

$$u_J^{(\alpha )}{\succ }_{\mathrm{TOP},\mathrm{lex}}{u}_K^{(\beta )}:\iff J{\succ }_{\mathrm{lex}}K\text{mbox}or(J=K\text{mbox}and\alpha <\beta )$$

$$u_J^{(\alpha )}{\succ }_{\mathrm{POT},\mathrm{lex}}{u}_K^{(\beta )}:\iff \alpha <\beta \text{mbox}or(\alpha =\beta \text{mbox}andJ{\succ }_{\mathrm{lex}}K),$$

$$G_{\alpha }:=\{\theta \in \mathrm{Mon}(\Delta )\mid \theta u^{\alpha }\in \{\mathrm{ld}(f_1),\dots ,\mathrm{ld}(f_s)\}\}.$$

$$\mathrm{NF}(\partial f_i,T,\succ )=0\text{mbox}forall\partial \in {\overline{\mu }}_i=\Delta \setminus {\mu }_i,i=1,\dots ,s.$$

$$E:q^{\infty }:=\{f\in \mathcal{R}\mid q^{r}f\in E\text{mbox}forsomer\in {\mathbb{Z}}_{\ge 0}\}$$

$${\tilde{f}}_1=0,\dots ,{\tilde{f}}_s=0,{\tilde{f}}_{s+1}\ne 0,\dots ,{\tilde{f}}_{s+t}\ne 0(s,t\in {\mathbb{Z}}_{\ge 0})$$

$$E:Q:=\{\tilde{f}\in \tilde{\mathcal{R}}\mid q\tilde{f}\in E\text{mbox}forsomeq\in Q\}.$$

Full text

Partial differential equations, Finite difference approximations, Consistency, Thomas decomposition

††isbn: ††copyright: acmcopyright††conference: International Symposium on Symbolic and Algebraic Computation; July 15-18, 2019; Beijing, China††journalyear: 2019††price: 15.00††ccs: Computing methodologies Algebraic algorithms††ccs: Mathematics of computing Partial differential equations††ccs: Mathematics of computing Nonlinear equations††ccs: Mathematics of computing Discretization

1. Introduction

Except very special cases, partial differential equations (PDE) admit numerical integration only. Historically first and one of the most-used numerical methods is finite difference method (Samarskii'01, ) based on approximation of PDE by difference equations defined on a chosen solution grid. To construct a numerical solution, the obtained finite difference approximation (FDA) to PDE is augmented with an appropriate discretization of initial or/and boundary condition(s) providing uniqueness of solution. As this takes place, the quality of numerical​ solution to PDE​ is​ determined by the quality of its​ FDA.

Any reasonable discretization must provide the convergence of a numerical solution to a solution of PDE in the limit when the grid spacings tend to zero. However, except for a very limited class of problems, convergence cannot be directly established. In practice, for a given FDA, its consistency and stability are analyzed as the necessary conditions for convergence. Consistency implies reduction of the FDA to the original PDE when the grid spacings tend to zero and stability provides boundedness of the error in the solution under 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 structure (Christiansen'11, ) of the PDE. Such mimetic or algebraic structure preserving FDA are more likely to produce highly accurate and stable numerical results (cf. (JCP'14, )). In (GR'10, ; G'12, ), for polynomially nonlinear PDE systems and regular solution grids, we introduced the novel concept of strong consistency, or s-consistency, which strengthens the concept of consistency and means that any element of the perfect difference ideal generated by the polynomials in FDA approximates an element in the radical differential ideal generated by the polynomials in PDE. In the subsequent papers (ABGLS'13, ; ABGLS'17, ), by computational experiments with two-dimensional incompressible Navier-Stokes equations, it was shown that s-consistent FDA have much better numerical behavior than FDA which are not s-consistent.

For linear PDE one can algorithmically verify (GR'10, ) s-consistency of their FDA. In the nonlinear case such verification (G'12, ) required computation of a difference Gröbner basis for FDA. Since difference polynomial rings (Levin'08, ) are non-Noetherian, the difference Gröbner basis algorithms (G'12, ; GLS'15, ) do not terminate in general. In comparison to differential algebra, fewer computational results have been obtained in difference algebra. A decomposition technique was developed only for binomial perfect difference ideals (BinomialDifference, ). More generally, in the present paper, a difference analogue of the differential Thomas decomposition (Thomas'37-62, ; BGLHR'12, ; Robertz6, ; GLHR'18, ) is obtained (see Section 6), which provides an algorithmic tool for s-consistency analysis of FDA to simple PDE subsystems on Cartesian grids (see Section 7). In particular, given an FDA to the momentum and continuity equations in the Navier-Stokes PDE system for incompressible flow, our approach derives an s-consistent approximation containing the pressure Poisson equation (see Section 9).

Completion to involution is the cornerstone of the differential Thomas decomposition (Thomas'37-62, ; BGLHR'12, ; Robertz6, ; GLHR'18, ). The underlying completion algorithm (G'05, ) is based on the theory of Janet division and Janet bases (G'05, ; Seiler'10, ; Robertz6, ) which stemmed from the Riquier-Janet theory (Riquier'10, ; Janet'29, ) of orthonomic PDE. Joseph M. Thomas (Thomas'37-62, ) generalized the Riquier-Janet theory to non-orthonomic polynomially nonlinear PDE and showed how to decompose them into the triangular subsystems with disjoint solution sets. Janet bases are Gröbner ones with additional structure, and Wu Wen-tsun was the first who showed (Wu'90, ) that the Riquier-Janet theory can be used for algorithmic construction of algebraic Gröbner bases. We dedicate this paper to commemoration of his Centennial Birthday.

2. Consistency

In the given paper we consider PDE systems of the form

[TABLE]

where ${\mathcal{R}}:={\mathcal{K}}\{\mathbf{u}\}$ is the ring of polynomials in the dependent variables $\mathbf{u}:=\{u^{(1)},\ldots,u^{(m)}\}$ and their partial derivatives obtained from the operator power products in $\{\partial_{1},\ldots,\partial_{n}\}$ $(\partial_{j}=\partial_{x_{j}})$. We shall assume that coefficients of the polynomials are rational functions in $\mathbf{a}:=\{a_{1},\ldots,a_{l}\}$, finitely many parameters (constants), over $\mathbb{Q}$, i.e. ${\mathcal{K}}:={\mathbb{Q}}(\mathbf{a})$. One can also extend the last field to ${\mathbb{Q}}(\mathbf{a},\mathbf{x})$, where $\mathbf{x}:=\{x_{1},\ldots,x_{n}\}$ is the set of independent variables. In this case we shall assume that coefficients of the differential polynomials in $F$ do not vanish in the grid points defined in (2) below.

To approximate (1) by a difference system we define a Cartesian computational grid (mesh) with spacing $0<h\in\mathbb{R}$ and fixed $\mathbf{x}$ by

[TABLE]

If the actual solution to (1) is $\mathbf{{u}}(\mathbf{x})$, then its approximation at the grid nodes will be denoted by $\mathbf{\tilde{u}}_{k_{1},\ldots,k_{n}}\approx\mathbf{u}(x_{1}+k_{1}h,\ldots,x_{n}+k_{n}h)$.

Let $\tilde{\mathcal{K}}:=\mathbb{Q}\,(\mathbf{a},h)$ and $\tilde{\mathcal{R}}$ be the difference polynomial ring over $\tilde{\mathcal{K}}$, where $\tilde{\mathcal{K}}$ is a difference field of constants (Levin'08, ) with differences $\Sigma:=\{\sigma_{1},\ldots,\sigma_{n}\}$ acting on a grid function $\tilde{u}^{(\alpha)}_{k_{1},\ldots,k_{n}}$ as the shift operators

[TABLE]

The elements in $\tilde{\mathcal{R}}$ are polynomials in the dependent variables $\tilde{u}^{(\alpha)}$ ($\alpha=1,...,m$) defined on the grid and in their shifts $\sigma_{1}^{i_{1}}...\sigma_{n}^{i_{n}}\tilde{u}^{(\alpha)}$ $(i_{j}\in\mathbb{Z})$. However, to provide termination of the decomposition algorithm of Sect. 6, we shall consider difference polynomials with non-negative shifts only. We denote by $\operatorname{Mon}(\Sigma)$ the set of monomials in $\sigma_{1}$, …, $\sigma_{n}$. The coefficients of the polynomials are in $\tilde{\mathcal{K}}$.

The standard method to obtain FDA of such type to the differential system (1) is replacement of the partial derivatives occurring in (1) by finite differences and application of appropriate power product of the forward-shift operators in (3) to eliminate negative shifts in indices which may come out of expressions like

[TABLE]

Furthermore, the difference system

[TABLE]

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

[TABLE]

This is a universally adopted notion of consistency for a finite difference discretization of PDE system (1) (cf. (Str'04, ), Ch.7) and means that Eq. (4) reduces to Eq. (1) when the mesh step $h$ tends to zero.

Definition 2.2.

(GR'10, )  We say that a difference equation $\tilde{f}(\mathbf{\tilde{u}})=0$, $\tilde{f}\in\tilde{\mathcal{R}}$, implies the differential equation $f(\mathbf{\mathbf{u}})=0$, $f\in{\mathcal{R}}$, and write $\tilde{f}\rhd f$, if the Taylor expansion of $\tilde{f}$ about the grid point $\mathbf{x}$, after clearing denominators containing $h$, yields

[TABLE]

and $\mathcal{O}\,(\,h^{d+1}\,)$ denotes terms whose degree in $h$ is at least $d+1$.

Remark 2.3.

Given $\tilde{f}(\tilde{\mathbf{u}})$, computation of $f(\mathbf{u})$ is straightforward and has been implemented as routine ContinuousLimit in the Maple package LDA (GR'12, ; GLS'15, ) (Linear Difference Algebra).

Definition 2.4.

(G'12, ) FDA (4) to PDE system (1) is strongly consistent or s-consistent if

[TABLE]

Here $\llbracket\tilde{F}\rrbracket$ and $\llbracket F\rrbracket$ denote the perfect difference ideal generated by $\tilde{F}$ in $\tilde{\mathcal{R}}$ and the radical differential ideal generated by $F$ in ${\mathcal{R}}$.

Remark 2.5.

It is clear that if condition (5) holds, then

[TABLE]

that is, $\tilde{f}(\mathbf{\tilde{u}})/h^{d}$ approximates $f(\mathbf{u})$. Accordingly, condition (6) means that, after clearing denominators, each element of $\llbracket\tilde{F}\rrbracket$ approximates an element of $\llbracket F\rrbracket$ in the sense of (7).

Lemma 2.6.

Let ${\mathcal{I}}=[{F}]$ be a differential ideal of ${\mathcal{R}}$ and $\tilde{\mathcal{I}}=[\tilde{F}]$ a difference ideal of $\tilde{\mathcal{R}}$ such that

[TABLE]

Then for the perfect closure $\llbracket\tilde{\mathcal{I}}\rrbracket$ of $\tilde{\mathcal{I}}$ in $\tilde{\mathcal{R}}$ the condition (6) holds.

Proof.

Let $\tilde{G}$ be a (possibly infinite) reduced Gröbner basis of $\tilde{\mathcal{I}}$ for an admissible monomial ordering $\succ$ (cf. (G'12, )). Then

[TABLE]

Here $\tilde{f}\in\tilde{\mathcal{I}}$, $\tilde{G}_{1}$ is a finite subset of $\tilde{G}$, $\operatorname{lm}$ denotes the leading monomial of its argument, and we use the multi-index notation

[TABLE]

In the continuous limit $\tilde{f}$ implies the differential polynomial

[TABLE]

where $\tilde{G}_{1}\rhd G_{1}$. Therefore, $\tilde{f}\rhd f\in[F]\subseteq\llbracket F\rrbracket$.

Let now $\tilde{p}\in\llbracket\tilde{F}\rrbracket\setminus[\tilde{F}]$ and $\theta_{1},\ldots,\theta_{r}\in\operatorname{Mon}(\Sigma)$ be such that

[TABLE]

From Eq. (8) it follows that $\tilde{q}\,\rhd\,q=p^{k_{1}+\cdots+k_{r}}$ where $\tilde{p}\,\rhd\,p$. Hence, $p\in\llbracket F\rrbracket$. The perfect ideal $\llbracket\tilde{F}\rrbracket$ can be constructed (Levin'08, ) from $[\tilde{F}]$ by the procedure in the form called shuffling and based on enlargement of the generator set $\tilde{F}$ with all polynomials $\tilde{p}$ occurring in $[\tilde{F}]$ in the form of Eq. (8) and on repetition of such enlargement. It is clear that each such enlargement of the intermediate ideals yields in the continuous limit a subset of $\llbracket F\rrbracket$. ∎

The criterion of s-consistency is given by the following theorem.

Theorem 2.7.

(G'12, )* A difference approximation (4) to a differential system (1) is s-consistent if and only if a reduced Gröbner basis $\tilde{G}\subset\tilde{\mathcal{R}}$ of the difference ideal $[\tilde{F}]\subset\tilde{\mathcal{R}}$ generated by $\tilde{F}$ satisfies*

[TABLE]

3. Janet division

We recall the concept of Janet division. For details we refer to, e.g., (Robertz6, , Subsect. 2.1.1), (G'05, ), (Seiler'10, , Ch. 3).

Let $K$ be a field and $R:=K[y_{1},\ldots,y_{n}]$ the commutative polynomial algebra over $K$ with indeterminates $y_{1}$, …, $y_{n}$. We denote by $\operatorname{Mon}(R)$ the set of monomials in $y_{1}$, …, $y_{n}$ and for a subset $\mu\subseteq\{y_{1},...,y_{n}\}$ we define $\operatorname{Mon}(\mu)$ to be the subset of $\operatorname{Mon}(R)$ consisting of the monomials involving only indeterminates from $\mu$.

If a term ordering on $R$ is fixed and $I$ is an ideal of $R$, then the set of leading monomials of non-zero polynomials in $I$ are known to form a set with the following property:

Definition 3.1.

A set $M\subseteq\operatorname{Mon}(R)$ is said to be $\operatorname{Mon}(R)$-multiple-closed if we have $rm\in M$ for all $m\in M$ and all $r\in\operatorname{Mon}(R)$.

The smallest $\operatorname{Mon}(R)$-multiple-closed set in $\operatorname{Mon}(R)$ containing a given set $G\subseteq\operatorname{Mon}(R)$ is denoted by $\langle\,G\,\rangle$. It is well known that every $\operatorname{Mon}(R)$-multiple-closed set in $\operatorname{Mon}(R)$ is finitely generated in that sense and that it has a unique minimal generating set.

We adopt Janet’s approach (Janet'29, ) of partitioning a $\operatorname{Mon}(R)$-multiple-closed set $M$ into finitely many subsets of the form $\operatorname{Mon}(\mu)\,m$, where $m\in M$ and $\mu=\mu(m,M)\subseteq\operatorname{Mon}(R)$ (referred to as Janet division).

Definition 3.2.

Let $G\subset\operatorname{Mon}(R)$ be finite and $m=y_{1}^{i_{1}}\cdots y_{n}^{i_{n}}\in G$. Then $y_{k}$ is said to be a multiplicative variable for $m$ if and only if

[TABLE]

This yields a partition $\{y_{1},\ldots,y_{n}\}=\mu(m,G)\uplus\overline{\mu}(m,G)$, where the elements of $\mu(m,G)$ (resp. $\overline{\mu}(m,G)$) are the multiplicative (resp. non-multiplicative) variables for $m$. The set $G$ is Janet complete if

[TABLE]

Proposition 3.3.

For every $\operatorname{Mon}(R)$-multiple-closed set $M$ there exists a finite Janet complete set $J\subset\operatorname{Mon}(R)$ such that $M=\langle\,J\,\rangle$.

If $G\subset\operatorname{Mon}(R)$ is finite, we call the minimal Janet complete set $J\supset G$ such that $\langle\,J\,\rangle=\langle\,G\,\rangle$ the Janet completion of $G$. It is obtained algorithmically by adding certain multiples of elements of $G$ to $G$ (which also proves Proposition 3.3), cf., e.g., (Robertz6, , Algorithm 2.1.6).

4. Simple Algebraic Systems

Fundamental for both the differential Thomas decomposition (recalled in Sect. 5) as well as its difference analogue to be introduced in Sect. 6 is the Thomas decomposition of an algebraic system $S$

[TABLE]

where $p_{1}$, …, $p_{s+t}\in R:=K[z_{1},\ldots,z_{n}]$. Here $K$ is a field of characteristic zero with algebraic closure $\overline{K}$, and $R$ is the commutative polynomial algebra over $K$ with indeterminates $z_{1}$, …, $z_{n}$. The solution set of the algebraic system $S$ in (9) is defined to be

[TABLE]

Assuming the indeterminates are ordered as in $z_{1}\succ z_{2}\succ\ldots\succ z_{n}$, a sequence of projections from $\overline{K}^{n}$ is defined correspondingly by

[TABLE]

For each $p\in R\setminus K$, this ordering defines the greatest indeterminate $\operatorname{ld}(p)$ occurring in $p$, referred to as leader, the coefficient $\operatorname{init}(p)$ of the highest power of $\operatorname{ld}(p)$ in $p$, called initial, and the discriminant $\operatorname{disc}(p):=(-1)^{d(d-1)/2}\operatorname{res}(p,\partial p/\partial\operatorname{ld}(p),\operatorname{ld}(p))\,/\,\operatorname{init}(p)$, where $d$ is the degree of $p$ in $\operatorname{ld}(p)$ and where $\operatorname{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. (1)

None of $p_{1}$, …, $p_{s}$, $p_{s+1}$, …, $p_{s+t}$ is constant. 2. (2)

The leaders of $p_{1}$, …, $p_{s}$, $p_{s+1}$, …, $p_{s+t}$ are pairwise distinct. 3. (3)

For every $r\in\{p_{1},\ldots,p_{s},p_{s+1},\ldots,p_{s+t}\}$, if $\operatorname{ld}(r)=z_{k}$, then the equation $\operatorname{init}(r)=0$ has no solution in $\pi_{k}(\operatorname{Sol}_{\overline{K}}(S))$. 4. (4)

For every $r\in\{p_{1},\ldots,p_{s},p_{s+1},\ldots,p_{s+t}\}$, if $\operatorname{ld}(r)=z_{k}$, then the equation $\operatorname{disc}(r)=0$ has no solution in $\pi_{k}(\operatorname{Sol}_{\overline{K}}(S))$.

(In (3) and (4), we have $\operatorname{init}(r)$, $\operatorname{disc}(r)\in K[z_{k+1},\ldots,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) is a finite collection of simple algebraic systems $S_{1}$, …, $S_{r}$ such that $\operatorname{Sol}_{\overline{K}}(S)\!=\!\operatorname{Sol}_{\overline{K}}(S_{1})\uplus...\uplus\operatorname{Sol}_{\overline{K}}(S_{r})$. It can be computed by an algorithm combining Euclidean pseudo-reduction and case distinctions. For details we refer to (BGLHR'12, ), (Robertz6, , Subsect. 2.2.1), (Wang'00, , Sect. 3.3).

5. Differential Thomas Decomposition

A system​ of polynomial partial differential equations and inequations

[TABLE]

is given by elements $f_{1}$, …, $f_{s+t}$ of the differential polynomial ring ${\mathcal{R}}$ in $u^{(1)}$, …, $u^{(m)}$ with commuting derivations $\Delta:=\{\partial_{1},\ldots,\partial_{n}\}$. For $\alpha\in\{1,\ldots,m\}$, $J\in(\mathbb{Z}_{\geq 0})^{n}$ we identify $u^{(\alpha)}_{J}$ and $\partial_{1}^{J_{1}}\cdots\partial_{n}^{J_{n}}u^{(\alpha)}$. Let $\Omega\subseteq\mathbb{R}^{n}$ be open and connected. The solution set of $S$ on $\Omega$ is

[TABLE]

Definition 5.1.

A ranking $\succ$ on ${\mathcal{R}}$ is a total ordering on the set

[TABLE]

such that for all $j\in\{1,\ldots,n\}$, $\alpha$, $\beta\in\{1,\ldots,m\}$, $J$, $K\in(\mathbb{Z}_{\geq 0})^{n}$ we have $\partial_{j}u^{(\alpha)}\succ u^{(\alpha)}$ and, if $u^{(\alpha)}_{J}\succ u^{(\beta)}_{K}$, then $\partial_{j}u^{(\alpha)}_{J}\succ\partial_{j}u^{(\beta)}_{K}$. A ranking $\succ$ is orderly if for all $\alpha$, $\beta\in\{1,\ldots,m\}$, $J$, $K\in(\mathbb{Z}_{\geq 0})^{n}$, $J_{1}+\cdots+J_{n}>K_{1}+\cdots+K_{n}$ implies $u^{(\alpha)}_{J}\succ u^{(\beta)}_{K}$.

Example 5.2.

Rankings $\succ_{{\rm TOP},{\rm lex}}$ and $\succ_{{\rm POT},{\rm lex}}$ on ${\mathcal{R}}$ are given by

[TABLE]

and

[TABLE]

respectively, where $\succ_{{\rm lex}}$ compares multi-indices lexicographically.

If a ranking $\succ$ on ${\mathcal{R}}$ is fixed, then for each $f\in{\mathcal{R}}\setminus{\mathcal{K}}$ the leader, initial and discriminant of $f$ are defined as in Section 4. Moreover, $\operatorname{sep}(f):=\partial f/\partial\operatorname{ld}(f)$ is called the separant of $f$.

Janet division associates (with respect to a total ordering of $\Delta$) to each $f_{i}=0$ with $\operatorname{ld}(f_{i})=\theta_{i}u^{(\alpha)}$ the set $\mu_{i}:=\mu(\theta_{i},G_{\alpha})\subseteq\Delta$ (resp. $\overline{\mu}_{i}:=\Delta\setminus\mu_{i}$) of admissible (resp. non-admissible) derivations, where

[TABLE]

We call $\{f_{1}\!=\!0,...,f_{s}\!=\!0\}$ or $T:=\{(f_{1},\mu_{1}),...,(f_{s},\mu_{s})\}$ Janet complete if each $G_{\alpha}$ equals its Janet completion, $\alpha=1$, …, $m$. Let $r\in{\mathcal{R}}$. If some $v\in\operatorname{Mon}(\Delta)u$ occurs in $r$ for which there exists $(f,\mu)\in T$ such that $v=\theta\operatorname{ld}(f)$ for some $\theta\in\operatorname{Mon}(\mu)$ and $\deg_{v}(r)\geq\deg_{v}(\theta f)$, then $r$ is Janet reducible modulo $T$. In this case, $(f,\mu)$ 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 $\operatorname{NF}(r,T,\succ)$, a Janet reduced differential polynomial, as explained in (Robertz6, , Algorithm 2.2.45).

Definition 5.3.

Let $T=\{\,(f_{1},\mu_{1}),\ldots,(f_{s},\mu_{s})\,\}$ be Janet complete. Then $\{\,f_{1}=0,\ldots,f_{s}=0\,\}$ or $T$ is said to be passive, if

[TABLE]

Definition 5.4.

Let a ranking $\succ$ on ${\mathcal{R}}$ and a total ordering on $\Delta$ be fixed. A differential system $S$ as in (10) is said to be simple if the following three conditions hold.

1. (1)

$S$ is simple as an algebraic system (in the finitely many indeterminates occurring in it, ordered by the ranking $\succ$). 2. (2)

$\{\,f_{1}=0,\,\ldots,\,f_{s}=0\,\}$ is passive. 3. (3)

The left hand sides $f_{s+1}$, …, $f_{s+t}$ are Janet reduced modulo the passive differential system $\{\,f_{1}=0,\,\ldots,\,f_{s}=0\,\}$.

Proposition 5.5 ((Robertz6, ), Prop. 2.2.50).

Let $S$ be a simple differential system, defined over ${\mathcal{R}}$, as in (10). Let $E$ be the differential ideal of ${\mathcal{R}}$ which is generated by $f_{1}$, …, $f_{s}$ and let $q$ be the product of the initials and separants of all $f_{1}$, …, $f_{s}$. Then the differential ideal

[TABLE]

is radical. Given $f\in{\mathcal{R}}$, we have $f\in E:q^{\infty}$ if and only if the Janet normal form of $f$ modulo $f_{1}$, …, $f_{s}$ is zero.

Definition 5.6.

A Thomas decomposition of a differential system $S$ as in (10) (with respect to $\succ$) is a finite collection of simple differential systems $S_{1}$, …, $S_{r}$ such that $\operatorname{Sol}_{\Omega}(S)=\operatorname{Sol}_{\Omega}(S_{1})\uplus...\uplus\operatorname{Sol}_{\Omega}(S_{r})$.

For any differential system $S$ as in (10) and any ranking $\succ$ on ${\mathcal{R}}$ a Thomas decomposition of $S$ can be computed algorithmically. For more details we refer to, e.g., (BGLHR'12, ), (Robertz6, , Subsection 2.2.2), (GLHR'18, ).

6. Decomposition of difference systems

A system $\tilde{S}$ of polynomial partial difference equations and inequations

[TABLE]

is given by elements $\tilde{f}_{1}$, …, $\tilde{f}_{s+t}$ of the difference polynomial ring $\tilde{\mathcal{R}}$ in $\tilde{u}^{(1)}$, …, $\tilde{u}^{(m)}$ with commuting automorphisms $\Sigma=\{\sigma_{1},\ldots,\sigma_{n}\}$. For $\alpha\in\{1,\ldots,m\}$, $J\in(\mathbb{Z}_{\geq 0})^{n}$ we identify $\tilde{u}^{(\alpha)}_{J}$ and $\sigma_{1}^{J_{1}}\cdots\sigma_{n}^{J_{n}}\tilde{u}^{(\alpha)}$. We denote by $\tilde{S}^{=}$ (resp. $\tilde{S}^{\neq}$) the set $\{\tilde{f}_{1},...,\tilde{f}_{s}\}$ (resp. $\{\tilde{f}_{s+1},...,\tilde{f}_{s+t}\}$).

A ranking on $\tilde{\mathcal{R}}$ is defined in the same way as in Definition 5.1 by replacing the action of $\partial_{i}$ by the action of $\sigma_{i}$ and $\Delta$ by $\Sigma$.

For a subset $L$ of $\tilde{\mathcal{R}}$ we denote by $[L]$ the difference ideal of $\tilde{\mathcal{R}}$ generated by $L$. Let $E$ be a difference ideal of $\tilde{\mathcal{R}}$ and $\emptyset\neq Q\subseteq\tilde{\mathcal{R}}$ be multiplicatively closed and closed under $\sigma_{1}$, …, $\sigma_{n}$. Then define

[TABLE]

Moreover, for $U\subseteq\operatorname{Mon}(\Sigma)\,\tilde{u}$ and $v\in\operatorname{Mon}(\Sigma)\,\tilde{u}$ we define

[TABLE]

The first algorithm to be introduced performs an auto-reduction of a finite set of difference polynomials.

Since leaders are dealt with in decreasing order with respect to $\succ$, 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 $\Sigma$) to each $\tilde{f}_{i}=0$ with $\operatorname{ld}(\tilde{f}_{i})=\theta_{i}\tilde{u}^{(\alpha)}$ the set $\mu_{i}:=\mu(\theta_{i},\tilde{G}_{\alpha})\subseteq\Sigma$ (resp. $\overline{\mu}_{i}:=\Sigma\setminus\mu_{i}$) of admissible (non-admissible) automorphisms, where

[TABLE]

We call $\{\tilde{f}_{1}\!=\!0,...,\tilde{f}_{s}\!=\!0\}$ or $T:=\{(\tilde{f}_{1},\mu_{1}),...,(\tilde{f}_{s},\mu_{s})\}$ Janet complete if each $\tilde{G}_{\alpha}$ equals its Janet completion, $\alpha=1$, …, $m$. Let $\tilde{r}\in\tilde{\mathcal{R}}$. If some $v\in\operatorname{Mon}(\Sigma)\tilde{u}$ occurs in $\tilde{r}$ for which there exists $(\tilde{f},\mu)\in T$ such that $v=\theta\operatorname{ld}(\tilde{f})$ for some $\theta\in\operatorname{Mon}(\mu)$ and $\deg_{v}(\tilde{r})\geq\deg_{v}(\theta\tilde{f})$, then $\tilde{r}$ is Janet reducible modulo $T$. In this case, $(\tilde{f},\mu)$ is called a Janet divisor of $\tilde{r}$. If $\tilde{r}$ is not Janet reducible modulo $T$, then $\tilde{r}$ is also said to be Janet reduced modulo $T$. Iterated pseudo-reductions of $\tilde{r}$ modulo $T$ yield its Janet normal form $\operatorname{NF}(\tilde{r},T,\succ)$, which is the Janet reduced difference polynomial $\tilde{r}^{\prime}$ returned by Algorithm 2.

Algorithm 2 terminates because each coefficient $\tilde{c}$ of $\tilde{r}^{\prime}$ is either constant or has a leader which is smaller than $\operatorname{ld}(\tilde{r}^{\prime})$ with respect to $\succ$, and a ranking $\succ$ does not allow infinitely decreasing chains. Correctness of the algorithm is clear.

Definition 6.1.

Let $T=\{\,(\tilde{f}_{1},\mu_{1}),\ldots,(\tilde{f}_{s},\mu_{s})\,\}$ be Janet complete. Then $\{\,\tilde{f}_{1}=0,\ldots,\tilde{f}_{s}=0\,\}$ or $T$ is said to be passive, if

[TABLE]

Definition 6.2.

Let a ranking $\succ$ on $\tilde{\mathcal{R}}$ and a total ordering on $\Sigma$ be fixed. A difference system $\tilde{S}$ as in (11) is said to be simple (resp., quasi-simple) if the following three conditions hold.

1. (1)

$\tilde{S}$ is simple (resp., quasi-simple) as an algebraic system (in the finitely many occurring indeterminates, ordered by $\succ$). 2. (2)

$\{\,\tilde{f}_{1}=0,\,\ldots,\,\tilde{f}_{s}=0\,\}$ is passive. 3. (3)

The left hand sides $\tilde{f}_{s+1}$, …, $\tilde{f}_{s+t}$ are Janet reduced modulo the passive difference system $\{\,\tilde{f}_{1}=0,\,\ldots,\,\tilde{f}_{s}=0\,\}$.

Proposition 6.3.

Let $\tilde{S}$ be a quasi-simple difference system over $\tilde{\mathcal{R}}$ as in (11). Let $E$ be the difference ideal of $\tilde{\mathcal{R}}$ generated by $\tilde{f}_{1}$, …, $\tilde{f}_{s}$ and let $Q$ be the smallest subset of $\tilde{\mathcal{R}}$ which is multiplicatively closed, closed under $\sigma_{1}$, …, $\sigma_{n}$ and contains the initials $q_{i}:=\operatorname{init}(\tilde{f}_{i})$ for all $i=1$, …, $s$. Then a difference polynomial $\tilde{f}\in\tilde{\mathcal{R}}$ is an element of

[TABLE]

if and only if the Janet normal form of $\tilde{f}$ modulo $\tilde{f}_{1}$, …, $\tilde{f}_{s}$ is zero.

Proof.

By definition of $E:Q$, every element $\tilde{f}\in\tilde{\mathcal{R}}$ for which Algorithm 2 yields Janet normal form zero is an element of $E:Q$.

Let $\tilde{f}\in E:Q$, $\tilde{f}\neq 0$. Then there exist $q\in Q$ and $k_{1}$, …, $k_{s}\in\mathbb{Z}_{\geq 0}$ and $c_{i,j}\in\tilde{\mathcal{R}}\setminus\{0\}$, $\alpha_{i,j}\in\operatorname{Mon}(\Sigma)$, $j=1$, …, $k_{i}$, $i=1$, …, $s$, such that

[TABLE]

Among all pairs $(i,j)$ for which $\alpha_{i,j}$ involves a non-admissible automorphism for $\tilde{f}_{i}=0$ let the pair $(i^{\star},j^{\star})$ be such that $\alpha_{i^{\star},j^{\star}}(\operatorname{ld}(\tilde{f}_{i^{\star}}))$ is maximal with respect to the ranking $\succ$. Let $\sigma$ be a non-admissible automorphism for $\tilde{f}_{i^{\star}}=0$ which divides the monomial $\alpha_{i^{\star},j^{\star}}$. Since $\{\tilde{f}_{1}=0,...,\tilde{f}_{s}=0\}$ is passive, there exist $b\in Q$, $l_{1}$, …, $l_{s}\in\mathbb{Z}_{\geq 0}$ and $d_{i,j}\in\tilde{\mathcal{R}}\setminus\{0\}$ and $\beta_{i,j}\in\operatorname{Mon}(\Sigma)$, $j=1$, …, $l_{i}$, $i=1$, …, $s$, such that

[TABLE]

where each $\beta_{i,j}$ involves only admissible automorphisms for $\tilde{f}_{i}=0$. Let $\gamma_{i^{\star},j^{\star}}:=\alpha_{i^{\star},j^{\star}}/\sigma$ and multiply (12) by $\gamma_{i^{\star},j^{\star}}(b)$ to obtain

[TABLE]

In this equation we replace

[TABLE]

by

[TABLE]

Since $\gamma_{i^{\star},j^{\star}}\,\beta_{i^{\star},j^{\star}}$ involves fewer non-admissible automorphisms for $\tilde{f}_{i}=0$ than $\alpha_{i^{\star},j^{\star}}$, iteration of this substitution process will rewrite equation (12) in such a way that every $\alpha_{i,j}(\operatorname{ld}(\tilde{f}_{i}))$ involving non-admissible automorphisms for $\tilde{f}_{i}=0$ will be less than $\alpha_{i^{\star},j^{\star}}(\operatorname{ld}(\tilde{f}_{i^{\star}}))$ with respect to $\succ$. A further iteration of this substitution process will therefore produce an equation as (12) with no $\alpha_{i,j}$ involving any non-admissible automorphisms for $\tilde{f}_{i}=0$.

This shows that for every $\tilde{f}\in(E:Q)\setminus\{0\}$ there exists a Janet divisor of $\operatorname{ld}(\tilde{f})$ in the passive set defined by $\tilde{f}_{1}=0$, …, $\tilde{f}_{s}=0$. ∎

Let $\Omega\subseteq\mathbb{R}^{n}$ be open and connected and fix $\mathbf{x}\in\Omega$. Denoting the grid in (2) by $\Gamma_{\mathbf{x},h}$, we define

[TABLE]

and for a system $\tilde{S}$ of partial difference equations and inequations as in (11) we define the solution set

[TABLE]

Definition 6.4.

Let $\tilde{S}$ be a finite difference system over $\tilde{\mathcal{R}}$ and $\succ$ a ranking on $\tilde{\mathcal{R}}$. A difference decomposition of $\tilde{S}$ is a finite collection of quasi-simple difference systems $\tilde{S}_{1}$, …, $\tilde{S}_{r}$ over $\tilde{\mathcal{R}}$ such that $\operatorname{Sol}_{\Omega,\mathbf{x},h}(\tilde{S})=\operatorname{Sol}_{\Omega,\mathbf{x},h}(\tilde{S}_{1})\uplus\ldots\uplus\operatorname{Sol}_{\Omega,\mathbf{x},h}(\tilde{S}_{r})$.

In the following algorithm, Decompose in step 11 refers to an algorithm which computes a smallest superset of $G=\{\tilde{f}_{1},\ldots,\tilde{f}_{s}\}$ in $\tilde{\mathcal{R}}$ that is Janet complete as defined on page 1 (see also Section 3).

Theorem 6.5.

Algorithm 3 terminates and is correct.

Proof.

Algorithm 3 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 3 is equivalent to the condition that $Q=\emptyset$ holds after finitely many steps.

Apart from step 1, new systems are inserted into $Q$ in steps 18 and 20. We consider the systems that are at some point an element of $Q$ as the vertices of a tree. The root of this tree is the input system $\tilde{S}$. The systems which are inserted into $Q$ in steps 18 and 20 are the vertices of the tree whose ancestor is the system $L$ that was extracted from $Q$ in step 3 which in the following steps produced these new systems. Since the for loop beginning in step 5 terminates, the degree of each vertex in the tree is finite. We claim that every 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 resulted from a non-trivial difference reduction in step 9. When this new system will be extracted from $Q$ in a later round, a decomposition into quasi-simple algebraic systems will be computed in step 4. This may produce new branches of the tree, but along any of these branches, after finitely many steps the condition $a=$ true in step 10 will hold, because the order of the shifts in leaders of the arising equations is bounded by the maximum order of shifts in leaders of the ancestor system $L$.

In case of step 18 we are going to show that after finitely many steps a difference equation is obtained whose leader has not shown up as a leader of an equation in any preceding system in the current branch of the tree. First of all, the passivity check (step 12) yielded an equation $\tilde{f}=0$, $\tilde{f}\in P\setminus\tilde{\mathcal{K}}$, which is Janet reduced modulo $J$. Hence, either $\operatorname{ld}(\tilde{f})$ is not contained in the multiple-closed set generated by $\operatorname{ld}(G)$, or there exists $(\tilde{f}^{\prime},\mu^{\prime})\in J$ such that $\operatorname{ld}(\tilde{f}^{\prime})$ is a Janet divisor of $\operatorname{ld}(\tilde{f})$, but the degree of $\tilde{f}$ in $\operatorname{ld}(\tilde{f})$ is smaller than the degree of $\tilde{f}^{\prime}$ in $\operatorname{ld}(\tilde{f}^{\prime})$. In the first case the above claim holds. The second case cannot repeat indefinitely: First of all, if $\operatorname{ld}(\tilde{f})=\operatorname{ld}(\tilde{f}^{\prime})$, then in a later round, either a pseudo-reduction of $\tilde{f}^{\prime}$ modulo $\tilde{f}$ will be performed if the initial of $\tilde{f}$ does not vanish, or $\operatorname{init}(\tilde{f})=0$ has been added as a new equation (with lower ranked leader). Since this leads to a sequence in $\operatorname{Mon}(\Sigma)$ which strictly decreases, infinite chains are excluded in this situation. If case $\operatorname{ld}(\tilde{f})\neq\operatorname{ld}(\tilde{f}^{\prime})$ occurs repeatedly, then a sequence $((\theta_{i}\,\tilde{u}^{(\alpha)})^{e_{i}})_{i=1,2,3,\ldots}$ of leaders of newly inserted equations arises, where $\theta_{i}\in\operatorname{Mon}(\Sigma)$, $\alpha\in\{1,\ldots,m\}$, $e_{i}\in\mathbb{Z}_{\geq 0}$, such that $e_{i+1}<e_{i}$ holds (and where also $\theta_{i}\mid\theta_{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$ if step 12 confirmed passivity. Such a system is quasi-simple as an algebraic system because (up to 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 6.2 is ensured by step 14. Hence, all difference systems in $T$ are quasi-simple. Splittings of systems only arise in step 4 by adding an equation $\operatorname{init}(\tilde{f})=0$ and the corresponding inequation $\operatorname{init}(\tilde{f})\neq 0$, respectively, to the two new systems replacing the given one. Since no solutions are lost or gained, this leads to a partition as required by Definition 6.4. ∎

7. s-consistency check

Recall that $\tilde{S}^{=}$ (resp. $\tilde{S}^{\neq}$) denotes the set of left hand sides of equations (resp. inequations) in a difference system $\tilde{S}$. We shall use the same notation for differential systems.

Clearly, if one approximates the partial derivatives occurring in a simple differential system $S$ by appropriate finite differences, then one obtains a w-consistent approximation $\tilde{S}$ to $S$ (cf. Sect. 2 and 9).

The following algorithm verifies s-consistency of such FDA.

Correctness of the algorithm follows from Definition 2.1 (extended to inequations) and from passivity of the output subsystems of Algorithm 3. Their solution spaces partition the solution space of the input FDA. Thereby, any subsystem $\tilde{L}_{i}$ in the output with $b_{i}=$ true is s-consistent with $L_{i}$, where $\tilde{L}_{i}\xrightarrow[h\rightarrow 0]{}{L}_{i}$ and w-consistent if $b_{i}=$ false. If $b_{i}=$ true for all $i$, then $\tilde{S}$ is s-consistent with $S$. Termination follows from that of the subalgorithms.

8. Illustrative Example

Example 8.1.

We consider the system of nonlinear PDEs

[TABLE]

which is a simple differential system, as it is easily checked that the cross-derivative $\partial_{y}(u_{x}-u^{2})-\partial_{x}(u_{y}+u^{2})$ reduces to zero modulo (13). We investigate the discretized system which is obtained by replacing $\partial_{x}$ and $\partial_{y}$ by the forward differences $D_{1}^{+}$, $D^{+}_{2}$, respectively:

[TABLE]

This system of nonlinear difference equations is simple as an algebraic system, but the passivity check reveals the consequence

[TABLE]

The continuous limit of $\tilde{u}_{i,j}^{4}$ for $h\to 0$ is the differential polynomial $u^{4}$, 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 $\partial_{x}$ and $\partial_{y}$ by $D_{1}^{+}$ as before and the backward difference $D_{2}^{-}$, respectively:

[TABLE]

In order to avoid negative shifts, we replace equation $(E)$ by $\sigma_{2}(E)$. Then this system of nonlinear difference equations is simple because it is algebraically simple and the passivity check yields

[TABLE]

We conclude that FDA (15) is s-consistent with system (13).

9. Navier-Stokes Equations

Example 9.1.

The Navier-Stokes equations for a three-dimensional incompressible viscous flow in vector notation are

[TABLE]

where $\mathbf{x}=(x_{1},x_{2},x_{3})$, $\mathbf{u}(\mathbf{x},t)$ is the velocity vector $\mathbf{u}=(u,v,w)$, $p(\mathbf{x},t)$ is the pressure and $\mathrm{Re}$ is the Reynolds number. For the ranking $\succ_{{\rm TOP},{\rm lex}}$ (Example 5.2) such that

[TABLE]

the (non-admissible) prolongation $\nabla\cdot\partial_{t}\mathbf{u}=0$ of the right (continuity) equation in (16) and its reduction modulo the left (momentum) equation yields the pressure Poisson equation

[TABLE]

which is the integrability condition (cf. (Seiler'10, ), p.50) to (16). Clearly, the differential system (16) and (18) satisfies the simplicity conditions (1)–(4) in Definition 4.1. Now we consider the following class of FDA to (16) defined on the four-dimensional grid (2)

[TABLE]

where $D_{t}$ approximates $\partial_{t}$, $\mathbf{D}=(D_{1},D_{2},D_{3})$ approximates $\nabla$ and ${\tilde{\Delta}}$ approximates $\Delta$. It is clear that system (19) is w-consistent with (16). If one considers the difference analogue of ranking (17) satisfying

[TABLE]

then completion of (19) to a passive form by Algorithm 3 is equivalent to enlargement of this system with the integrability condition

[TABLE]

Eq. (21) approximates Eq. (18) and can be obtained, in the full analogy with the differential case, by the prolongation $\mathbf{D}\cdot D_{t}\tilde{u}=0$ of the discrete continuity equation in system (19) and its reduction by the discrete momentum equation.

The left-hand sides of Eqs. (16) and (18) form a difference Gröbner basis of the ideal generated by Eqs. (19) in $\mathbb{Q}(\mathrm{Re,h})\{\mathbf{\tilde{u}},\tilde{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 the Navier-Stokes and pressure Poisson equations in the two-dimensional case as well. Examples of such FDA were studied in (ABGLS'13, ). One more s-consistent two-dimensional FDA was derived in (ABGLS'17, ). In its approximation of Eq. (18) the redundant to zero term $-\frac{\tilde{\Delta}\,(\nabla\cdot\mathbf{\tilde{u}})}{\mathrm{Re}}$ was included in the left-hand side of (21). This inclusion improves the numerical behavior of FDA (cf. (Rempfer'06, ), Sect.3.2).

Example 9.3.

For the two-dimensional system (16), (18) with grid velocities $(u,v)$ and pressure $p$ we consider the discretization

[TABLE]

where

[TABLE]

and $i\in\{1,2\}$. Then FDA (22) is w-consistent with (16) and (18). However, it is s-inconsistent since $\tilde{e}^{(4)}\not\in\llbracket\tilde{\mathcal{I}}\rrbracket$ where $\llbracket\tilde{\mathcal{I}}\rrbracket\subset\tilde{\mathcal{R}}$ is the perfect closure (see Lemma 2.6) of the ideal generated by $\{\tilde{e}^{(1)},\tilde{e}^{(2)},\tilde{e}^{(3)}\}$. It follows, as modulo $\llbracket\tilde{\mathcal{I}}\rrbracket$ the equality holds

[TABLE]

whereas the difference operator $\mathbf{D}\cdot\mathbf{D}$ in (21) is not equal to $\tilde{\Delta}$:

[TABLE]

10. Conclusions

In this paper, for the first time, we devised a universal algorithmic approach to check s(trong)-consistency of a system of finite difference equations that approximates a polynomially nonlinear PDE system on a Cartesian solution grid. In our earlier paper (GR'10, ) we studied this problem for linear PDE systems and showed how to check their s-consistency by using differential and difference Gröbner bases of ideals generated by the polynomials in PDE and FDA. As this takes place, all related computations can be done, for example, with the Maple packages LDA (GR'12, ) and Janet (Maple-Janet'03, ).

Extension of the Gröbner basis method to the nonlinear case is not algorithmic due to the non-Noetherity of differential and difference polynomial rings. On the other hand, the differential Thomas decomposition (Def. 5.6) and its difference analogue (Def. 6.4) are fully algorithmic (cf. (GLHR'18, ; BGLHR'12, ; Robertz6, ) and Alg. 3). These decompositions are essentials of the s-consistency check (Alg. 4). The differential Thomas decomposition is built into Maple 2018 and its implementation 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 system and for a (w-consistent) FDA Algorithm 4 returns false, as it takes place in Example 9.3, then we have to try another FDA and check the s-consistency again. In doing so, if we know a minimal generating set for the radical differential ideal generated by the input simple differential system, then its FDA should be tried as an input for Algorithm 3. Such is indeed the case with the Navier-Stokes equations (Ex. 9.1), for which Algorithm 3 returns s-consistent discretization (19), (21) if it is applied to Eqs. (16) and ranking (20).

However, the choice of FDA to the minimal generating set for the simple differential system as an input for Algorithm 3 not always yields s-consistent FDA, as demonstrated by Example 8.1. In addition, designing an algorithm for construction of a minimal generating set for an ideal is an open problem for commutative polynomial rings and is probably unsolvable in the differential case.

In applications of finite difference methods to PDE systems which have integrability conditions, it is important not only to preserve these conditions at the discrete level, but to ensure also that FDA is s-consistent with the PDE system. FDA (19), (21) to the Navier-Stokes equations (16) satisfies this requirement and for this reason it 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 supported by the Russian Foundation for Basic Research (grant No. 18-51-18005) and by the RUDN University Program (5-100).