arXiv:1801.09942v1 [physics.comp-ph] 30 Jan 2018

The MAPLE package TDDS for computing
Thomas decompositions of systems of nonlinear
PDEs
Vladimir P. Gerdt
gerdt@jinr.ru

Markus Lange-Hegermann
markus.lange.hegermann@rwth-aachen.de

Daniel Robertz
daniel.robertz@plymouth.ac.uk
Abstract
We present the Maple package TDDS (Thomas Decomposition of Differential Systems) for decomposition of polynomially nonlinear differential
systems, which in addition to equations may contain inequations, into a
finite set of differentially triangular and algebraically simple subsystems
whose subsets of equations are involutive. Usually the decomposed system
is substantially easier to investigate and solve both analytically and numerically. The distinctive property of a Thomas decomposition is disjointness of the solution sets of the output subsystems. Thereby, a solution of a
well-posed initial problem belongs to one and only one output subsystem.
The Thomas decomposition is fully algorithmic. It allows to perform important elements of algebraic analysis of an input differential system such
as: verifying consistency, i.e., the existence of solutions; detecting the arbitrariness in the general analytic solution; given an additional equation,
checking whether this equation is satisfied by all common solutions of the
input system; eliminating a part of dependent variables from the system
if such elimination is possible; revealing hidden constraints on dependent
variables, etc. Examples illustrating the use of the package are given.

1

Introduction

The Maple package TDDS (Thomas Decomposition of Differential Systems) is
applicable to a set of finite-order partial differential equations (PDEs) of the
form (cf. Section 3)
!
∂ j1 +···+jn uk
pi x1 , . . . , xn ; u1 , . . . , um , . . . , j1
, . . . = 0 , i = 1, . . . , s , (1)
∂x1 · · · ∂xjnn
where k = 1, . . . , m, uk = uk (x1 , . . . , xn ). It is assumed that the left hand
sides pi in (1) are polynomials in their arguments. The package also allows
1

enlargement of (1) with a set of inequations {q1 6= 0, . . . , qt 6= 0} where qk
(k = 1, . . . , t) are also polynomials in the independent variables x1 , . . . , xn ,
dependent variables u1 , . . . , um and their partial derivatives.
A constructive algebraic approach to study systems of the form (1) goes back
to the following classical theorem proved by Kovalevskaya [1] (cf. [2]).
Theorem 1 (Cauchy-Kovalevskaya theorem). Let the left hand sides in system (1) read
!
∂ j1 +···+jn uk
∂ ni ui
pi =
− Fi x1 , . . . , xn ; u1 , . . . , um , . . . , j1
,...
(2)
∂xn1 i
∂x1 · · · ∂xjnn
where j1 + · · · + jn = j ≤ ni , j1 < ni and all the functions Fi (not necessarily
polynomial) are analytic in a neighborhood of the point
xi = x0i ,

uk = u0k ,

0
rk;
j1 ,...,jn :=

∂ j1 +···+jn uk
∂xj11 · · · ∂xjnn

(i, k = 1, . . . , m = s) .
x1 =x01 ,...,xn =x0n

(3)
Then in some neighborhood of the point (x01 , . . . , x0n ) the PDE system pi = 0
(i = 1, . . . , s) has a unique analytic solution satisfying the initial conditions

uk = φk (x2 , x3 , . . . , xn ) ,




 ∂uk = φk; 1 (x2 , x3 , . . . , xn ) ,
∂x1
for x1 = x01 , k = 1, . . . , m ,
(4)
..........................................



n
−1
k

uk
 ∂ n −1
= φk; nk −1 (x2 , x3 , . . . , xn ) ,
∂x1 k

where all φ are arbitrary analytic functions of their arguments in a neighborhood of the point (x02 , . . . , x0n ) such that they take at this point the initial values
satisfying (3) and (4).

Riquier [3], Janet [4] and Thomas [5, 6] developed a framework for generalization of Theorem 1. Riquier introduced a ranking ≻ on partial derivatives (see
Section 3, Definition 4), called Riquier ranking (cf. Definition 5 or [7]). Given a
Riquier ranking ≻, a PDE system (1) is orthonomic [3] if each of its equations
is solved with respect to the highest ranked partial derivative occurring in the
equation, and hence has the form
pi = 0 ,

pi := δi uki − Fi ,

i ∈ {1, . . . , s} ,

ki ∈ {1, . . . , m} ,

(5)

where the highest ranked derivative in pi , i.e. its leader (see Section 3), is written
as δi uki with a differential operator δi . This derivative is called principal and the
derivatives occurring in Fi are called parametric. The orthonomic system (5)
is called passive if all its differential and algebraic consequences do not lead
to constraints on the parametric derivatives. It should be noted that a PDE
system in the Kovalevskaya form (2) is orthonomic and passive for a certain
Riquier ranking [5].
2

Riquier [3] proved the existence of analytic solutions for orthonomic and
passive systems of PDEs. Janet [4] designed algebraic criteria of passivity for
orthonomic systems in terms of monomials associated with the principal derivatives in accordance to the mapping
∂ j1 +···+jn
∂xj11 · · · ∂xjnn

7−→ xj11 · · · xjnn .

(6)

For k ∈ {1, 2, . . . , m} let Mk be the leading monomial set, the set of monomials
associated by (6) with the principal partial derivatives of uk . The Janet criteria
of passivity are based on a certain partition of variables for each w ∈ Mk :
{x1 , . . . , xn } = M(w, Mk ) ∪ N M(w, Mk ) .
Moreover, for a linear PDE system the algebraic criteria allow to transform it
algorithmically into a passive form.
By using the results of Riquier and Janet, Thomas [5] formulated a Cauchy
problem providing the uniqueness and existence of an analytic solution in terms
of the monomials associated with the parametric derivatives (complementary
monomials). Furthermore, Thomas [6, 8] generalized the ideas and methods of
Riquier-Janet theory [3, 4] to PDE systems of form (1). He showed that such
differential systems can be decomposed into finitely many passive subsystems.
In so doing, each of these subsystems which we call simple differential system1
has certain triangular structure and can be solved with respect to its leader such
that the solved system is passive and orthonomic. It admits [6, 8] posing of a
Cauchy problem with initial data generalizing those in (3), (4) and providing
the uniqueness and (for a Riquier ranking) existence of an analytic solution.
Based on the ideas of Janet and Thomas, the foundations of differential algebra were developed by Ritt [9]. Then Wu [10] (cf. also [11]) further developed
the characteristic set method introduced by Ritt. The first implementation of
Thomas decomposition for systems of algebraic equations and ordinary differential systems was developed by Wang [12, 13].
In our papers [14, 15] (see also the book [16]) the Thomas approach was
algorithmized and implemented in full generality in Maple. It should be emphasized that the Thomas decomposition is different from two other triangular
decompositions into triangular and passive differential subsystems based on Rosenfeld-Gröbner [17] and rif [18] algorithms, which are implemented in packages
provided in Maple. The principal distinction of the Thomas decomposition is
disjointness of the solution sets of the output subsystems. Moreover, the strategy for completing differential systems to passive ones is based on Janet’s criteria
mentioned above.
The present paper is organized as follows. In Section 2 we illustrate by an
explicit example of nonlinear PDE systems some features of the differential
Thomas decomposition. Basic theoretical aspects of simple differential systems
1 Thomas

in [6] called them passive standard systems.

3

are described in Section 3 including underlying definitions and statements. Section 4 presents a list of commands of the package TDDS, and some of them are
illustrated by examples in Section 5.

2

Thomas decomposition of nonlinear PDE systems

Among the features of a Thomas decomposition of a system of nonlinear PDEs
is the possibility to determine power series solutions of the system around a sufficiently generic point in a straightforward way, to decide whether another PDE
is a consequence of the system, and to solve differential elimination problems.
We illustrate these features with explicit examples.
As a first example we consider the following system of nonlinear PDEs for
one unknown function u(t, x)

ut − 6 u ux + ux,x,x = 0 ,
(7)
u ut,x − ut ux = 0 ,
which is a combination of the Korteweg-de Vries equation and a Wronskian
determinant expressing that u(t, x) is a product of a function of t and a function
of x (cf. also [16, Ex. 2.2.61]). We would like to determine all power series
solutions (around (0, 0))
u(t, x) =

∞ X
∞
X

aj,k

j=0 k=0

tj xk
,
j! k!

aj,k ∈ C ,

(8)

of system (7). The idea is to partition the set of unknown coefficients aj,k into
two subsets: Taylor coefficients whose values can be chosen arbitrarily (up to
certain genericity assumptions, cf. below), and Taylor coefficients whose values
are determined from the chosen values (as roots of certain univariate polynomials in general) by substituting (8) into the system of PDEs and comparing
coefficients. However, this procedure requires preparatory symbolic manipulations of (7) in order to ensure passivity (or formal integrability), i.e. to ensure
that comparing coefficients results in a complete set of conditions for aj,k .
The Thomas decomposition method transforms any given system of polynomially nonlinear partial differential equations (and/or inequations) into an
equivalent finite collection of so-called simple differential systems. The solution
sets of the resulting simple differential systems form a partition of the solution
set of the original PDE system. Moreover, each simple differential system is formally integrable and allows to solve for the coefficients aj,k in (8) successively.
Different strategies, e.g., concerning the order in which partial derivatives of
u(t, x) and the corresponding Taylor coefficients aj,k are dealt with are possible. More precisely, the left hand sides of equations (and inequations) can be
expressed as univariate polynomials in their highest ranked partial derivatives
with coefficients that are polynomials in lower ranked partial derivatives, which
4

can be recursively represented in the same way. The Thomas decomposition
method depends on a given ranking, i.e., a total order ≻ on the set of all partial
derivatives

 j+k
∂
u
j,
k
∈
Z
≥0
∂tj ∂xk
which respects differentiation, i.e., for all j1 , k1 , j2 , k2 , m, n ∈ Z≥0 we have
∂ j1 +k1 u
∂ j2 +k2 u
≻
k
j
∂t 1 ∂x 1
∂tj2 ∂xk2

∂ j2 +m+k2 +n u
∂ j1 +m+k1 +n u
≻
,
j
+m
k
+n
∂t 1
∂x 1
∂tj2 +m ∂xk2 +n

=⇒

and which does not admit infinitely descending chains, i.e.,
∂ j+k u
≻u
∂tj ∂xk

for all j, k ∈ Z≥0 ,

j + k > 0.

For the given example we choose the ranking ≻ which is defined by


or
 j1 + k1 > j2 + k2
j2 +k2
j1 +k1
u
∂
u
∂
( j1 + k1 = j2 + k2 and ( k1 < k2
≻ j2 k2
⇐⇒

∂tj1 ∂xk1
∂t ∂x
 or k = k and j < j )),
1
2
1
2

(9)

which is analogous to the degree-reverse lexicographical term ordering used for
Gröbner basis computations [19]. More generally, when dealing with n independent variables, the degree-reverse lexicographical ranking is defined by comparing the total differentiation order and in case of equality defining the greater
derivative to be the one involving least differentiation when comparing lexicographically from ∂n backwards to ∂1 .
Given a system of polynomially nonlinear partial differential equations (and/or
inequations) and a ranking ≻, the Thomas decomposition method computes
linear combinations of partial derivatives of equations so as to eliminate highly
ranked partial derivatives as much as possible. The polynomial reductions which
are applied for this purpose assume that the coefficient of the highest power of
the highest ranked derivative in the divisor does not vanish on the solution set
of the system. In general this requires case distinctions. The Thomas decomposition method splits systems by introducing new equations and the complementary inequations as necessary for performing reductions and ensuring the
equivalence of the resulting decomposition to the original system. Hence, one
may think of the construction of a Thomas decomposition as a binary tree whose
leaves correspond to the resulting simple differential systems.
For system (7) the resulting Thomas decomposition (with respect to ≻)
consists of three simple differential systems:

ut = 0 ,




ux,x,x − 6 u ux = 0 ,
(10)

u 6= 0 ,



ux,x 6= 0 ,
5



 ut − 6 u ux
ux,x


u

u =

=
=
6
=

0,
0,
0,

(11)

0.

(12)

In the above systems the highest ranked partial derivatives are underlined
(where not obvious). We conclude that the set of power series solutions (8)
of (7) is the disjoint union of the sets of power series solutions of the simple systems (10), (11), (12). For (10) we can choose the values of a0,0 , a0,1 , a0,2 ∈ C
subject to the condition that then u(t, x) and its second derivative with respect
to x are not identically zero. Then all other Taylor coefficients are uniquely
determined through the equations in (10):
for all j ≥ 1 , k ≥ 0 ,

aj,k = 0
a0,3+n = 6

n 
X
k=0


n
a0,k a0,n−k+1
k

for all n ≥ 0 .

For (11) we can choose the values of a0,0 , a0,1 ∈ C subject to the condition
that then u(t, x) is not identically zero. Then all other Taylor coefficients are
uniquely determined by
for all j ≥ 0 , k ≥ 2 ,

aj,k = 0
a1+m,n

m X
n   
X
m n
= 6
aj,k am−j,n−k+1
j
k
j=0
k=0

= 6

3

1+m

1+m+n
(1 + m)! a1−n
0,0 a0,1

for all m ≥ 0 , n ∈ {0, 1} .

The theory of simple differential systems

We recall the basic principles of simple differential systems and Thomas decomposition. For more details, we refer to [16].
We consider partial differential equations of the form p = 0 as well as inequations of the form q 6= 0, where p and q are polynomials in unknown functions u1 ,
u2 , . . . , um of independent variables x1 , x2 , . . . , xn and their partial derivatives.
Since the differentiation order of the equations to be dealt with by the Thomas
decomposition method is not precisely known before applying the method to
the given PDE system, the set of all expressions potentially occurring as left
hand sides is the smallest polynomial ring containing all partial derivatives of
u1 , u2 , . . . , um , viz., the differential polynomial ring
R := Q{u1 , . . . , um } := Q[ (uk )J | k ∈ {1, . . . , m}, J ∈ (Z≥0 )n ] ,
where the differential indeterminate uk = (uk )(0,...,0) represents the unknown
function uk (x1 , . . . , xn ) with the same name and, more generally, (uk )J for the
6

multi-index J = (j1 , j2 , . . . , jn ) ∈ (Z≥0 )n represents the partial derivative
∂ j1 +j2 +...+jn uk
∂xj11 ∂xj22 . . . ∂xjnn

.

The polynomial ring R is closed under the derivations ∂1 , . . . , ∂n acting as
∂i (uk )J := (uk )J+1i ,

J + 1i := (j1 , . . . , ji−1 , ji + 1, ji+1 , . . . , jn ) ,

with additive extension to R respecting the product rule of differentiation.
Hence, ∂i acts as the partial differential operator ∂/∂xi . The coefficient field
Q of rational numbers can also be replaced with a larger field containing Q
admitting n derivations of which then the derivations ∂i are extensions to R.
For example, the coefficient field can be chosen to be the field Q(x1 , . . . , xn ) of
rational functions in x1 , . . . , xn with the usual derivations, if the system to be
dealt with consists of PDEs with rational function coefficients.
Given a system of partial differential equations
p1 = 0 ,

p2 = 0 ,

...,

ps = 0 ,

(13)

where pi ∈ R, the consequences of (13) obtained by taking linear combinations
and derivatives form a differential ideal of R, viz., a (non-empty) subset of R
which is closed under taking linear combinations of its elements with coefficients
in R and under differentiation (i.e., application of ∂1 , ∂2 , . . . , ∂n ). The differential ideal in question is said to be generated by the left hand sides p1 , p2 , . . . ,
ps .
Not every differential ideal of R admits a finite set of generators. However,
the differential ideals corresponding to differential systems in a Thomas decomposition are finitely generated (cf. Theorems 3 and 8 below) which provides
termination of the decomposition algorithm (cf. [16]). Differential ideals of this
kind arise as follows.
Let p be a differential polynomial in the differential indeterminate u and f
be a complex analytic function of x1 , . . . , xn . We write p(f ) for the result of
substituting f and its partial derivatives ∂ J f for u and uJ in p, respectively,
J ∈ (Z≥0 )n . Similar notation can be introduced when substituting m analytic
functions f1 , . . . , fm for m differential indeterminates u1 , . . . , um .
Now, (pr )(f ) = (p(f ))r = 0 for some positive integer r implies p(f ) = 0. A
differential ideal of R which contains for each of its elements p also all differential
polynomials in R of which a power is equal to p is said to be radical. The differential ideal of all p ∈ R which vanish under substitution of any analytic solution
of a PDE system is a radical differential ideal. The following important theorem establishes a one-to-one correspondence between radical differential ideals
of R and solutions sets (with complex analytic functions on suitable domains)
of PDE systems which are defined over R.
Theorem 2 (Nullstellensatz of Ritt-Raudenbush, [9], Sects. II.7–11, IX.27). Let
I be the differential ideal of R generated by the left hand sides p1 , p2 , . . . , ps of
7

a PDE system. If a differential polynomial p ∈ R vanishes under substitution of
any analytic solution of p1 = 0, . . . , ps = 0, then some power of p is an element
of I.
Radical differential ideals are finitely generated in the following sense.
Theorem 3 (Basis Theorem of Ritt-Raudenbush, [9], Sects. I.12–14). For every
radical differential ideal I of R there exists a finite subset B of I such that I is
the smallest radical differential ideal of R which contains B.
Each differential system in a Thomas decomposition corresponds to a radical
differential ideal (cf. Theorem 8 below). The notions of simple differential system
and Thomas decomposition require a choice of a total ordering on the symbols
representing the unknown functions and their partial derivatives.
Definition 4. A ranking ≻ on R is a total ordering on
Θu := { (uk )J | k ∈ {1, . . . , m}, J ∈ (Z≥0 )n }
such that the following two conditions are satisfied, which express that derivatives of an unknown function are ranked higher than the unknown function itself
and that ≻ is compatible with applying derivations.
1. For all i ∈ {1, . . . , n} and k ∈ {1, . . . , m} we have ∂i uk ≻ uk .
2. For all i ∈ {1, . . . , n}, k, l ∈ {1, . . . , m}, K, L ∈ (Z≥0 )n , the implication
(uk )K ≻ (ul )L ⇒ ∂i (uk )K ≻ ∂i (ul )L holds.
Examples of rankings were given in Section 2. Of special interest to us are
Riquier rankings.
Definition 5. We define a Riquier ranking to be a ranking ≻ such that
(
)
∂ j1 +···+jn
∀ δ1 , δ2 ∈ Θ :=
j1 , . . . , jn ∈ Z≥0
∂xj11 · · · ∂xjnn
we have δ1 u ≻ δ2 u if the total differentiation order in δ1 u is greater than the
one in δ2 u for any dependent variable u, and the following condition holds for
all dependent variables v, w:
δ1 v ≻ δ2 v

=⇒

δ1 w ≻ δ2 w .

Assume that a ranking ≻ on R has been chosen. Then every non-constant
polynomial p ∈ R involves a symbol (uk )J which is maximal with respect to
≻. It is referred to as the leader ld(p) of p. The differential polynomial p can
be represented recursively as a polynomial in ld(p) with coefficients which are
polynomials in their leaders, etc. The coefficient of the highest power of ld(p)
in p is called initial of p and denoted by init(p). The formal derivative of p with
respect to ld(p) is called the separant of p and denoted by sep(p).
8

A reduction process for differential polynomials can now be introduced as
follows. Let p1 and p2 be non-constant differential polynomials in R. In case
ld(p1 ) = ld(p2 ) we consider the degrees d1 and d2 of p1 and p2 in v := ld(p1 ),
respectively. If d1 ≥ d2 , then
init(p2 ) p1 − init(p1 ) v d1 −d2 p2
is either constant, or has a leader which is ranked lower than v, or has the same
leader, but has smaller degree in ld(p1 ) than p1 . If d1 < d2 , no reduction of p1
modulo p2 is possible. If ld(p1 ) 6= ld(p2 ), but there exists J ∈ (Z≥0 )n such that
v := ld(p1 ) = ∂ J ld(p2 ), then
sep(p2 ) p1 − init(p1 ) v d1 −1 ∂ J p2
is either constant or has a leader which is ranked lower than v (because ∂ J p2
has degree one in v.) If no J ∈ (Z≥0 )n exists such that ld(p1 ) = ∂ J ld(p2 ), then
no reduction of p1 modulo p2 is possible.
This reduction process can be adapted so as to eliminate any occurrence (in
sufficiently high degree) of symbols (uk )J which are leaders or derivatives of
leaders of p1 , . . . , ps ∈ R in a given differential polynomial p ∈ R. We say that
p ∈ R reduces to zero modulo p1 , . . . , ps ∈ R and their derivatives if p can be
reduced to the zero polynomial in this way.
Definition 6. A system of partial differential equations and inequations
p1 = 0 ,

... ,

ps = 0 ,

q1 6= 0 ,

... ,

qt 6= 0 ,

(s, t ∈ Z≥0 )

(14)

where p1 , . . . , ps , q1 , . . . , qt are non-constant differential polynomials in R, is
said to be simple (with respect to ≻) if the following conditions are satisfied.
1. The leaders of p1 , p2 , . . . , ps , q1 , q2 , . . . , qt are pairwise different.
2. Let v1 ≻ v2 ≻ . . . ≻ vk be the elements of Θu which effectively occur in the differential polynomials p1 , . . . , ps , q1 , . . . , qt . We consider
(14) as a system of polynomial equations and inequations in v1 , v2 , . . . ,
vk . If r ∈ {p1 , . . . , ps , q1 , . . . , qt } has leader vℓ , thus is a polynomial
r(vℓ , vℓ+1 , . . . , vk ), then we require that for every solution (a1 , a2 , . . . , ak ) ∈
Ck of (14) the polynomial r(vℓ , aℓ+1 , aℓ+2 , . . . , ak ) has the same degree in
vℓ as r(vℓ , vℓ+1 , . . . , vk ) and has no multiple roots. (Equivalently, the initial of r and the discriminant of r with respect to vℓ do not vanish on the
solution set of (14) in Ck .)
3. The differential consequences of p1 = 0, p2 = 0, . . . , ps = 0 contain all
integrability conditions of this PDE system, i.e., the cross-derivative of
each pair of distinct equations whose leaders involve the same unknown
function reduces to zero modulo p1 , . . . , ps and their derivatives (passivity
or formal integrability).
4. No reduction of q1 , q2 , . . . , qt is possible modulo p1 , p2 , . . . , ps and their
derivatives.
9

Definition 7. Let S be a system of partial differential equations and inequations, defined over R. A Thomas decomposition of S (with respect to ≻) is a
finite collection of simple differential systems S1 , . . . , Sr , defined over R, such
that the solution sets of S1 , . . . , Sr form a partition of the solution set of S.
The method outlined above allows to compute a Thomas decomposition
for any differential system S as considered above, with respect to any ranking
≻, in finitely many steps [15], [16]. However, a Thomas decomposition of a
differential system is not uniquely determined in general. The relevance of
simple differential systems and the decomposition of a general differential system
into simple differential systems is explained by the following theorem.
Theorem 8 ([16], Prop. 2.2.50). Let a simple differential system S be given by
p1 = 0 ,

... ,

ps = 0 ,

q1 6= 0 ,

... ,

qt 6= 0 ,

(s, t ∈ Z≥0 )

(15)

where p1 , . . . , ps , q1 , . . . , qt ∈ R. Let E be the differential ideal of R which is
generated by p1 , p2 , . . . , ps . Moreover, let q ∈ R be the product of the initials
and separants of p1 , p2 , . . . , ps . Then the differential ideal
E : q ∞ := { p ∈ R | q r p ∈ E for some r ∈ Z≥0 }
is equal to the set of differential polynomials in R which vanish under substitution of any analytic solution of S. In particular, it is a radical differential ideal.
A differential polynomial p is an element of E : q ∞ if and only if p reduces to
zero modulo p1 , . . . , ps and their derivatives.
Remark 9. Let S be a simple differential system as in Theorem 8. For each
k ∈ {1, . . . , m} let ∂ Jk,1 , ∂ Jk,2 , . . . , ∂ Jk,nk , where Jk,i ∈ (Z≥0 )n , nk ∈ Z≥0 , be
the differential operators such that ∂ Jk,i uk = ld(pji ) for some ji ∈ {1, . . . , s}.
Due to the characterization of the vanishing ideal given by Theorem 8, the set
of principal derivatives
P :=

nk
m [
[

{ ∂ J ∂ Jk,i uk | J ∈ (Z≥0 )n }

k=1 i=1

consists of those elements v ∈ Θu for which there exists an equation with leader
v that is a consequence of S. We refer to the elements of the complement Θu \ P
as the parametric derivatives. Note that ld(qj ) ∈ Θu \ P for all j = 1, . . . , t
because of conditions 1 and 4 of Definition 6.
The differential Thomas decomposition is an appropriate tool to study and
construct formal and convergent power series solutions (cf. [20, §1.4]).
For a formal statement, we need two, generically satisfied, regularity assumptions. Let S be a simple differential system as in Theorem 8 over a coefficient field of meromorphic functions in n complex variables x1 , . . . , xn . We call
x0 = (x01 , . . . , x0n ) ∈ Cn a regular point of S if all coefficients of pi and qj are
10

holomorphic and no init(pi ), sep(pi ), qj vanishes at x0 for all 1 ≤ i ≤ s, 1 ≤ j ≤ t.
Consider formal power series f = (f1 , . . . , fm ) given by
X
(xn − x0n )Jn
(x1 − x01 )J1
...
fk :=
ck,J
J1 !
Jn !
n
J∈(Z≥0 )

around a regular point x0 of S. We call any assignment ck,J ∈ C where ∂ J uk are
the parametric derivatives with qj (f1 , . . . , fm )(x01 , . . . , x0n ) 6= 0 for j = 1, . . . , t a
regular initial condition of S at the point x0 .
Theorem 10. Let S be a simple differential system and x0 a regular point of
S. Any regular initial condition of S at x0 completes to a formal power series
solution of S around x0 . For each such regular initial condition there are exactly
as many completions as the product of the degrees of the equations of S in their
leaders. In particular, if all equations in S are quasilinear, then there exists a
unique such completion.
There may exist additional formal power series solutions that stem from
non-regular initial conditions (cf. [20, §1.4], [21, Ex. 4.8, Ex. 4.9]).
By the results of Riquier [3], this theorem holds for convergent power series
solutions and, hence, generalizes Theorem 1 of Cauchy and Kovalevskaya (cf.
also [5], [22], [23], [24]). We call a regular initial condition of S at x0 C-regular,
if all fk are (locally) convergent when setting all ck,J to zero for all principal
derivatives ∂ J uk .
Corollary 11. Let S be a simple differential system w.r.t. a Riquier ranking
and x0 a regular point of S. If in Theorem 10 the regular initial condition is
C-regular, then any completion to a formal power series solution is also (locally)
convergent.
A ranking ≻ on R = Q{u1, . . . , um } is called an elimination ranking with
blocks B1 , B2 , . . . , Br if B1 ∪ B2 ∪ . . . ∪ Br is a partition of the set {u1 , . . . , um }
and we have
∂ K ui ≻ ∂ L uj ,
K, L ∈ (Z≥0 )n ,
whenever ui ∈ Bk and uj ∈ Bl with k < l.
By computing a Thomas decomposition with respect to an elimination ranking one may find all consequences of a differential system involving only differential indeterminates which are in lower ranked blocks. More precisely, we have
the following theorem.
Theorem 12 ([16], Prop. 3.1.36). Let S be a simple differential system with
respect to an elimination ranking ≻ with blocks B1 , B2 , . . . , Br . Let E and q
be defined as in Theorem 8, and for all i ∈ {1, . . . , r} let Ei be the differential
ideal of Q{ uj | uj ∈ B1 ∪ B2 ∪ . . . ∪ Bi } which is generated by the left hand
sides of the equations in S which only involve differential indeterminates in
B1 ∪ B2 ∪ . . . ∪ Bi , and let qi be the product of the initials and separants of these
differential polynomials. Then for every i ∈ {1, . . . , r} we have
(E : q ∞ ) ∩ Q{ uj | uj ∈ B1 ∪ B2 ∪ . . . ∪ Bi } = Ei : qi∞ .
11

4

Package commands

The ranking, the total order on the set of differential indeterminates, is globally determined by the command ComputeRanking(ivar,dvar). Here, ivar is
list of (lists of) independent variables, and dvar is a list of (lists of) unknown
functions. The unknown functions can be grouped in lists to use a block elimination ranking, e.g. dvar=[[u1],[u2,u3]] represents u1 and all its derivatives
being bigger than u2 or u3 or any of their derivatives. In case two unknown
functions are equally big, we can order their derivatives by grouping the independent variables in lists, e.g. ivar=[[x],[y,z]] means that a differential
indeterminate ranks higher if its order w.r.t. x is bigger; in case of a tie the oder
in y and z is considered. If no additional grouping is done, the degree-reverse
lexicographical ranking (cf. (9)) is chosen.
For the remainder of this section, we assume the following ranking.
ComputeRanking([x,y], [u,v]);
Derivatives of the unknown functions are encoded by indexing them with the
3
order vector as index, e.g. u[2, 1] ≡ ux,x,y ≡ ∂x∂2 ∂y u(x, y). The coefficients can
be any valid Maple expressions. In particular, functions in the independent
variables are correctly differentiated. Conversion from and to the diff notation
of Maple works via Diff2JetList and JetList2Diff.
p := Diff2JetList(u(x,y)+v(x,y)*diff(u(x,y),x,x,y));
p := v0,0 u2,1 + u0,0
JetList2Diff(p);

 3
u(x, y) + v(x, y) ∂x∂2 ∂y u(x, y)

We can compute (partial) derivatives, initial, leader, and separant of a differential polynomial.
PartialDerivative(p, x, y);
u2,1 v1,1 + u2,2 v1,0 + u3,1 v0,1 + u3,2 v0,0 + u1,1
Leader(p);
u2,1
Initial(p);
v0,0
Separant(p);
v0,0
The command DifferentialThomasDecomposition computes a Thomas decomposition. Its first parameter is a list of (left hand sides of) differential
equations and its second parameter is a list of (left hand sides of) differential
inequations. It returns a Thomas decomposition as a list of differential systems.

12

res:=DifferentialThomasDecomposition([p], []);
res := [DifferentialSystem, DifferentialSystem]
The commands DifferentialSystemEquations resp. DifferentialSystemInequations
resp. PrettyPrintDifferentialSystem extract equations resp. inequations
resp. both from a given differential system.
PrettyPrintDifferentialSystem(res[1]);
 3

[ ∂x∂2 ∂y u(x, y) v(x, y) + u(x, y) = 0, v(x, y) 6= 0 ]
PrettyPrintDifferentialSystem(res[2]);
[ u(x, y) = 0, v(x, y) = 0 ]
A differential system can be passed to the differential equation solvers of Maple
via MyPDSolve.
map(MyPDSolve, res);

[ {u(x, y) = 0, v(x, y) = v(x, y)} , u(x, y) = u(x, y), v(x, y) = −

u(x,y)
∂3
u(x,y)
∂y∂x2

{u(x, y) = 0, v(x, y) = 0} ]



,

The command IntersectDecompositions computes the intersection of two lists
of disjoint systems.
IntersectDecompositions(res,
DifferentialThomasDecomposition([u[0,0]-v[0,0]], []));
[DifferentialSystem, DifferentialSystem]
PrettyPrintDifferentialSystem(%);
[ [ u(x, y) − v(x, y) = 0,

∂3
∂y∂x2

v(x, y) + 1 = 0, v(x, y) 6= 0 ],

[ u(x, y) = 0, v (x, y) = 0 ] ]
Both DifferentialSystemReduce and DifferentialSystemNormalForm bring
(list of) differential equation(s) in a reduced form w.r.t. (list of) differential
system(s). This reduced form is zero iff the differential equation is a consequence of the system. The output is a (list of (lists of)) differential equation(s), as each differential equation is treated with each differential system.
DifferentialSystemReduce applies a differential reduction, hence the resulting
differential equation is equivalent to the input differential equation up to a factor, which is implied to be non-zero by the system. DifferentialSystemNormalForm
additionally divides by this factor to have an equivalent differential equation,
at the cost of denominators.
DifferentialSystemReduce(
[PartialDerivative(p,x),PartialDerivative(p,y),u[1,1]],
res);
13

[ [0, 0], [0, 0], [u1,1, 0] ]
The dimension polynomial DifferentialSystemDimensionPolynomial is a univariate polynomial function [25]. When evaluated on any high enough order,
it yields the number of generically freely choosable power series coefficients of
a solution of a simple differential system up to that order. It is a numerical
polynomial and, hence, can uniquely be represented as a linear combination of
i
the binomial coefficients Cs+i
for 0 ≤ i ≤ n with positive interger coefficients
by DifferentialSystemDimensionPolynomialCanonicalBase.
DifferentialSystemDimensionPolynomial(res[1]);
1 2
2s

+ 92 s + 1

DifferentialSystemDimensionPolynomialCanonicalBase(res[1]);

2+s
+ 3s
s

The command PowerSeriesSolution allows to compute the set of generic power
series solutions up to a given order around a given point.
PowerSeriesSolution(res[1], 2, [a,b]);
[u0,0 + u1,0 (x − a) + u0,1 (y − b) , v0,0 + v1,0 (x − a) + v0,1 (y − b)]

5

Examples

Example 13 (Consistency check). This example illustrates the verification of
consistency for differential systems. Consider the following systems of PDEs
with two independent variables x, y and the dependent variable u(x, y)
ux uy + ux + 1 = 0 ,

uux,x − ux − u2y + au = 0 ,

a ∈ R,

and detect all values of the parameter ’a’ that provide providing consistency
of the system. In doing so, we treat parameter ’a’ as a function of x, y with
identically zero partial derivatives.
with(DifferentialThomas):
ComputeRanking([x,y], [[u],[a]]):
DS := [u[1,0]*u[0,1]+u[1,0]+1,u[0,0]*u[2,0]-u[1,0]-u[0,1]^2+
a[0,0]*u[0,0],a[1,0],a[0,1]]:
TD := DifferentialThomasDecomposition([DS, []);
TD := [DifferentialSystem]
PrettyPrintDifferentialSystem(TD[1]);
3 
2


∂
∂
∂
[
u(x, y) +
u(x, y) + 2
u(x, y) + 1 = 0,
∂x
∂x
∂x

 


∂
∂
∂
u(x, y)
u(x, y) +
u(x, y) + 1 = 0, a(x, y) = 0, u(x, y) 6= 0 ]
∂y
∂x
∂x


14

Thus, the differential system under consideration is consistent if and only if
a = 0. Since every simple system has a solution [6, 15, 16], a PDE system is
consistent if and only if its Thomas decomposition is nonempty.
Example 14. As a second example we consider the Navier-Stokes equations
for an incompressible flow of a constant viscosity fluid:

ut + (u · ∇) u + ∇p − µ ∆u = 0 ,
(16)
∇ ·u = 0,
where u = (u, v, w) is the velocity vector, p the pressure, and µ = 1/Re for
the Reynolds number Re. In Cartesian coordinates x, y, z of R3 we have,
equivalently,

ut + u ux + v uy + w uz + px − µ (ux,x + uy,y + uz,z ) = 0 ,




vt + u vx + v vy + w vz + py − µ (vx,x + vy,y + vz,z ) = 0 ,
(17)

w
+
u wx + v wy + w wz + pz − µ (wx,x + wy,y + wz,z ) = 0 ,

t


ux + vy + wz = 0 .

We would like to verify algorithmically the well known fact that the Poisson
pressure equation
∆p + ∇ · ((u · ∇) u) = 0
(18)

is a consequence of the Navier-Stokes equations (16).
We choose the degree-reverse lexicographical ranking ≻ on the set of partial
derivatives


∂ i+j+k+l f
f
∈
{u,
v,
w,
p},
i,
j,
k,
l
∈
Z
≥0
∂ti ∂xj ∂y k ∂z l
which first compares the differential operators and in case of equality compares
the unknown functions to which the differential operators are applied according
to u ≻ v ≻ w ≻ p. For system (17) the resulting Thomas decomposition (with
respect to ≻) consists of only one simple differential system:














ux + vy + wz
µ vx,x + µ vy,y + µ vz,z − uvx − vvy − wvz − py − vt
µ uy,y + µ uz,z − µ vx,y − µ wx,z + uvy + uwz − vuy − wuz − px − ut
µ wx,x + µ wy,y + µ wz,z − uwx − vwy − wwz − pz − wt
px,x + py,y + pz,z + 2 uy vx + 2 uz wx + 2 vy2 + 2 vy wz + 2 vz wy + 2 wz2

=
=
=
=
=

0,
0,
0,
0,
0,

where some reductions have been performed on the first equation in (17) to
yield the second equation in the obtained simple system. The last equation was
obtained as


∂
∂
∂
∂
∂
∂
∂
A1 +
A2 +
A3 + v ∆ −
−u
−v
−w
+ 2 vy + 2 wz A4 ,
∂x
∂y
∂z
∂t
∂x
∂y
∂z
15

where A1 , A2 , A3 , A4 are the left hand sides of the equations in (17) and ∆ is
the Laplace operator. In fact, modulo the other equations in the system, the
last equation is equivalent to the Poisson pressure equation (18).
The above simple system of the Thomas decomposition allows to enumerate
the Taylor coefficients of u(t, x, y, z), v(t, x, y, z), w(t, x, y, z), p(t, x, y, z) whose
values can be chosen arbitrarily to obtain a power series solution of (16). These
correspond to the partial derivatives of u, v, w, p which are not highest ranked
derivatives in any differential equation among the consequences of the simple
system. For u(t, x, y, z) this enumeration is obtained from the non-zero terms
of the formal power series expansion of
1
,
(1 − ∂t ) (1 − ∂y ) (1 − ∂z )
where ∂t = ∂/∂t, ∂y = ∂/∂y, ∂z = ∂/∂z. In other words, an arbitrary analytic
function of t, y, z can be prescribed in terms of boundary conditions to determine
u(t, x, y, z). For v(t, x, y, z) the enumeration is given by
1
∂x
+
(1 − ∂t ) (1 − ∂y ) (1 − ∂z ) (1 − ∂t ) (1 − ∂z )
and for w(t, x, y, z) and for p(t, x, y, z) by
1
∂x
+
.
(1 − ∂t ) (1 − ∂y ) (1 − ∂z ) (1 − ∂t ) (1 − ∂y ) (1 − ∂z )
Hence, extending the Cauchy-Kovaleskaya Theorem 1 we may pose the Cauchy
problem for the Navier-Stokes equations (16) around an arbitrary point (t0 , x0 , y0 , z0 )
as follows:

u(t, x0 , y, z) = f1 (t, y, z) ,






v(t, x0 , y, z) = f2 (t, y, z) ,





∂v


(t, x0 , y0 , z) = f3 (t, z) ,



∂x



w(t, x0 , y, z) = f4 (t, y, z) ,


∂w


(t, x0 , y, z) = f5 (t, y, z) ,


 ∂x





p(t, x0 , y, z) = f6 (t, y, z) ,





∂p


(t, x0 , y, z) = f7 (t, y, z) ,
∂x
where f1 , f2 , . . . , f7 are arbitrary functions of their arguments which are analytic
around the point (t0 , x0 , y0 , z0 ). The arbitrariness of analytic solutions to (16)
is determined by f1 , f2 , . . . , f7 .
Example 15. We demonstrate how to study the classical Cole-Hopf transformation by using the differential Thomas decomposition (cf. also [15, Ex. 3.8]).
The claim is that for every non-zero analytic solution η(t, x) of the heat equation
ut + ux,x = 0
16

(19)

the function
ζ(t, x) :=

ηx (t, x)
η(t, x)

(20)

is a solution of Burgers’ equation
ut + ux,x + 2 u ux = 0 .

(21)

We define a ranking ≻ on the ring of differential polynomials in η and ζ such
that any partial derivative of η is ranked higher with respect to ≻ than any
partial derivative of ζ. According to Theorem 12 this allows determining all
differential equations which are satisfied by ζ from a Thomas decomposition
with respect to ≻.
with(DifferentialThomas):
ComputeRanking([t,x], [[eta],[zeta]]):
We define the differential system which combines (19) and (20).
CH := [eta[1,0]+eta[0,2], eta[0,0]*zeta[0,0]-eta[0,1]];
CH := [ η1,0 + η0,2 , η0,0 ζ0,0 − η0,1 ]
We also include the assumption η 6= 0 as an inequation.
TD := DifferentialThomasDecomposition(CH, [eta[0,0]]);
TD := [DifferentialSystem]
PrettyPrintDifferentialSystem(TD[1]);


∂
∂
∂
2
η(t, x) + (ζ(t, x)) η(t, x) +
ζ(t, x) η(t, x) = 0,
η(t, x) − ζ(t, x) η(t, x) = 0,
[
∂t
∂x
∂x
∂
∂
∂2
ζ(t, x) +
ζ(t, x) + 2 ζ(t, x)
ζ(t, x) = 0, η(t, x) 6= 0 ]
∂x2
∂t
∂x
The simple system of the resulting Thomas decomposition allows to read off
that ζ as defined by (20) is a solution of (21) if η is a solution of (19), which
proves the original claim. Conversely, since the above simple differential system
is consistent with the heat equation (19) for η by construction, we conclude
that for any solution ζ of (21) there exists a solution η of (19) such that the
Cole-Hopf transformation of η is ζ.
Example 16 (Detection of hidden constraints for the system with singular
Lagrangian). We apply the differential Thomas decomposition to compute the
hidden Lagrangian constraints for the problem taken from [26, Ex. 6.3], [27,
Eq. (8.1)]. The model Lagrangian reads
L = q22 (q1 )2t + q12 (q2 )2t + 2 q1 q2 (q1 )t (q2 )t + q12 + q22 .
We define L and the corresponding Euler-Lagrange equations in Maple:

17

L := q2[0]^2*q1[1]^2+q1[0]^2*q2[1]^2+2*q1[0]*q2[0]*q1[1]*q2[1]
+q1[0]^2+q2[0]^2;
EL := map(a->PartialDerivative(diff(L, a[1]), t)-diff(L, a[0]),
dvar);
L := q1 0 2 q2 1 2 + 2 q1 0 q1 1 q2 0 q2 1 + q1 1 2 q2 0 2 + q1 0 2 + q2 0 2
EL := [ 2 q1 0 q2 0 q2 2 + 4 q1 1 q2 0 q2 1 + 2 q1 2 q2 0 2 − 2 q1 0 ,
2 q1 0 2 q2 2 + 4 q1 0 q1 1 q2 1 + 2 q1 0 q1 2 q2 0 − 2 q2 0 ]
with(DifferentialThomas):
ivar := [t]: dvar := [q1,q2]:
We choose the ranking which compares the differentiation order of partial derivatives first and in case of equality gives priority to q1 over q2 .
ComputeRanking(ivar, dvar);
TD := DifferentialThomasDecomposition(EL, []);
TD := [DifferentialSystem, DifferentialSystem, DifferentialSystem]
PrettyPrintDifferentialSystem(TD[1]);
 2


2
d
d
[ q1 (t) + q2 (t) = 0, 2 q2 (t)
q2 (t) + 2
q2 (t) − 1 = 0, q2 (t) 6= 0 ]
dt2
dt
The first expression q1 (t) + q2 (t) = 0 is a new constraint in the simple system
TD[1].
PrettyPrintDifferentialSystem(TD[2]);


2
 2
d
d
q2 (t) + 2
q2 (t) − 1 = 0, q2 (t) 6= 0 ]
[q1 (t) − q2 (t) = 0, 2 q2 (t)
dt2
dt
The first expression q1 (t) − q2 (t) = 0 is a new constraint in the simple system
TD[2].
PrettyPrintDifferentialSystem(TD[3]);
[ q1 (t) = 0, q2 (t) = 0 ]
These constraints are complementary to the first two constraints.
Example 17. The Landau-Lifshitz-Gilbert equations describe a unit magnetization vector


u(t)
m :=  v(t) 
w(t)
m:=<u[0],v[0],w[0]>:
with three dependent variables u, v, and w as entries. We automatically
replicate results from [28]. Therefore, we assume a constant external magnetic field
18

h_eff:=<0,0,h3-lambda*m[3]>:
in the direction of the w-axis, a self interaction of the magnetic vector in
the same direction for λ ∈ {−1, 0, 1}, and assume an additional spin torque
term
j:=<0,0,j3>:
aligned in direction to the w-axis, hence counteracting the damping.
Under these assumptions the Landau-Lifshitz-Gilbert equations are given by
the equations

α2 + 1 mt = αm × (m × heff ) − αm × j − m × heff − m × (m × j)
LLG:=(alpha^2+1)*<u[1],v[1],w[1]>
-alpha*CrossProduct(m,CrossProduct(m,h_eff))
+alpha*CrossProduct(m,j)
+CrossProduct(m,h_eff)
+CrossProduct(m,CrossProduct(m,j)):

where α is a positive real number indicating the strength of the damping.
The alignedness assumption implies symmetry w.r.t. rotations around the waxis. Hence, all periodic solutions must be parallel to the equator, which implies
d
dt w(t) = 0. We furthermore encode the unit length of the magnetization vector
and model the parameters by differential indeterminates with derivative zero.
LLG:=[op(convert(LLG,list)),w[1],Diff2JetList(u(t)^2+v(t)^2+w(t)^2-1),
alpha[1],j3[1],h3[1],lambda[1]]:
We use the inequation α 6= 0 to model 0 < α. To remove superfluous complex
solutions, add the inequation α2 +1 6= 0. Furthermore, assume that the external
field and the spin torque are non-degenerate, i.e., that h3 and j3 are non-zero.
ineq:=[alpha[0],alpha[0]^2+1,h3[0],j3[0]]:
A Thomas decomposition of this system consists of 24 simple systems:
ComputeRanking([t],[u,v,w,alpha,j3,h3,lambda]);
res:=DifferentialThomasDecomposition(LLG,ineq):
nops(res);
24
However, many of these systems contain only non-real (complex) solutions. Remove the systems that imply one of the following four pairs of constraints:
u(t)2 + v(t)2 = 0 6= v(t)
j32 + (h3 − λ)2 = 0 6= h3 − λ
j32 + (h3 + λ)2 = 0 6= h3 + λ
j32 + 4λ2 = 0 6= λ

19

l:=[[u[0]^2+v[0]^2,v[0]],
[j3[0]^2+(h3[0]-lambda[0])^2,h3[0]-lambda[0]],
[j3[0]^2+(h3[0]+lambda[0])^2,h3[0]+lambda[0]],
[j3[0]^2+4*lambda[0]^2,lambda[0]]]:
for c in l do
res:=remove(a->DifferentialSystemReduce(a,c[1])=0
and c[2] in DifferentialSystemInequations(a), res):
od:
nops(res);
6
The first system describes the generic case. It results from incorporating
some inequations and transforming the equations. Hence, the genericly correct
information can be read off this system. First we observe the equation
PrettyPrintDifferentialSystem(res[1])[3];
GenericEquation:=Diff2JetList(lhs(%)):
−α (t) λ (t) w (t) + α (t) h3 (t) − j3 (t) = 0
This equation holds whenenever w(t) is not the constant function 1 or −1:
map(b->DifferentialSystemReduce(b,w),
select(a->0<>DifferentialSystemReduce(a,GenericEquation),res));
[−1, 1]
We describe the real solutions
with(RealDomain):
of the first system and make some assumptions about signs of trigonometric
functions to get a simplified form. Making the opposite assumption only changes
signs in the solution set.
sol:=MyPDSolve(res[1]):
l:=[subs(sol,alpha(t))=alpha,subs(sol,lambda(t))=lambda,
subs(sol,j3(t))=j3,subs(sol,h3(t))=h3]:
u(t)=simplify(factor(subs(l,subs(sol,u(t)))),trig) assuming alpha>0:
simplify(%) assuming cos(j3*signum(lambda)*(-t+_C5)/alpha)>0:
subs([-j3+alpha*h3+alpha*lambda=c1,-j3+alpha*h3-alpha*lambda=c2],%):
simplify(subs([c2=c3/c1],%)):
subs(c3=(-j3+alpha*h3+alpha*lambda)*(-j3+alpha*h3-alpha*lambda),%):
subs(signum(lambda)=lambda,%);


C5 −t)
((h3 − λ) · α − j3)((h3 + λ) · α − j3) cos j3λ( α
p
u(t) = −
α|λ| −((h3 − λ) · α − j3)((h3 + λ) · α − j3)
v(t)=simplify(factor(subs(l,subs(sol,v(t)))),trig) assuming alpha>0:
simplify(%) assuming cos(j3*signum(lambda)*(-t+_C5)/alpha)>0:
subs(signum(lambda)=lambda,%);


p
j3 λ ( C5 − t)
−1
(α · λ)
v(t) = − −((h3 − λ) · α − j3)((h3 + λ) · α − j3) sin
α
20

w(t)=subs(l,subs(sol,w(t)));
w(t) =

αh3 − j3
αλ

In this case u and v are phase shifted, and the angular velocity is ± j3
α , as
λ ∈ {−1, 0, 1}. The solutions are real if and only if
factor(subs([j3=j3_divided_by_alpha*alpha,lambda=1],
-(-j3+alpha*h3-alpha*lambda)*(-j3+alpha*h3+alpha*lambda)))>0:
solve(subs(alpha=1,%),j3_divided_by_alpha) for j3_divided_by_alpha>0:
subs(j3_divided_by_alpha=j3/alpha,%);


j3
j3
< h3 + 1, h3 − 1 <
α
α
holds. (This last computation yields the same results for λ = −1.)
The solutions in the second, third, resp. fourth system are the same, with
the additional constraints h3 = λ, h3 = −λ, resp. λ = 0. The last two systems
contain the constant solution at the south pole and north pole.

6

Conclusion

In this paper we presented an overview of the differential Thomas decomposition
method and of an implemention in Maple. (Systems of) differential equations of
polynomial type arise everywhere in physics. The differential Thomas decomposition provides a universal algorithmic tool for investigating their algebraic
properties and for constructing solutions by splitting the original equation(s)
into a finite set of differential subsystems with disjoint solution sets. These
new systems have additional algebraic properties, e.g. a triangular structure,
the absence of multiple solutions, they disclose hidden constraints, they generate explicit conditions on parameters, and they detecte the arbitrariness in
analytic solutions. The Thomas decomposition is applicable both to a single
polynomially nonlinear differential equation and to a system of such equations.
Moreover, one can also add inequation(s), which removes unwanted solutions.
We illustrated these properties of the output of differential Thomas decomposition by a number of examples (Sect.5). All these examples, except the first
one, arise in the context of mathematical physics.
It should be noted that for an input nonlinear differential system with several independent and dependent variables, especially with a high degree of nonlinearity, the intermediate computations can be very tedious and/or require a
large amount of computer memory. More often than not, this happens when an
elimination ranking is used.

7

Acknowledgments

The contribution of the first author (V.P.G.) has been partially supported
by the Russian Foundation for Basic Research (grant No.16-01-00080) and by
21

the Ministry of Education and Science of the Russian Federation (agreement
02.a03.21.0008)

References
[1] S. Kowalevsky, Zur Theorie der partiellen Differentialgleichung, J. Reine &
Angewandte Mathematik 80 (1875) 1–32.
[2] I. G. Petrovsky, Lectures on partial differential equations, Interscience Publishers, New York-London, 1954, translated by A. Shenitzer.
[3] C. Riquier, Les systèmes d’équations aux dérivées partielles, GauthierVillars, Paris, 1910.
[4] M. Janet, Les systèmes d’équations aux dérivées partielles, J. de Math., 8e
série, 3 (1920) 65–151.
[5] J. M. Thomas, Riquier’s existence theorems, Ann. of Math. (2) 30 (1-4)
(1928/29) 285–310.
[6] J. M. Thomas, Differential Systems, American Mathematical Society Colloquium Publications, Vol. XXI, American Mathematical Society, New York,
N. Y., 1937.
[7] C. J. Rust, G. J. Reid, Ranking of Partial Derivatives, in: W. W. Küchlin
(Ed.), Proceedings of the 1997 International Symposium on Symbolic and
Algebraic Computation, Maui, Hawaii (USA), ACM Press, 1997, pp. 9–16.
[8] J. M. Thomas, Systems and Roots, William Byrd Press, Richmond, VA,
1962.
[9] J. F. Ritt, Differential Algebra, American Mathematical Society Colloquium Publications, Vol. XXXIII, American Mathematical Society, New
York, N. Y., 1950.

[10] W. T. Wu, A constructive theory of differential algebraic geometry based on works of J. F. Ritt with part
in: Differential geometry and differential equations (Shanghai, 1985), Vol.
1255 of Lecture Notes in Math., Springer, Berlin, 1987, pp. 173–189.
URL https://doi.org/10.1007/BFb0077689
[11] Z. Li, D. Wang, Coherent, regular and simple systems in zero decompositions of partial differential systems, System Science and Mathematical
Sciences 12 (1999) 43–60.
[12] D. Wang, Elimination methods, Texts and Monographs in Symbolic Computation, Springer-Verlag, Vienna, 2001.
[13] D. Wang, Elimination practice, Imperial College Press, London, 2004, software tools and applications. With 1 CD-ROM (UNIX/LINUX, Windows).
http://www-salsa.lip6.fr/~wang/epsilon.
22

[14] V. P. Gerdt, On decomposition of algebraic PDE systems into simple subsystems, Acta Appl. Math. 101 (1-3) (2008) 39–51.
[15] T. Bächler, V. P. Gerdt, M. Lange-Hegermann, D. Robertz, Algorithmic
Thomas decomposition of algebraic and differential systems, J. Symbolic
Comput. 47 (10) (2012) 1233–1266.
[16] D. Robertz, Formal Algorithmic Elimination for PDEs, Vol. 2121 of Lecture
Notes in Mathematics, Springer, Cham, 2014.
[17] F. Boulier, D. Lazard, F. Ollivier, M. Petitot, Computing representations
for radicals of finitely generated differential ideals, Appl. Algebra Engrg.
Comm. Comput. 20 (1) (2009) 73–121.
[18] G. J. Reid, A. D. Wittkopf, A. Boulton, Reduction of systems of nonlinear
partial differential equations to simplified involutive forms, European J.
Appl. Math. 7 (6) (1996) 635–666.
[19] B. Sturmfels, What is . . . a Gröbner basis?, Notices Amer. Math. Soc.
52 (10) (2005) 1199–1200.
[20] M. Lange-Hegermann, Counting Solutions of Differential Equations,
Ph.D. thesis, RWTH Aachen University, Germany, available online at
http://publications.rwth-aachen.de/record/229056 (2014).
[21] M.
Lange-Hegermann,
The
differential counting polynomial,
Foundations of Computational Mathematics, available online at
https://doi.org/10.1007/s10208-018-9378-z,
arXiv:1407.5838
(2018).
[22] J. M. Thomas, Riquier’s existence theorems, Ann. of Math. (2) 35 (2)
(1934) 306–311.
[23] V. P. Gerdt, Algebraically simple involutive differential systems and Cauchy
problem, Zap. Nauchn. Sem. S.-Peterburg. Otdel. Mat. Inst. Steklov.
(POMI) 373 (Teoriya Predstavleniı̆, Dinamicheskie Sistemy, Kombinatornye Metody. XVII) (2009) 94–103, 347.
[24] C. J. Rust, G. J. Reid, A. D. Wittkopf, Existence and uniqueness theorems
for formal power series solutions of analytic differential systems, in: Proceedings of the 1999 International Symposium on Symbolic and Algebraic
Computation (Vancouver, BC), ACM, New York, 1999, pp. 105–112.
[25] M. Lange-Hegermann, The differential dimension polynomial for characterizable differential ideals, in: Algorithmic and Experimental Methods in
Algebra, Geometry, and Number Theory, Springer, 2017, accepted for publication, arXiv:1401.5959.
[26] V. P. Gerdt, D. Robertz, Lagrangian Constraints and Differential Thomas
Decomposition, Advances in Applied Mathematics 72 (2016) 113–138.
23

[27] A. Deriglazov, Classical mechanics, Hamiltonian and Lagrangian formalism, Springer, Heidelberg, 2010.
[28] F.
Müller,
Hopf
bifurcation
in
Landau-Lifshitz-Gilbert
dynamics
with
a
spin
torque,
Available
online
at
https://opus4.kobv.de/opus4-matheon/frontdoor/index/index/docId/402
(2007).

24

