arXiv:1401.5438v1 [math.CA] 21 Jan 2014

On the Reduction of Singularly-Perturbed
Linear Differential Systems∗
Moulay Barkatou, Suzy S. Maddah
XLIM UMR 7252 ; DMI
University of Limoges; CNRS
123, Avenue Albert Thomas
87060 Limoges, France
moulay.barkatou@unilim.fr
suzy.maddah@etu.unilim.fr

†

Hassan Abbas
Laboratory of Mathematics
Lebanese University, Beirut, Lebanon
habbas@ul.edu.lb
February 17, 2018

Abstract
In this article, we recover singularly-perturbed linear differential
systems from their turning points and reduce their parameter singularity’s rank to its minimal integer value. Our treatment is Moser-based;
that is to say it is based on the reduction criterion introduced for linear
singular differential systems in [Moser, 1960]. Such algorithms have
proved their utility in the symbolic resolution of the systems of linear functional equations [Barkatou, 1989; Barkatou et al., 2008, 2009],
giving rise to the package ISOLDE [Barkatou et al., 2013], as well as in
the perturbed algebraic eigenvalue problem [Jeannerod et al., 1999].
In particular, we generalize the Moser-based algorithm described in
[Barkatou, 1995]. Our algorithm, implemented in the computer algebra system Maple, paves the way for efficient symbolic resolution of
singularly-perturbed linear differential systems as well as further applications of Moser-based reduction over bivariate (differential) fields
[Barkatou et al., 2014].
∗
†

Submitted to ISSAC’14, Kobe, Japan
Enrolled under a joint PhD program with the Lebanese University

1

Keywords: Moser-based reduction, Perturbed linear differential systems,
Turning points, Computer algebra.

1

Introduction

Let K be a commutative field of characteristic zero equipped with a derivation δ, that is a map δ : K → K satisfying
δ(f + g) = δf + δg and δ(f g) = δ(f )g + f δ(g) for all f, g ∈ K.
We denote by C its field of constants, V a K-vector space of dimension n, and
A an element in K n×n . Set ∆ = δ − A, then ∆ is a δ- differential operator
acting on V , that is, a C-linear endomorphism of V satisfying the Leibniz
condition:
∀f ∈ K, v ∈ V ∆(f v) = δ(f )v + f ∆(v).
Let O = C[[x, ǫ]], the ring of formal power series in x and ǫ over the field of
complex numbers, where x is a complex variable and ǫ is a small parameter.
d
Denote by K its field of fractions equipped with δ = ǫ dx
. Let Ox = C[[x]],
Kx be its fields of fractions, and h, p be integers (h ≥ 0). Thereby, ∆ is the
differential operator associated to a singularly-perturbed linear differential
system. Denoting by Y an unknown n-dimensional column vector, the former
is written as
∞
X
dY
−h −p
Ak (x)ǫk Y.
(1)
= A(x, ǫ)Y = ǫ x
ǫ
dx
k=0
Such systems have countless applications which are traced back to the year
1817 and their study encompasses a vast body of literature (see, e.g. [Chen,
1984; Lin, 1966; McHugh, 1971; Wasow, 1985] and references therein). However, their symbolic resolution is still open to investigation.
Clearly, system (1) is a singular perturbation of the widely studied linear
singular system of differential equations (see, e.g., [Balser, 2000; Hsieh et al.,
d
1999] and references therein). The latter is obtained by setting δ = x dx
and
K = C((x)), the univariate field of formal Laurent power series. Hence ∆ is
associated to
∞
X
dY
−p
x
Ak xk Y.
(2)
= A(x)Y = x
dx
k=0
Contrary to system (1), there exist efficient algorithms contributing to the
formal reduction of system (2) (see, e.g [Barkatou, 1997; Pflugel , 2000]),
that is the algorithmic procedure that computes a change of basis w.r.t.
which A(x) has normal form facilitating the construction of formal solutions.
Among which, rank reduction algorithms play an eminent role as they reduce
the nonnegative integer p, called the Poincaré rank, to its minimal integer
value, the true Poincaré rank (see, e.g., [Barkatou et al., 2009; Levelt, 1991]).
In particular, if the latter is null then system (2) is said to be regular singular,
2

that is to say, in any small sector, its solutions grow at most as an algebraic
function. Otherwise, it is irregular singular.
However, the classical formal simplification of system (2) begins with the
reduction of its leading coefficient matrix A0 to its Jordan form. Hence, in
addition to the usual difficulties encountered within the formal reduction itself, additional ones arise for system (1) since its leading coefficient matrix
A0 (x) is a matrix-valued function rather than a constant one. In particular, if it has mulitple eigenvalues then they might coalesce (see, e.g., the
example page 223, [Wasow, 1985]). In classical textbooks, points where the
Jordan form of A0 (x) is unstable, that is to say either the multiplicity of the
eigenvalues or the degrees of the elementary divisors are not constant in a
neighborhood of such points, are referred to as turning points (see, e.g., page
57, [Wasow, 1979]). The behavior of solutions of differential systems around
such points is far more complicated than that of system (2) around an irregular singularity. In fact, the neighborhood of such a point is decomposed into
a finite number of x-dependent neighborhoods in each of which the solution
behaves quite differently though it may still be asymptotically accomplished
(see, e.g., [Iwano et al., 1963]). In particular, if h ≤ −1 then the solution of
system (1) can be sought upon presenting the solution as a power series in
ǫ. The latter can then be inserted into (1) and the like powers of ǫ equated.
This reduces the problem to solving recursively a set of non-homogeneous
linear singular differential systems over Kx , the first of which is system (2)
(see, e.g. page 52, [Chen, 1984] and page 101, [Wasow, 1979]). Hence, it
seems plausible to investigate wether h rather than p can be diminished for
such systems.
The methods proposed in the literature of system (1) either exclude turning point [Chen, 1984] or are not algorithmic throughout. Moreover, they
make an essential use of the so-called Arnold-Wasow form. For the univariate
case, system (2), the research advanced profoundly in the last two decades
making use of methods of modern algebra and topology. The former classical approach is substituted by efficient algorithms (see, e.g., [Barkatou, 1997;
Barkatou et al., 2013, 2009; Pflugel , 2000] and references therein). It was
the hope of Wasow [Wasow, 1985], in his 1985 treatise summing up contemporary research directions and results on system (1), that techniques of
system (2) be generalized to tackle the problems of system (1).
In this article, we are interested in recovering system (1) from its turning
points and decoupling it into a set of systems of lower dimensions for each of
which h has a minimal integer value.
Given system (2), Moser defined two rational numbers:
m(A) =
µ(A) =

rank(A0 )
) and
n
min { m(T −1 ∆T ) | T ∈ GL(V )}.
max (0, p +

(3)
(4)

It follows that system (2) is regular whenever µ(A) ≤ 1. For m(A) > 1, he
3

proved that m(A) > µ(A) if and only if the polynomial
θ(λ) := xrank(A0 ) det(λI +

A0
+ A1 )|x=0
x

vanishes identically in λ. In this case, system (2) (resp. A(x)) is said to
be Moser-reducible and m(A) can be diminished by applying a coordinate
transformation Y = T Z where T ∈ GL(V ) of the form
T = (P0 + P1 x) diag(1, . . . , 1, x, . . . , x)
where P0 , P1 are constant matrices and det(P0 ) 6= 0 [Theorems 1 and 2,
page 381, [Moser, 1960]]. This notion and algorithms developed (see, e.g.
[Barkatou et al., 2009]), were generalized as well to linear functional matrix
equations in [Barkatou et al., 2008], a particular case of which is system (2).
This gave rise to the package ISOLDE [Barkatou et al., 2013], written in the
computer algebra system Maple and dedicated to their symbolic resolution.
Moser’s reduction criterion was also borrowed from the theory of differential systems to investigate efficient algorithmic resolution of the perturbed
algebraic eigenvalue-eigenvector problem in [Jeannerod et al., 1999]. This
problem is another prominent example in the univariate case, for which δ is
the zero map and K = C((ǫ)). Thereby, ∆ is just a linear operator in the
standard way and A(ǫ) is the widely studied perturbation of the constant
matrix A0 (see, e.g., [Kato, 1980]).
However, despite their utility and efficiency in these univariate cases,
Moser-based algorithms are not considered yet over bivariate fields. In particular, the second author of this article developed in [Barkatou, 1995] a
d
Moser-based algorithm for the differential systems associated to δ = x dx
and
K = Q(x), the field of rational functions in x. This algorithm is the one we
generalize here to system (1). Although Moser’s motivation for introducing
these rational numbers for system (2) was to distinguish regular singular systems in the sense of Poincaré, we show hereby that they serve as well for the
reduction of h to its minimal integer value.
This article is organized as follows: In Section 2, we give preliminaries. In
Section 3, we present a Moser-based method of recovery from turning points.
In Section 4, we introduce a Moser-based reduction algorithm to diminish h
to its minimal integer value. Our main results are Theorem 1 and Proposition 1. In Subection 4.3, we give an outline of a generalization of Levelt’s
rank reduction algorithm of [Levelt, 1991] and illustrate by an example its
comparison to the Moser-based one. Finally, we conclude in Section 5 and
point out some prospects of further research. We remark however that as
our interest is formal, reference to asymptotic behavior and/or convergence is
dropped. One may consult in the latter direction, e.g., [Balser, 2000; Wasow,
1985].

4

2

Preliminaries

Consider system (1). Without loss of generality, we can assume that Ak (x) ∈
Kxn×n and the leading coefficient matrix A0 (x) is a nonzero element of Oxn×n
otherwise h and p can be readjusted. We refer to A00 := A(0, 0) as the
leading constant matrix. The nonnegative integer h will be referred to as
the rank of singularity in ǫ (ǫ-rank, in short). We assume that system (1)
has at most one turning point otherwise the region of study can be shrinked.
Moreover, this turning point is placed at the origin, otherwise, a translation
in the independent variable can be performed. We denote by Ik and Ok,l the
identity and zero matrices respectively of prescribed dimensions.
Let T ∈ GLn (K) then the change of the dependent variable Y = T Z in
(1) gives
∞
X
dZ
Ak (x)ǫk Z.
(5)
ǫ
= Ã(x, ǫ)Z = ǫ−h̃ x−p̃
dx
k=0
where Ã(x, ǫ) ∈ K n×n and
T [A] := Ã(x, ǫ) = T −1 A(x, ǫ)T − ǫ T −1

d
T.
dx

(6)

Systems (1) and (5) are called equivalent. In fact, T is a change of basis and
Ã(x, ǫ) is the matrix of T −1 ∆T .
Remark 1 Given system (1) and an equivalent system (5).
• If h > 0 and T (x) ∈ GLn (Kx ) then it follows from (6)
x−p̃ (Ã0 (x) + Ã1 (x)ǫ) = T −1 x−p (A0 (x) + A1 (x)ǫ)T.
Hence, since we’ll be interested solely in Ã0 (x) and/or Ã1 (x), it suffices
to investigate T −1 AT , the similarity term of T [A]. In particular, if
T (x) ∈ GLn (Ox ) then p = p̃ as well.
• If T = diag(ǫα1 , . . . , ǫαn ), where α1 , . . . , αn are integers, then δT is a
zero matrix. Hence, T [A] = T −1 AT .
A classical tool in perturbation theory is the so-called splitting which separates off the existing distinct coalescence patterns. Whenever the leading
constant matrix A00 admits at least two distinct eigenvalues, there exists
a coordinate transformation T ∈ GLn (O) which block-diagonalizes System
(1) so that it can be decoupled into subsystems of lower dimensions whose
leading constant matrices have a unique distinct eigenvalue each (see, e.g.,
Theorem XII − 4 − 1 page 381, [Hsieh et al., 1999]).
Remark 2 In (Theorem 2.3 − 1, page 17, [Wasow, 1985]), we have a generalization of the splitting to the well-behaved case which refers to A0 (x) being
5

holomorphically similar to its holomorphic Jordan form. In particular, its
eigenvalues don’t coalesce in the region of study. However, we limit ourselves
to the weaker version given above as we aim to give a general discussion
which doesn’t exclude turning points.
We now consider one of these resulting subsystems and assume its leading
constant matrix is in jordan normal form with a unique eigenvalue γ ∈ C.
This can be always attained by a constant transformation. Upon applying
the so-called eigenvalue shifting, i.e.
Z
−h−1
γx−p dx)Z,
Y = exp(ǫ
it is easy to verify from (6) that the resulting system has a nilpotent leading
constant matrix. Hence, without loss of generality, we can assume that system (1) is such that A00 is nilpotent. Clearly, it doesn’t follow that A0 (x) is
nilpotent and we deviate here from the classical treatment of system (2) as
we may encounter turning points.
The subscripts are to be ommitted and x, ǫ dropped from the notation
whenever ambiguity is not likely to arise, e.g. A(x, ǫ) , A0 (x), A1 (x) will be
denoted by A, A0 , A1 respectively. We set r = rank(A0 (x)).

3

Recovery from Turning Points

We arrive at this section with A0 (x) of system (1) being nonnilpotent in
contrary to the leading constant matrix A00 . We show that by Moser-based
reduction of A0 (x) itself and ramification in x, that is a readjustment x = ts
with s a positive integer, we can modify radically the nilpotency of A00 , i.e.
arrive at a system for which splitting can be applied.
The motivation behind considering such a general form of system (1)
rather than that for which p ≤ 0 is to be justified in this section. In fact,
system (1) will undergo a sequence of transformations which might introduce
or elevate the order of the pole in x. This elevation is inevitable and is
introduced identically in the classical treatment of such systems. However,
the order of poles of Ak (x), k = 0, 1, 2, . . . grows at worst linearly with k
which maintains an asymptotic validity of the formal construction [page 60,
[Wasow, 1979]].
By Remark (1), the discussion is restricted to the similarity term T −1 AT
of T [A]. Hence, A0 (x) can be treated as a perturbation of A00 .
Proposition 1 Given system (1)
ǫ

∞
X
dY
= A(x, ǫ)Y = ǫ−h x−p
Ak (x)ǫk Y
dx
k=0

6

where A0 (x) = A00 + A01 x + A02 x2 + . . . such that A00 is a nilpotent constant
matrix. If A0 (x) has a nonzero eigenvalue then there exist a positive integer
s and a T ∈ GLn (Kx ) such that by setting x = ts , Y = T Z, we have in (5)
Ã0 (t) := T −1 A0 (t)T = Ã00 + Ã01 t + Ã02 t2 + . . .
where Ã00 has a nonzero constant eigenvalue.
Proof. The eigenvalues of A0 (x) admit a formal expansion in the fractional
powers of x in the neighborhood of x = 0 (see, e.g. [Kato,P1980]). We are
j/s
interested only in their first nonzero terms. Let µ(x) = ∞
be a
j=0 µj x
nonzero eigenvalue of A0 (x) whose leading exponent, i.e. smallest j/s for
which µj 6= 0 and j, s are coprime, is minimal among the other nonzero
eigenvalues. Then, without loss of generality, we can assume that s = 1
otherwise we set x = ts . Now, let T ∈ GLn (Kx ) such that Ã0 (x) is Moserirreducible in x i.e.
θ(λ) = xrank(Ã00 ) det(λI +

Ã00
+ Ã01 )|x=0
x

doesn’t vanish identically in λ. Then there are n − degθ eigenvalues of A0 (x)
whose leading exponents lie in [0, 1[, and deg θ eigenvalues for which the
leading exponent is greater or equal 1 [?]in [Jeannerod et al., 1999]]. But
µ(x) is an eigenvalue of Ã0 (x) with a minimal leading exponent and hence
it is among those whose leading exponents lie in [0, 1[. By our assumption,
s = 1 and hence the leading exponent of µ(x) is zero and µ0 6= 0. Since µ0 is
an eigenvalue of Ã00 , it follows that the latter is nonnilpotent.
Remark 3 The eigenvalues of A0 (x) are the roots of the algebraic scalar
equation f (x, µ) = det(A(x) − µIn ) = 0 and can be computed by NewtonPuiseux algorithm. The linear transformation T ∈ GLn (Kx ) can be computed
efficiently via ISOLDE [Barkatou et al., 2013].
We illustrate our approach by an Example from [page 88, [Wasow, 1979]].
We remark however that the technique proposed in the former, adapted
from [Turrittin, 1952] and [Iwano et al., 1963], debutes with the reduction
of A(x, ǫ) to its Arnold-Wasow form as mentioned above and then constructing the so-called characteristic polygon. This particular example was given
initially in that for simplicity and hence, the transformations computed by
our algorithm, which doesn’t require reduction to this form, do not deviate
hereby from those in the former.


0 1 0
Example 1 Let ǫ dY
= ǫ−1 0 0 1  Y . Clearly, we have a truning point
dx
ǫ 0 x
at x = 0 and


0 1 0
A0 (x) = 0 0 1  .
0 0 x
7

A00 is nilpotent in Jordan form while the eigenvalues of A0 (x) are 0, 0, x
whence s = 1. A simple calculation shows that A0 (x) is Moser-reducible in
x (θ(λ) ≡ 0). Let T = diag(1, x, x2 ) then upon setting Y = T Z we get

 

0 1 0
0 0 0
dZ
= ǫ−1 x{0 0 1 + 0 0 0 x−3 ǫ +
ǫ
dx
0 0 1
1 0 0


0 0
0
0 −1 0  x−2 ǫ2 }Z.
0 0 −2
The constant leading coefficient matrix of the resulting matrix is no longer
nilpotent. Hence the system can be decoupled into two subsystems upon setting
Z = T W where

 

1 0 1
−1 −1 0
T = 0 1 1 + −1 −1 −2 x−3 ǫ + O(x−6 ǫ2 ).
0 0 1
−1 −1 0
The resulting equivalent system then consists of the two uncoupled lower dimension systems where W = [W1 , W2 ]T :

 

dW1
0 1
−1 0 −3
−1
ǫ
= ǫ x{
+
x ǫ+
0 0
−1 0
dx


1
−1
x−6 ǫ2 + O(x−9 ǫ3 )}W1 .
1 −1 + x4
dW2
= ǫ−1 x{1 + x−3 ǫ + (1 + 2x4 )x−6 ǫ2 + O(x−9 ǫ3 )}W2 .
dx
For the former subsystem, A0 (x) and A00 are now simultaneously nilpotent.
As for the latter, the exponential part is 21 ǫ−2 x2 (1 + O(x−3 ǫ)) which exhibits
that the eigenvalue x of the leading coefficient matrix of the system we started
with is recovered as expected.
ǫ

4

ǫ-Rank Reduction

Following Section 3, we can now assume without loss of generality that A0 (x)
and A00 are simultaneously nilpotent. The system is recovered from its turning points and no further reduction can be attained via Splitting. In analogy
to modern techniques for system (2), we investigate ǫ-rank reduction of system (1).

8

4.1

ǫ-Reduction Criteria

Following (3) and (4), we define respectively the ǫ-Moser rank and ǫ-Moser
Invariant of system (1) as the rational numbers:
mǫ (A) =
µǫ (A) =

r
) and
n
min {mǫ (T [A])|T ∈ GLn (K) }.
max (0, h +

Definition 1 System (1) (the matrix A respectively) is called ǫ-reducible in
the sense of Moser (ǫ-reducible, in short) if mǫ (A) > µǫ (A), otherwise it is
said to be ǫ-irreducible.
Remark 4 This definition is not to be mixed neither with the usual sense
of reduced system in the literature, i.e. the system (2) obtained from system
(1) (ǫ = 0); nor with the usual sense of Moser-reduced systems as for system
(2), i.e. the systems whose p is minimal.
In particular, if mǫ (A) ≤ 1 then h = 0 and no further reduction is to be
sought via this approach.

4.2

ǫ-Reduction Algorithm

We follow the algorithmic description of [Barkatou, 1995]. The main result
of this section is the following theorem which is to be proved after giving its
useful building blocks.
Theorem 1 Given system (1)
∞
X
dY
−h −p
ǫ
Ak (x)ǫk Y
= A(x, ǫ)Y = ǫ x
dx
k=0

such that rank(A0 (x)) = r and mǫ (A) > 1. A necessary and sufficient
condition for the system to be ǫ-reducible (in the sense of Moser), i.e. the
existence of a T (x, ǫ) ∈ Gln (K) such that r(Ã0 (x)) < r, is that the polynomial
θ(λ) := ǫr det(λI +

A0
+ A1 )|ǫ = 0
ǫ

vanishes identically in λ.
Moreover, T (x, ǫ) can always be chosen to be a product of unimodular transformations in GLn (Ox ) and polynomial transformations of the form diag(ǫα1 , . . . , ǫαn )
where α1 , . . . , αn are nonnegative integers.

9

Lemma 1 There exists a unimodular matrix U(x) in GLn (Ox ) such that the
leading coefficient matrix Ã0 (x) of U[A] have the following form
 11

Ã0 O
Ã0 =
(7)
Ã21
O
0
 11 
Ã0
11
where Ã0 is a square matrix of dimension r and
is a n × r matrix of
Ã21
0
full column rank r.
Proof. By Remark 1, Ã0 = U −1 A0 U hence it suffices to search a similarity
transformation U(x). Since Ox is a Principal Ideal Domain (the ideals of Ox
are of the form xk O), it is well known that one can construct unimodular
transformations Q(x), R(x) lying in GLn (Ox ) such that the matrix QA0 R
has the Smith Normal Form
Q A0 R = diag(xβ1 , . . . , xβr , 0, . . . , 0))
where detR(0) 6= 0, detQ(0) 6= 0 , and β1 , . . . , βr in Z with 0 ≤ β1 ≤ β2 ≤
· · · ≤ βr .
It follows that we can compute a unimodular matrix R(x) in GLn (Ox ) so that
its n − r last columns form a basis of ker(A0 ). Hence, we set U(x) = R(x).
Remark 5 In practice, U(x) can be obtained by performing gaussian elimination on the columns of A0 (x) taking as pivots the elements of minimum
order (valuation) in x as already suggested in [Barkatou et al., 2006].
Hence, we can suppose without loss of generality that A0 (x) consists of r independent columnsand (n−r) zero columns. We partition A1 (x) conformally
A11
A12
1
1
, and consider
with A0 , i.e. A1 =
21
A1 A22
1

 11
A0
A12
1
.
(8)
Gλ (A) =
A21
A22
0
1 + λIn−r
We illustrate our progress with the following simple example.
Example 2 Given ǫ dY
= A(x, ǫ)Y where
dx


ǫ −x3 ǫ (1 + x)ǫ
0
 x2
xǫ
0
−2xǫ
.
A = ǫ−2 
−x
0
0
2ǫ 
0
2
0
ǫ2
Clearly, A0 (x) is nilpotent of rank 2

0
 x2
Gλ (A) = 
−x
0

and

0 x+1 0
0
0
−2x
.
0
λ
2 
2
0
λ
10

(9)

This consideration of Gλ (A) gives an ǫ-reduction criteria equivalent to θ(λ)
as demonstrated in the following lemma.
Lemma 2 Det(Gλ (A) ≡ 0 vanishes identically in λ if and only if θ(λ) does.
Proof. Let D(ǫ) = diag(ǫIr , In−r ) then we can write ǫ−1 A = ND −1 where
N := N(x, ǫ) ∈ K n×n has no poles in ǫ. Set D0 = D(0) and N0 = N(x, 0).
Then we have
det(Gλ (A) = det(N0 + λD0 ) = det(N + λD)|ǫ=0
A
= (det( + λIn ) det(D))|ǫ=0
ǫ
A0
+ A1 + λIn ) ǫr )|ǫ=0 = θ(λ).
= (det(
ǫ

Proposition 2 If mǫ (A) > 1 and det(Gλ (A)) ≡ 0 then there exists a unimodular matrix Q(x) in GLn (Ox ) with det Q(x) = ±1 such that the matrix
Gλ (Q[A]) has the form
 11

A0
U1
U2
W2  ,
Gλ (Q[A]) =  V1 W1 + λIn−r−ρ
(10)
M1
M2
W3 + λIρ
where 0 ≤ ρ ≤ n − r, W1 , W3 are square matrices of order (n − r − ρ) and
ρ respectively , and

 11


A0 U1
U1 ),
(11)
) = rank( A11
rank(
0
M1 M2


U1 ) < r.
rank( A11
(12)
0

Our procedure is that of [Barkatou, 1995] except for the properties of M1 ,
M2 , and W3 . The nullity of M1 and M2 in the former is replaced by the
weaker condition (11) here. Otherwise, the unimodularity of Q(x) cannot
be guaranteed. Moreover, this refinement avoids unnecessary computations.
The bridge between both descriptions is established in Remark 6 before proceeding to the proof of the Proposition.
Remark 6 Suppose that Gλ (Q[A]) has the form (10) and there exists a
transformation T (x) ∈ GLn (Kx ) such that Gλ (T [Q[A]]) has the form

 11
A0
U1
U2
W2  ,
Gλ (T [Q[A]]) =  V1 W1 + λIn−r−ρ
O
O
W̃3 + λIρ

11

where 0 ≤ ρ ≤ n − r, W1 , W̃3 are square matrices of order (n − r − ρ) and
ρ respectively, and W̃3 is upper triangular with zero diagonal. Then,
 11

A0
U1
ρ
det Gλ (T [Q[A]]) = λ det
.
V1 W1 + λIn−r−ρ
If det(Gλ (A) ≡ 0) then we have as well det(Gλ (Q[A]) = det(Gλ (T [Q[A]]) ≡
0. Hence,
 11

A0
U1
det
= 0.
(13)
V1 W1 + λIn−r−ρ

 11
A0 U1
(ρ)
(ρ)
.
We shall denote by G0 the matrix G0 =
V1 W1
Proof. (Proposition 2) Since det(Gλ (A)) ≡ 0 then in particular, det(Gλ=0 (A)) =
0. The matrix Gλ=0 (A) is thus singular. Let E (respectively F ) be the vector
space spanned by the first r (resp. last n − r) rows of Gλ=0 (A). We have
dim(E + F ) = rank(Gλ=0 (A)) < n.
If dim(E) < r then set ρ = 0. Otherwise, since
dim(E + F ) = dim(E) + dim(F ) − dim(E ∩ F ) < n,
it follows that either dim(F ) < n − r or dim(E ∩ F ) > 0. In both cases,
there exists at least a row vector W1 (x) = (w1 , . . . , wn ) with entries in C((x))
in the left null space of Gλ=0 (A), such that wi 6= 0 for some r + 1 ≤ i ≤ n.
We can assume without loss of generality that W1 (x) has its entries in Ox .
Indeed, this assumption can be guaranteed by a construction as in Lemma 1.
Let V1 (x) be the nonzero row vector whose entries vi are such that vi = wi
for r + 1 ≤ i ≤ n and zero elsewhere. It can be attained that vn = 1 upon
replacing A by P [A] where P is a permutation matrix which exchanges the
rows of A so that the pivot is taken to have the minimum order (valuation) in
x. Let R1 (x) be the matrix whose nth row is formed by V1 (x), the other rows
being those of the identity matrix In . Note that R1 (x) is a lower triangular
matrix in GLn (Ox ) with det R1 = 1. Now put A(1) = R1−1 [A]. Thus, Gλ (A(1) )
has the form (10) with (11) and ρ = 1.
(1)
By (13), the matrix G0 (x) is singular so if the condition (12) does not
occur then one can find, by same argument as above, a vector V2 (x) in Ox
(1)
whose (n − 1)th component is 1 such that V2 G0 = 0. If R2 (x) denotes the
matrix whose n − 1th row is formed by (V2 (x), 0), the others being those of
the identity matrix In , and A(2) = R2−1 [A(1) ] then the matrix Gλ (A(2) ) has
the form (10) with (11) and ρ = 2.
A finite number of applications of this process yields an equivalent matrix
(ρ)
A with (11) for which either (12) occurs or ρ = n−r. But in the latter case
one has, again by Remark 6 and (13), that det A11
0 = 0, and so (12) occurs.
The matrix Q(x) we are seeking is obtained as a product of permutation
12

matrices or lower triangular matrices whose determinant is 1, and so its
determinant is ±1.
We consider again Example 2.
Example* 1 A simple calculation shows that det(Gλ (A)) ≡ 0 hence A is
ǫ-reducible. From (9), for λ = 0, we have the singular matrix




1 0 0 0
0 0 x+1 0


 x2 0
0
−2x
 . Let Q = 0 1 0 0 ,
Gλ=0 (A) = 


−x 0
0 0 0 1
0
2
0 0 1 0
0 2
0
0


ǫ −x3 ǫ
0
(1 + x)ǫ
 x2
xǫ −2xǫ
0 
.
hence Q[A] = ǫ−2 
2
 0
2
ǫ
0 
−x
0
2ǫ
0
Moreover, Gλ (Q[A])
0 0
 x2 0
Gλ (Q[A]) = 
 0 2
−x 0

has the form
 (10) with ρ = 1 and r = 2. In fact,
0 (1 + x)
0
0 
.
λ
0 
2
λ

Proposition 3 If mǫ (A) > 1 and Gλ (A) has the form (10) with conditions
(11) and (12) satisfied, then system (1) (resp. A) is ǫ-reducible and reduction can be carried out with the so-called shearing transformations S =
diag(ǫIr , In−r−ρ, ǫIρ ) whenever ρ 6= 0 and S = diag(ǫIr , In−r ) otherwise.
Proof. We partition A(x, ǫ) as follows

 11
A
A12 A13
A = A21 A22 A23 
A31 A32 A33
where A11 , A22 , A33 are of dimensions r, n − r − ρ, and ρ respectively. It is
easy to verify then that

 11 −1 12
A
ǫ A
A13
A22
ǫA23  .
Ã(x, ǫ) = S −1 AS − S −1 δS = S −1 AS = ǫA21
31
−1 32
A
ǫ A
A33
Hence, the new leading coefficient matrix is
 11

A0 U1 O
O O
Ã0 =  O
M1 M2 O
and rank(Ã0 ) = rank(A11
0

U1 ) < r.
13

Example* 2 Let S = diag(ǫ, ǫ, 1, ǫ) then


ǫ −x3 ǫ
0 (1 + x)ǫ
 x2
xǫ −2x
0 
.
S[Q[A]] = ǫ−2 
2
 0
2ǫ
ǫ
0 
−x
0
2
0

It is clear that the leading term has rank 1 < r.

Now, we are ready to give the proof of Theorem 1.
Proof. (Theorem 1) We first prove that the condition is necessary. In
fact, suppose that such a T exists then we have det(λI + Aǫ ) = det(λI + Ãǫ ).
It then suffices to see that
A
Ã
ǫr det(λI + )|ǫ=0 = ǫr det(λI + )|ǫ=0 = ǫr−m P (x, ǫ)|ǫ=0
ǫ
ǫ
n×n
where m ≤ r(Ã0 ) < r and P (x, ǫ) ∈ K
has no poles in ǫ. Moreover, it is
evident that θ(λ) = 0 since
∞
X
A
A0
r
r
ǫ det(λI + )|ǫ=0 = ǫ det(λI +
+ A1 +
Ak ǫk−1 )|ǫ=0 .
ǫ
ǫ
k=2
As for the sufficiency of the theorem, by Lemma 1, we can assume that A0 (x)
is in the form (7). Then, Gλ (A) is constructed as in (8). Since θ(λ) = 0
and mǫ (A) > 1, it follows from Lemma 2 that det(Gλ (A)) = 0 and A is
ǫ-reducible. Hence, the matrix S[Q[A]] where S, Q are as in Propositions 2
and 3 respectively, has the desired property.
Remark 7 The ǫ-reducibility of A implies that the rank of the leading coefficient matrix can be reduced without necessarily reducing the ǫ-rank of the
system. If the ǫ-reduction criteria holds for a sufficient number of equivalent
systems then a repetitive application of such a transformation results in an
equivalent system whose leading coefficient matrix has a null rank, hence h is
reduced by one (e.g. Example 3). At this point, the discussion restarts from
the nature of the eigenvalues of the leading constant matrix.
= A(x, ǫ)Y where
Example 3 Given ǫ dY
dx

 3
2xǫ 3x2 ǫ4 2xǫ2 (2x + 1)ǫ5

 0
ǫ4
0
0
.
A = ǫ−3 

 0
0
ǫ2
0
−2x
0
0
0

The gauge transformaiton Y = T Z computed by our algorithm results in an
equivalent ǫ-irreducible system whose ǫ-rank is diminished by two as illustrated below:




1
0
0
0
0 ǫ3 0 0

 0 0 ǫ 0
2xǫ 3x2 2x + 1
−1 2x
.

and
T
[A]
=
ǫ
T =
0
ǫ3 0 0 0
0
ǫ2
0 
0 −2xǫ 0
0
0 0 0 1
14

Algorithm 1 Moser-based ǫ-Rank Reduction of System (1)
Input: A(x, ǫ) of system (1)
Output: T (x, ǫ) ∈ GLn (K) and an equivalent system given by T [A] which
is ǫ-irreducible in the sense of Moser.
T ← In ;
h ← ǫ-rank of A;
U(x) ← unimodular transformation computed in Lemma 1 and Remark 5
such that U −1 A0 U has form (7);
T ← T U; A ← U −1 AU − ǫU −1 dU
;
dx
while Det(Gλ (A)) = 0 and h > 0 do
Q(x), ρ ← Proposition 2 ;
S(ǫ) ← Proposition 3;
P ← QS;
T ← TP;
;
A ← P −1 AP − ǫP −1 dP
dx
h ← ǫ-rank of A;
U(x) ← Lemma 1 and Remark 5;
T ← T U;
A ← U −1 AU − ǫU −1 dU
;
dx
end while.
return (T, A).

15

4.3

Example of Comparison with Levelt’s Approach

While Moser defined two rational numbers in [Moser, 1960], Levelt investigated the existence of stationary sequences of free lattices in [Levelt, 1991].
Since Moser-based and Levelt’s algorithms serve the same utility, i.e. rank
reduction of system (2), it is natural to question their comparison. The cost
analysis in the univariate case of the Moser-based algorithm described in
[Barkatou, 1995] and that of Levelt’s [Levelt, 1991] was given in [page 108,
[LeRoux, 2006]]. Both algorithms turn out to have an identical cost which
suggests a further experimental study so that they can be well-compared.
The latter was generalized in [Barkatou et al., 2006] to a certain class of
Pfaffian systems over bivariate fields. Consequently, following arguments
similar to those of [Barkatou et al., 2006] whenever the leading coefficient
matrix is nilpotent and to Section 3 of this article whenever it’s not, Levelt’s
algorithm seems adaptable to system (1) as well. Such an attempt would
give rise to Algorithm 2.
Algorithm 2 Generalization of Levelt’s Rank Reduction to System (1)
Input: A(x, ǫ) of system (1);
Output: T (x, ǫ) ∈ GLn (K) computed via Levelt’s approach and an
equivalent system given by T [A] which is ǫ-irreducible in the sense of Moser.
T ← In
i ← 0;
h ← ǫ-rank of A;
while i < n − 1 and h > 0 do r ← rank(A0 );
U(x) ← unimodular transformation of Lemma 1 and Remark 5 such
that U −1 A0 U has the form (7);
S(x) ← diag(ǫIr , In−r );
P ← US;
T ← TP;
A ← P −1 AP − ǫP −1 ∂P
;
∂x
h̃ ← ǫ-rank of A;
if h̃ < h then
i ← 0;
else
i ← i + 1;
end if
h ← h̃;
end while.
return (T, A).
It is clear that Algorithm 2 coincides with Algorithm 1 for ρ = 0. For
ρ > 0, which frequently occurs for matrices of dimension greater than 3, we
ran simple examples of systems (1). Despite the identical cost of both algo16

rithms, these examples exhibited that Algorithm 2 complicates dramatically
the coefficient matrices of the system under reduction. One factor in this
complication stems from the weak termination criterion of this algorithm.
However, even upon adjoining Moser’s termination criterion (given by θ(λ))
to this algorithm, as suggested in [Secion 5 of [Barkatou et al., 2006]], the
result remains less satisfying than that of Algorithm 1. We exhibit here a
selected example over O = C[[x, ǫ]]. Additional examples along with our
implementation are available online at [Maddah, 2014].
= ǫ−h x−p A(x, ǫ)Y where
Example 4 Let ǫ dY
dx


0
0 0 ǫx(ǫx + 3) ǫ3 (x + 9) ǫx2
 x
0 0
0
9ǫ2 x2
0 

 2
6
 x
−1 0 ǫx(ǫx + 1)
0
0 
.

A= 3
2

x
+
x
−x
0
ǫx
0
0


 0
0 x
−3ǫ
0
0 
x3 − 1 5 0
0
0
0
It is easily verified that A is Moser-reducible in ǫ. Furthermore, we have
mǫ (A) = h+ 36 and µǫ (A) = h+ 62 . Hence, it suffices to run only one reduction
step, for which the rank of the leading coefficient matrix is dropped by one.
We give in the following the transformation T (x, ǫ) and T [A] as computed
by our Moser-based algorithm, the generalization of Levelt’s algorithm with
Moser’s reducibility criterion, and Levelt’s algorithm respectively. For the
latter, T and T [A] are to be found at [Maddah, 2014] due to the lack of space
here. However we illustrate the dramatic growth of their coefficients by listing
one entry of each.


ǫ 0 0 0 0 0
0 ǫ 0 0 0 0


0 0 ǫ 0 0 0

• T =
0 0 0 0 0 ǫ  and


0 0 0 1 0 0
0 0 0 0 ǫ 0


0
0 0 ǫ2 (x + 9) x2 ǫ ǫx(ǫx + 3)


x
0 0
9x2 ǫ
0
0


2
6

x
−1 0
0
0 ǫx(ǫx + 1)


T [A] = 
2

0
0
ǫx
0
0
−3ǫ

 3

 x −1
5 0
0
0
0
2
2
x(x + 1) −x 0
0
0
xǫ

3
ǫ 0
0
0 0 −9xǫ2
 0 0 −1/3x2 ǫ2 ǫ2 0
0 



0 0
0
0
ǫ
0
 and
• T =
2

0 ǫ
−1/3ǫx
0
0
0


0 0
0
0 0
1 
0 0
ǫ
0 0
0
17

T [A] = [ãij ] including entries having 2-digit coefficients and degree 8 in
x, e.g.,

ã12 = −x(27ǫ2 − ǫx − 3)



1 2


ã13 = 3 x (−x + 27ǫ)
ã41 = ǫx4 − ǫx + 3



ã52 = ǫ2 x(ǫx6 + 1)



ã = −(1/3)ǫ2 x8
53

• T = [tij ], T [A] = [ãij ] have entries with 4-digit and 10-digit coefficients.
Degrees in x surpass 8, e.g.

2
2 5 3
2
11
10

t23 = 1458ǫ x (x + 3x + 2) /(243x + 1944x
+3627x9 + 1026x8 + 1944x7 − 8820x6 + 432x5


−4860x4 − 72x3 − 652x2 − 72x − 324)

25
24
23

ã13 = (39366ǫx + 1220346ǫx + 14565420ǫx +



83731482ǫx22 − 1003833x23 + 236309724ǫx21





−25095825x22 + 708588ǫ2x19 + 299723976ǫx20





−237858120x21 + 17714700ǫ2x18 + 191306610ǫx19





−1080181170x20 + 143134776ǫ2x17 + 140023890ǫx18




−2462664465x19 + 462261816ǫ2x16 − 236799612ǫx17





−3116776563x18 + 553879620ǫ2x15 + 98626896ǫx16





−4040830962x17 + 384763284ǫ2x14 − 334491552ǫx15





−5062097592x16 + 750333456ǫ2x13 + 149450184ǫx14





−3027609900x15 + 681679152ǫ2x12 − 469371996ǫx13



14
2 11
12


−4012634700x + 460735776ǫ x − 88155972ǫx
−1297170936x13 + 804260016ǫ2x10 − 332366796ǫx11



−786840174x12 + 245153952ǫ2x9 − 90553302ǫx10





−404146368x11 + 383522040ǫ2x8 − 73322280ǫx9





+451532556x10 + 111422304ǫ2x7 − 9084492ǫx8





−104014116x9 + 74305512ǫ2x6 − 14703120ǫx7





+222024672x8 + 17635968ǫ2x5 − 7039224ǫx6




−17230320x7 − 23184ǫ2 x4 − 8030664ǫx5 + 88635312x6





−414720ǫ2 x3 − 3070548ǫx4 − 294176x5





−1819584ǫ2 x2 + 32403312x4 + 419904ǫ2x + 1692576x3





+944784ǫ2 + 3895776x2 + 524880x + 944784)/(243x11





+1944x10 + 3627x9 + 1026x8 + 1944x7 − 8820x6



+432x5 − 4860x4 − 72x3 − 652x2 − 72x − 324)2 ;
18

5

Conclusion and Further Investigations

We proposed a Moser-based algorithm which recovers a singularly-perturbed
linear differential system from its turning points and reduces its ǫ-rank to
its minimal integer value. A complementary step to attain the full formal
reduction would be to find the ramification in the parameter which renders
the general case to the case discussed here in a recursive process. One approach is that based on analysis by a Newton polygon and applied to system
(2) in [Barkatou, 1997] and the scalar case of system (1) in [Macutan, 1999].
The sufficient number of coefficient matrices in computations is still to be
investigated.
In the usual treatment of turning points, a restraining index χ is defined
and updated at every reduction step to observe the growth of the order of
poles in the Ak (x)’s (see, e.g. [Wasow, 1985, 1979] and references therein).
This restraining index plays a role in the asymptotic interpretation of the
formal solutions. As demonstrated in Section 3, the transformations we apply
to recover from turning points are polynomial (shearings). Hence, the growth
of the poles order is bounded and can be expected apriori. The insight this
gives into the restraining index is to be investigated.
Examples comparing this algorithm with a generalization of Levelt’s favors the former. However, it suggests that Levelt’s algorithm be generalized
to system (1) and that a bit complexity study comparing both algorithms
be held alongside. Furthermore, it motivates generalization of the Moserbased algorithms over differential bivariate fields, e.g. Pfaffian systems in
two variables [Barkatou et al., 2014].
An additional field of investigation is the two-parameter algebraic eigenvalue problem as a generalization of the one parameter case investigated via
Moser-based approach in [Jeannerod et al., 1999]. In fact, the main role in
the reduction process is reserved to the similarity term of T [A]. Hence, the
discussion of such problems is not expected to deviate from the discussion
presented here in the differential case.

References
M. Barkatou, S.S. Maddah, and H. Abbas. Formal Reduction of a class of
Pfaffian Systems in Two Variables. Submitted to ISSAC’14. A preliminary
version is available at arXiv.
W. Balser. Formal Power Series and Linear Systems of Meromorphic Ordinary Differential Equations. Springer-Verlag, New York, 2000.
M. Barkatou. An algorithm to compute the exponential part of a formal
fundamental matrix solution of a linear differential system. Journal of App.
Alg. in Eng. Comm. and Comp., 8(1):1-23, 1997.
19

M. Barkatou. A Rational Version of Moser’s Algorithm. In Proceedings of the
International Symposium on Symbolic and Algebraic Computation, pages
297-302. ACM Press, July 1995.
M. Barkatou. On the reduction of Linear Systems of Difference Equations.
In Proceedings of the International Symposium on Symbolic and Algebraic
Computation, pages 1-6. ACM Press, USA, 1989.
M. Barkatou, G. Broughton, and E. Pflugel. Regular Systems of Linear Functional Equations and Applications. In Proceedings of the International
Symposium on Symbolic and Algebraic Computation, pages15-22. ACM
Press, 2008.
M. Barkatou and E. Pflugel. ISOLDE, Integration of Systems of Ordinary
Linear Differential Equations. Available at: http://isolde.sourceforge.net/
M. Barkatou and E. Pflugel. On the Moser-and super-reduction algorithms
of systems of linear differential equations and their complexity. Journal of
Sym. Comput., 44 (8), 1017-1036, 2009.
M. Barkatou and N. LeRoux. Rank Reduction of a class of Pfaffian Systems in
Two Variables. In Proceedings of the International Symposium on Symbolic
and Algebraic Computation, pages 204-211. ACM Press, 2006.
G. Chen. Solutions Formelles de Systèmes d’Equations Différentièlles
Linèaires Ordinaires Homogènes. PhD Thesis. Université Joseph Fourier.
Grenoble 1. 1984.
P.F. Hsieh and Y. Sibuya. Basic theory of Ordinary Differential Equations.
Springer. NewYork, USA, 1999.
M. Iwano and Y. Sibuya. Reduction of the Order of a Linear Ordinary Differential Equation Containing a Small Parameter. 1963
C.P. Jeannerod and E. Pflugel. A Reduction Algorithm for Matrices Depending on a Parameter. In Proceedings of the International Symposum
on Symbolic and Algebraic Computation, Pages 121-128. ACM Press, USA
1999.
T. Kato. Perturbation Theory for Linear Operators. Springer. Berlin. 1980.
N. LeRoux. Solutions formelles d’équations aux dérivées partielles. Ph.D.
Thesis. University of Limoges. 2006.
A.H.M. Levelt. Stabilizing Differential Operators: a method for Computing
Invariants at Irregular Singularities. Differential Equations and Computer
Algebra, M. Singer (ed.), pages 181-228, 1991.
C. C. Lin, The theory of hydrodynamic stability. Cambridge Univ. Press.
Cambridge, 1966.
20

Y.O. Macutan. Formal Solutions of Scalar Singularly-Perturbed Linear Differential Equations. In Proceedings of the International Symposum on Symbolic and Algebraic Computation, Pages 113-120. ACM Press, USA 1999.
Maddah S.S. http : //www.unilim.f r/pages perso/suzy.maddah/
J.A.M. McHugh. An historical Survey of Ordinary Linear Differential Equations with a Large Parameter and Turning Points. Archive for History of
Exact Sciences, 7(4): pp 277-324,1971.
J. Moser. The Order of a Singularity in Fuchs’ Theory. Mathematische
Zeitschrift, 72:379- 398, 1960.
E. Pflugel. Effective Formal Reduction of Linear Differential Systems. Journal
of App. Alg. in Eng. Comm. and Comp., 10, 153-187, 2000.
H. L. Turrittin. Asymptotic Expansions of Solutions of Systems of Ordinary
Differential Equations. Contributions to the Theory of Nonlinear Oscillations II. Ann. of Math. Studies. No 29: 81-116. 1952.
W. Wasow. Linear Turning Point Theory. Springer-Verlag. 1985.
W. Wasow. Topics in the Theory of Linear Ordinary Differential Equations
Having Singularities with respect to a Parameter. Institut de Recherche
Mathématique Avancée. Université Louis Pasteur. Strasbourg. 1979.

21

