<figure><img src="image_1.png" /></figure>

HAL
open science

# A Mathematica Package for Certifying the Nonexistence of Darboux Polynomials

Maxime Bridoux, Khalil Ghorbal

► To cite this version:

Maxime Bridoux, Khalil Ghorbal. A Mathematica Package for Certifying the Nonexistence of Darboux Polynomials. 2024. hal-04818282

HAL Id: hal-04818282

https://inria.hal.science/hal-04818282v1

Preprint submitted on 4 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Mathematica Package for Certifying the Nonexistence of Darboux Polynomials

Maxime Bridoux Khalil Ghorbal

Inria, France

firstname.lastname@inria.fr

## Abstract

We present a Mathematica package to help investigate the existence of Darboux polynomials for (polynomial) ordinary differential equations (ODE). If we denote the derivation encoded by a given ODE by $D$, the package infers and propagates necessary conditions for a *generic* ansatz $p$ to be a Darboux polynomial for $D$. By generic, we mean that not only the coefficients of $p$ are unknown, but also the exponent $d$ of its leading monomial. We present some examples for which the package produces a formal certificate proving the nonexistence of nontrivial Darboux polynomials. The package complements well standard algorithms that search for Darboux polynomials up to a fixed bound on the total degree.

## 1 Introduction

In his seminal work, Gaston Darboux introduced algebraic *particular integrals*, known today as Darboux polynomials [Gor01], as a mean to construct *general integrals* for polynomial ODEs of the standard form $\dot{x}_i = f_i(x_1, \dots, x_n)$, $i = 1, \dots, n$, where $f_1, \dots, f_n$ are multivariate polynomials in $x_1, \dots, x_n$ over some field, $\dot{x}_i$ denotes the derivative of $x_i$ with respect to an independent variable $t$. From a differential algebraic perspective, the ODE defines a *derivation* $D = \sum_{i=1}^n f_i \partial_i$, acting on the ring of polynomials where $\partial_i$ denotes the partial derivative with respect to $x_i$. A polynomial $p$ is *Darboux for D*, or simply *Darboux* when $D$ is clear from the context, whenever $D(p) = qp$ for some polynomial $q$, called the *cofactor* of $p$. (Polynomials of total degree zero are trivially Darboux.)

Darboux generation algorithms (e.g. [Man93]) are semi-decision procedures enumerating all Darboux polynomials up to a certain fixed bound on the total degree. The bound is eventually incremented until finding a (not necessarily irreducible) Darboux polynomial or reaching memory and/or time limits. The Mathematica package[^1] we present is primarily designed to prove the nonexistence of nontrivial Darboux polynomials. The necessary conditions it infers complements however generation algorithms by providing useful constraints on both the degree and support of a generic ansatz $p$ to be Darboux. For the proofs and further details about the implemented algorithms, we refer the reader to [GB24].

<footnote>[^1]: Available in https://gitlab.com/maximebridoux/darbouxcertification.</footnote>
## 2 Case studies

Let $x = (x_1, \dots, x_n)$ denote a set of variables and $\alpha = (\alpha_1, \dots, \alpha_n)$ be a vector of natural numbers. We use the shorthand notation $x^\alpha$ to denote the multivariate monomial $\prod_{i=1}^n x_i^{\alpha_i}$. Given a monomial order, the *multidegree* of a multivariate polynomial $p$ is the exponent of its leading monomial, denoted $\text{LM}(p)$. An ansatz $p = \sum_\alpha a_\alpha x^\alpha$ is said to be *generic* when its multidegree $d$ is not fixed. We represent a generic polynomial implicitly by a partial function 'a' which associates to each exponent $\alpha$ the formal symbol $a_\alpha$. Known expressions $a(\alpha)$ of coefficients $a_\alpha$ are stored in a finite list of *substitution rules* of the form $a_\alpha \to a(\alpha)$. Without loss of generality, we suppose that the generic ansatz $p$ is monic, that is the coefficient of $\text{LM}(p)$ is set to 1. This list of substitution rules is thus initially set to $\{a_d \to 1\}$.

For each example of the following section, our package provides a dedicated Mathematica notebook (version 13.3). All computations (for these benchmarks) take few seconds to run on a standard laptop running Ubuntu 22.04.01 and equipped with an i7-12700H CPU together with 32GB of memory. We refer to [GB24] for further details about the implemented algorithms.

### 2.1 Van der Pol Oscillator

The dynamic of the Van der Pol oscillator is defined by

$$ 
\begin{aligned}
\dot{x}_1 &= x_2 \\
\dot{x}_2 &= \mu(1 - x_1^2)x_2 - x_1
\end{aligned}
$$

which we implement as a key-value association. Monomial orders are encoded via their associated weight matrix [Rob85].

```
In[1]:= vdp = <| x1 <- x2, x2 <- µ (1 - x1^2) x2 - x1 |>;
lex = {{1,0}, {0,1}}; dlex21 = {{1,1}, {0,1}};
```

The routine `ComputeQuotient(ode, d, a, morder)` returns the simplified quotient $q$ of the division of $D(p)$ by $p$ with respect to the monomial order `morder`. The symbol `a`, together with the multidegree `d`, provide an implicit encoding of the generic ansatz $p$. Hence $a[\{1 + d[1], -1 + d[2]\}]$ represents the coefficient $a_{d_1+1,d_2-1}$ of $p$. The polynomial $D(p)$ is computed only partially on the fly: only monomials of $D(p)$ that are divisible by $\text{LM}(p)$ are relevant. We can effectively retrieve the result of [GB24, Proposition 5.6] (stated for $\mu = 1$) using the following command:

```
In[2]:= q = ComputeQuotient[vdp, {d[1], d[2]}, a, dlex21]
Out[2]= a[\{1 + d[1], -1 + d[2]\}] (1 + d[1]) + µ d[2] - x1^2 µ d[2]
```

Once a quotient $q$ is computed, the next step is to compute the coefficients of the remainder $r = D(p) - qp$. This is performed by the routine `Coeff(ode, d, a, q, alpha)`, which outputs the general expression of the coefficient of $x^\alpha$ in $r$. For $p$ to be Darboux, all the coefficients of $r$ have to vanish. By scanning the coefficients of $r$, one infers necessary conditions on $a_\alpha$ and $d$. The routine `LinearPropagation(ode, d, a, q, morder, iterationLimit, Substd*, Substp*, Conditions*)` implements the constants' propagation described in [GB24, Section 6]. It outputs a tuple of three lists of substitution rules for $p$ to be Darboux: (1) a list of necessary conditions on the multidegree $d$, (2) a list of necessary conditions on the coefficients $a_\alpha$, and (3) a formal *certificate* briefly described in the sequel. We run below the procedure for the Van der Pol dynamic with 6 as an iteration limit.
$$
\begin{align*}
\text{In[3]} &:= \{\text{substd, substp, certif}\} = \text{LinearPropagation}[\text{vdp}, \{\text{d}[1], \text{d}[2]\}, \text{a}, \text{q}, \text{dlex21}, 6] \\
\text{Out[3]} &:= \{\}, \{\text{a}[\{\text{d}[1], \text{d}[2]\}] \rightarrow 1\}, \{\}
\end{align*}
$$

This first run didn't infer any extra condition and simply returned the default values. In such
cases, one can either increase the iteration limit, or add extra initial assumptions on *d* or *a*<sub>α</sub> by
supplying a list of substitution rules via the optional arguments "Substd" and "Substp" respectively.
A closer look at the Van der Pol example reveals that the propagation requires some extra knowledge
on the parameter μ of the dynamic. Indeed, depending on μ, *D* might have a Darboux polynomial.
For instance, when μ = 0, the oscillator reduces to the standard pendulum for which the polynomial
*x*<sub>1</sub><sup>2</sup> + *x*<sub>2</sub><sup>2</sup> + *c* is Darboux for any constant *c*. The optional argument "Conditions" serves precisely this
purpose and allows the user to specify extra conditions on the parameters of the ODE. A second
run assuming μ ≠ 0 is shown below.

$$
\begin{align*}
\text{In[4]} &:= \{\text{substd, substp, certif}\} = \text{LinearPropagation}[\text{vdp}, \{\text{d}[1], \text{d}[2]\}, \text{a}, \text{q}, \text{dlex21}, 6, \\
& \quad \text{"Conditions"} \rightarrow \{\mu \neq 0\}] \\
\text{Out[4]} &:= \{\{\text{d}[1] \rightarrow 0, \text{d}[2] \rightarrow 0\}, \{\text{a}[\{0, 0\}] \rightarrow 1\}, \dots \}
\end{align*}
$$

The procedure outputs that for *p* to be Darboux, the multidegree *d* has to be (0, 0), reducing *p* to a constant. This concludes the proof of nonexistence of nontrivial Darboux polynomials. Note that 6 is the minimum iteration limit for which a complete proof is found. The corresponding certificate (which is truncated with dots in the previous output cell) is given as an ordered list of pairs (*d* + *η*, substitution rule), where *η* is a vector of integers. The pair means that the substitution rule was deduced from computing the coefficient *c*<sub>*d*+*η*</sub> of *x*<sup>*d*+*η*</sup> in *r* and solving the equation *c*<sub>*d*+*η*</sub> = 0. For instance, the first element of the certificate is:

$$
\begin{align*}
\text{In[5]} &:= \text{certif}[1] \\
\text{Out[5]} &= \{\{-1 + d[1], 1 + d[2]\}, a\{-3 + d[1], 1 + d[2]\} \rightarrow \frac{d[1]}{\mu}\}
\end{align*}
$$

To formally check the certificate a posteriori, one can compute the coefficient $c_{d+\eta}$ of $r$ and verify
that the associated substitution rule can indeed be deduced. For instance, for $c_{-1+d_1,1+d_2}$ one gets
the (linear) constraint $a_{-3+d_1,1+d_2} = \frac{d_1}{\mu}$.

$$
\begin{align*}
\text{In[6]} &:= \text{Expand}[\text{Coeff}[\text{vdp}, \{\text{d}[1], \text{d}[2]\}, \text{a}, \text{q}, \{-1 + \text{d}[1], 1 + \text{d}[2]\}, \text{dlex21}]] \\
\text{Out[6]} &= -\mu \, a[\{-3 + \text{d}[1], 1 + \text{d}[2]\} + \text{d}[1]
\end{align*}
$$

2.2 Liénard systems

Our implementation can also infer necessary conditions on the parameters of a given ODE for a
Darboux polynomial to exist. Given *f*, *g* univariate polynomials, the class of Liénard systems is
given by the following ODE (the Van der Pol oscillator is in particular a typical Liénard system):

$$
\begin{align*}
\dot{x}_1 &= x_2 \\
\dot{x}_2 &= -f(x_1)x_2 - g(x_1) \, .
\end{align*}
$$
We provide a notebook example for the case where $\deg f = \deg g = 3$. Computing the quotient using the monomial order $\mathbf{DLex}_{21}$ and performing the constraints' propagation infer the substitution $\{ \mathbf{d}[1] \to \mathbf{0} \}$. When executed with 7 as the iteration limit, the routine prompts a dialog box suggesting that an equation with several possible solutions was found, and asks the user to select one solution to explore. In this example, one can either choose the first option $\{\mathbf{d}[2] \to \mathbf{0}\}$, which completes the ongoing proof (returning $d_1 = d_2 = 0$), or choose to enforce extra conditions on $f$ and $g$ (as we did for $\mu$ in the Van der Pol system). Assuming that $d_2 > 0$ and both $f_0, g_0$ (the constant coefficients of $f$ and $g$) are distinct from 0, the procedure is unable to find a bound on $d_2$ but generates the substitutions shown below.

$$ \ln[7] := \{\text{substd, substp, certif} \} = \text{LinearPropagation}[\text{lienard}, \{\mathbf{d}[1], \mathbf{d}[2]\}, \text{a, q, dlex21, 7, "Conditions"} \to \{\mathbf{f}[3] \neq 0, \mathbf{g}[3] \neq 0\}]; \text{substd} $$

$$ \text{Out}[7] = \{\mathbf{d}[1] \to \mathbf{0}, \mathbf{g}[3] \to \frac{\mathbf{f}[3] \ \mathbf{g}[0]}{\mathbf{f}[0]}, \mathbf{g}[1] \to \frac{\mathbf{f}[1] \ \mathbf{g}[0]}{\mathbf{f}[0]}, \mathbf{g}[2] \to \frac{\mathbf{f}[2] \ \mathbf{g}[0]}{\mathbf{f}[0]}\} $$

Instantiating $f$ and $g$ as suggested, we immediately get $p = x_2 - \frac{g_0}{f_0}$ as a Darboux polynomial for $d_2 = 1$.

## 2.3 Rössler system

The Rössler system is given by the following ODE featuring three real parameters $\mu_1, \mu_2$, and $\mu_3$:

$$ \begin{aligned} \dot{x} &= -y - z \\ \dot{y} &= x + \mu_1 y \\ \dot{z} &= xz + \mu_2 - \mu_3 z \end{aligned} $$

In [LZ02, Theorem 1], Llibre and Zhang characterize the conditions on $\mu_1, \mu_2$ and $\mu_3$ under which the system has Darboux polynomials. For instance, when $\mu_2 = 0$, $z$ is a Darboux polynomial with cofactor $x - \mu_3$. An important, and particularly tedious, part of their characterization was to prove that, when $\mu_2 \neq 0$, the system has no Darboux polynomial. Our package proves automatically the same result and provides a formal certificate as a chain of 134 substitution rules for the longest proof (cf. the associated notebook).

## References

[GB24] Khalil Ghorbal and Maxime Bridoux. On bounding the degree of irreducible darboux polynomials. Submitted to ISSAC 2024., 2024.

[Gor01] Alain Goriely. *Integrability and Nonintegrability of Dynamical Systems*. Advanced series in nonlinear dynamics. World Scientific, 2001.

[LZ02] Jaume Llibre and Xiang Zhang. Darboux integrability for the rössler system. *International Journal of Bifurcation and Chaos*, 12(02):421–428, 2002.

[Man93] Yiu-Kwong Man. Computing closed form solutions of first order odes using the prelle-singer procedure. *J. Symb. Comput.*, 16(5):423–443, 1993.

[Rob85] Lorenzo Robbiano. Term orderings on the polynomial ring. In *EUROCAL*, volume 204, pages 513–517, 04 1985.