Arithmetic of D-Algebraic Functions

Let us now assume the univariate setting. To define our target ODE system, we use the multivariate differential polynomial ring $S_{\mathbf{y},z}$, where $\mathbf{y}=(y_{1},\ldots,y_{n})$ for some $n\in\mathbb{N}$. Our main interest is in systems of the form

where $\mathbf{A}=\left(A_{1},\ldots,A_{n}\right)\in\mathbb{K}(x)(y_{1},\ldots,y_{n})^{n},$ $B\in\mathbb{K}(x)(y_{1},\ldots,y_{n})$. In the context of differential algebra, we consider the following $n+1$ polynomials

where $Q$ is the least common multiple of the denominators in $(\mathcal{M})$, and $A_{i}=a_{i}/Q,$ $i=1,\ldots,n,B=b/Q$. With an appropriate monomial ordering over the differential polynomial ring $S_{\mathbf{y},z}$, multiplying by $Q$ in (2) is equivalent to clearing denominators and considering $Q$ as the least common multiple of the initials of the resulting differential polynomials. We make the above choice for convenience as it allows us to defer the selection of a monomial ordering to a later stage, which also appears more suitable for the algorithmic part (see also [12]). Next, we consider the differential ideal

where $:Q^{\infty}$ denotes the saturation with $\{Q\}$. We define $I_{\mathcal{M}}^{\leq j}$ as the algebraic ideal containing the polynomials in (2) and their first derivatives of order at most $j$ saturated with $\{Q\}$. The system $(\mathcal{M})$ may be seen as a single-output state-space model without input. General setting and computations of such models are discussed in [17, 34]. There it is also explained how so-called input-output equations of $(\mathcal{M})$, i.e. ADEs that relate the output variable $z$ and the input variable (usually denoted $u$), are deduced from the elimination ideal

Hence in our setting, the input-output equations correspond to ADEs satisfied by the rational function $B(\mathbf{y})$.

A system of the type of $(\mathcal{M})$ is called rational dynamical system of dimension $n$ (the dimension of $\mathbf{y}$). In their work [33], Pavlov and Pogudin developed an algorithmic approach to derive such systems from first-order ADEs of two dependent variables, which represent the output and the input. The relationship between these systems and the closure properties of D-algebraic functions is discussed in [1]. Specifically, when the dynamical system $(\mathcal{M})$ models ADEs fulfilled by some D-algebraic functions $f_{1},f_{2},\ldots,f_{N}$, such that $B$ in $(\mathcal{M})$ represents a rational expression in the dependent and independent variables, the corresponding input-output equations are ADEs satisfied by $B(f_{1},\ldots,f_{N})$. In particular, thanks to this construction and the non-trivial nature of the elimination ideal that generates the input-output equations of $(\mathcal{M})$, we can reestablish that the set of D-algebraic functions is a field.

We will say that a computed algebraic (partial or ordinary) differential equation is of the least order if its order is the minimal order differential polynomial in the corresponding differential ideal. The arithmetic of two D-algebraic functions in the l.h.o. and non-l.h.o. cases were treated in [1]. In this paper, we improve the computational bound given in that paper for non-l.h.o. algebraic ODEs. We demonstrate that the bound established for l.h.o. ADEs also applies to non-l.h.o. ADEs. In essence, the sum of the orders of some given ADEs serves as an upper bound for the minimal order of an ADE satisfied by rational expressions formed from generic solutions of those ADEs.

We here recall the result from [1] for the arithmetic of D-algebraic functions that are zeros of l.h.o. differential polynomials. Since each $p_{i}$ is l.h.o, we can write its highest derivative of the dependent variable in terms of the lower ones from the associated ADE. We write

Let us now construct the ODE system (here a dynamical system) that models the operation $r(f_{1},\ldots,f_{N})=:f$.

Define the indeterminates

This implies that

Furthermore, let

Set $M\coloneqq\sum_{i=1}^{N}n_{i}$. The following $M$ dimensional dynamical system models $f$:

We denote the corresponding differential ideal by $I_{\mathcal{M}_{f}}\in S_{\mathbf{w},z}$. The above construction sustains the following proposition, which extends [1, Proposition 4.9.] for arbitrarily many differential polynomials.

Suppose there are $N_{1}\in\mathbb{N}\setminus\{0\}$ differential polynomials that are non-l.h.o. among the $p_{i}$’s, $N_{1}\leq N$. In [1] this situation is overcome by replacing those differential polynomials with their first derivatives. This brings the problem back to the l.h.o. case, but increase the bound for the order of the differential equation sought by $N_{1}$. The main issue is that without taking derivatives, we do not easily relate the arithmetic in this case to a rational dynamical system.

Our main observation is that, provided some algorithmic changes and a little theoretical adaptation, the arguments over which the computations from the dynamical systems rely are still valid. Indeed, suppose that $p_{j}$, $1\leq j\leq N$ is a non-l.h.o. differential polynomial, and let $m_{j}$ be the degree of $y_{j}^{(n_{j})}$ in $p_{j}$. We can rewrite $p_{j}$ as

where $c_{m_{j}}\in L\coloneqq\mathbb{K}[x,y_{j}^{(\leq n-1)}]$, $p_{j,1}\in L[y_{j}^{(n_{j})}]$ such that $\deg_{y_{j}^{(n_{j})}}(p_{j,1})<m_{j}$. The coefficient $c_{m_{j}}$ is often called the initial of $p_{j}$. We make the same change of differential indeterminates as in Section 2.1 and build an ODE system of the form

where $m_{j}=\deg_{y_{j}^{(n_{j})}}(p_{j}),j=1,\ldots,N$. The rational functions $A_{j},j=1,\ldots,N$ result from solving $p_{j}=0$ for ${y_{j}^{(n_{j})}}^{m_{j}}$, i.e., with the same variable names in the expressions of $A_{j}$ and $p_{j}$, we can write

Observe that the denominators are free of the indeterminate representing $y_{j}^{(n_{j})}$. This allows us to define $Q$ in $(\mathcal{M}_{f}*)$ as in $(\mathcal{M})$ by using all the denominators of the system. To ease the theoretical understanding, let us simplify the notation by rewriting $(\mathcal{M}_{f}^{*})$ in the following abstract form

where

• •

  $\mu_{i},i=1,\ldots,M$, is either $1$ or one of the $m_{j},j=1,\ldots,N$ in $(\mathcal{M}_{f}^{*})$;

• •

  $\mathcal{A}_{i}$ (with numerator $a_{i}$) is either $w_{i+1}$ or the part of $A_{j}$, that is free of $w_{j}^{\prime}$ in $(\mathcal{M}_{f}^{*})$, $j=1,\ldots,N$;

• •

  $\mathcal{E}_{\mathcal{A}_{i}}(w_{i}^{\prime})$ (with numerator $e_{a_{i}}$) is either $0$ or the part of $A_{j}$ that contains $w_{j}^{\prime}$ (in its numerator) in $(\mathcal{M}_{f}^{*})$, $j=1,\ldots,N$;

• •

  $B(w_{1},\ldots,w_{M})=r(w_{1},w_{n_{1}+1},\ldots,w_{\sum_{i=1}^{N-1}n_{i}+1})$ with numerator $b$.

We call a system of the form of (10) a radical-rational dynamical system. Note that the shape of $(\mathcal{M}_{f}^{*})$ imposes certain constraints on the radical-rational dynamical systems we use. In the context of differential algebra, we consider the differential ideal

where $H$ is the least common multiple of $Q$ and the separants of $Q{\mathbf{w}^{\prime}}^{\mu}-\mathbf{a}(\mathbf{w})-\mathbf{e}_{\mathbf{a}}(\mathbf{w}^{\prime})$. These separants are explicitely given by

where derivations are taken component wise.

Let us now present our algorithmic steps for computing ADEs satisfied by rational expressions of univariate D-algebraic functions.

To prove the correctness of Algorithm 1, we must show that the elimination ideal $I_{\mathcal{M}_{f}^{*}}^{\leq M}\cap\mathbb{K}[x][z^{(\infty)}]$ in item 8 is non-trivial. This fact is established by the following theorem.

From the proof of Theorem 2.1, one should observe that $M$ is the minimal integer for which the used arguments hold. This relates to the minimality of the order of the ADE obtained after the computations. Moreover, the l.h.o and the non-l.h.o cases are unified in (10). We make $\mathbf{\mathcal{E}}$ explicit to illustrate the difference with the construction of $(\mathcal{M}_{f})$ in Section 2.1.

The main point behind Algorithm 1 is the systematic construction of ADEs for rational expressions of D-algebraic functions using ODE systems (or triangular sets). Note that the saturation with respect to the least common multiple of the denominators and the separants implies that the D-algebraic functions we deal with are the generic solutions of ADEs.

We have implemented Algorithm 1 in the NLDE package (see [43, 41]). It complements the implementation for arithmetic operations from [1]. For full details on the syntax, we refer the reader to the documentation at MathRepo NLDE Doc (see also [41]). We here present examples to illustrate what can be seen as an improvement of the implementation in [1]. The commands for arithmetic operations are NLDE:-arithmeticDalg and NLDE:-unaryDalg. The former is for several input ADEs, and the latter is for a single input ADE. By default, NLDE:-arithmeticDalg follows the second method in [1], which implies that in the non-l.h.o. case, the order of the returned ADE may be higher than expected. When some of the input ADEs are not l.h.o., one specifies lho=false to use Algorithm 1 and obtain an ADE of the smallest order possible. Regarding Gröbner bases elimination, by default,

• •

  when lho=true, we use the pure lexicographic ordering plex from the Maple package Groebner; this is the same ordering used in Theorem 2.1, which guarantees the desired result.

• •

  When lho=false, we use the lexdeg elimination order of Bayer and Stillman [4] implemented by the EliminationIdeal command of the Maple package PolynomialIdeals. This is a block order where blocks are ordered with plex. In our case, there are two blocks: the variables of the elimination ideal constitute one block, and the remaining variables constitute the other one. Both blocks are internally ordered with a suitable order like degrevlex for instance. Although this ordering generally returns the desired result, it can also fail to find it sometimes. Nevertheless, lexdeg tends to provide a better efficiency. To use plex when lho=false, one further specifies lhoplex=true.

To use the lexdeg ordering when lho=true, one specifies ordering=lexdeg. For unary operations, NLDE:-unaryDalg was updated to only implement the method of this paper.

Note that either of the implementations of the Rosendfeld-Gröbner algorithm and the Thomas decomposition could not find (21) within a reasonable time. The former runs out of memory, and the latter keeps running even after an hour of computations. They tend to behave this way for ADEs of higher degrees. We also mention that even in the situation where one of these algorithms would return a correct result, the exponent in (22) could not be found since those algorithms output generators of a radical differential ideal, unlike our algorithm which computes in the corresponding differential ideal. It would be interesting to know what these exponents mean.

We are going to introduce a total ordering on the set of derivative operators with the aim of using them in increasing order up to a given limit. Algorithmically, this will entail to defining a map that sends every non-negative integer to a derivative operator. In [25] and many references in differential algebra (see also [15]), the notion of order for partial differential rings is different from the one we use. The prevalent definition does not easily adapt to our computational goal. Indeed, given the set of derivations $\Delta_{2}=\{\partial_{x_{1}},\partial_{x_{2}}\}$, one considers the free multiplicative semigroup $\Theta_{2}$ generated by the elements of $\Delta_{2}$. To any derivative operator

one defines the order $s=n_{1}+n_{2}$. The derivative $\theta u$ of $u\in R$, is called the derivative of $u$ of order $s$ (supposedly w.r.t. $\theta$). The order of a differential polynomial $p\in R[y^{(\infty,\infty)}]$ (also called $\Delta_{2}$-polynomial ring) is then the maximum $n_{i}+n_{j}$ such that $y^{(n_{i},n_{j})}$ occurs in $p$. What we wish to add in this definition is a total ordering (or ranking) of the derivative operators in $\Theta_{2}$. Observe that for $\Delta_{1}=\{\partial_{x}\}$, there is a natural ranking of the derivative operators:

allowing to uniquely associate a non-negative integer with the order of a particular derivative. With the above definition of order for bivariate partial differential polynomial rings, there are two possible choices for the derivative of $y$ of order $1$, namely, $y^{(0,1)}$ and $y^{(1,0)}$. We combine total order and classical orderings (lexicographic or reverse-lexicographic).

The total-order is used to ensure the possibility to move from $y^{(n_{1},0)}$ to $y^{(0,n_{2})}$ by some applications of the chosen derivation rule. In fact, the bijective map between $\mathbb{N}$ and $\mathbb{N}^{2}$ should be “obvious” by construction. Thus what we call order for $y^{(n_{1},n_{2})}$ is either the non-negative integer whose image (or preimage) is $(n_{1},n_{2})$ through that bijection, or simply $(n_{1},n_{2})$, assuming $\Theta_{2}$ is well-ordered (might be implicit for non-computational purposes). We will use the latter, and define a derivation rule based on a bijective map between $\mathbb{N}^{2}$ and $\mathbb{N}$.

Our idea is to use a “very” classical map in set theory. Consider the following start of a ranking for the application of derivative operators

This corresponds to the following ranking of derivative orders

One can associate each ordered pair of integers with a unique positive integer. Indeed, the underlying map is the so-called Cantor pairing function, which is a bijection from $\mathbb{N}^{2}$ to $\mathbb{N}$. The underlying ordering of the ordered pairs is the graded (or degree) lexicographic ordering; however, to avoid confusion with the lexicographic ordering that we will use for monomials, we use the terminology of Cantor pairing function or Cantor $l$-tuple functions for arbitrarily many independent variables. As proven by Rudolf Fueter and George Pólya [16], the Cantor pairing function is the only quadratic polynomial that maps $\mathbb{N}^{2}$ to $\mathbb{N}$. Therefore for two variables, it appears as the most “natural and efficient” choice for processing speed on a regular computer.

We see that such a ranking is inherent to the definition of a total order on $\Theta_{2}$ to rank the variable in $S_{y}$. Note that this ranking is in line with the terminology introduced by Ritt and should not be confused with monomial orderings. Let $\sigma_{2}$ be the Cantor pairing function, $\sigma_{2}^{-1}$ its functional inverse such that $\sigma_{2}^{-1}(k)=(\sigma_{2,1}^{-1}(k),\sigma_{2,2}^{-1}(k))$. For $k\in\mathbb{N}$, the $k$th derivation, denoted $\theta_{x_{1},x_{2}}^{k}$, is:

We view $R[y^{(\infty,\infty)}]$ as

This implicitly defines an operator which we denote by $\theta_{x_{1},x_{2}}$.