Isolde, a Package for Computing
Invariants of Systems of Ordinary Linear Differential Equations

Eckhard Pfl�gel

LMC-IMAG, Grenoble

Algorithms Seminar

October 6, 1997

[summary by Fr�d�ric Chyzak]

A properly typeset version of this document is available in postscript and in pdf.

If some fonts do not look right on your screen, this might be fixed by configuring your browser (see the documentation here).

The Maple package Isolde is a package for studying and solving systems of linear differential equations. More specifically, it deals with two main kinds of problems: The approach followed is a direct treatment of the system, avoiding any method akin to that of cyclic vectors. Formal invariants of linear first-order differential systems are introduced in the next section, where we also briefly list the operations available in Isolde. Then, we focus in the last two sections on an efficient algorithm due to E. Pfl�gel to search for exponential solutions [4].

The package Isolde is developed by A. Barkatou and E. Pfl�gel and is available at
http://www-lmc.imag.fr/CF/logiciel.html.

1   Formal Invariants, Formal Solutions, Closed Form Solutions

For a subfield K of the field C of complex numbers with algebraic closure K, and a matrix A� Mn(K(x)), consider the linear first-order differential system
Y'=AY.     (1)
One either looks for vector solutions or for matrix solutions, i.e., whose columns are vector solutions. A formal fundamental matrix at x0�K�{�} is a matrix solution of rank n of the form [6]
F(t)=H(t)t
L
 
eQ(t)     where tr=x-x0 for a positive integer r,     (2)
for a matrix of formal power series H� Mn(K)[[t]], a constant matrix L� Mn(K) and a diagonal matrix Q with Laurent polynomial entries in t-1K[t-1]. Note that
F(t)=H(t)exp �
�		 �
� W(t) dt		 �
�		     for W=L t-1+Q' � t-1 Mn(K)[t-1].
The expression exp(� W(t) dt)=tLeQ(t) is called the exponential part of F. A formal invariant is any quantity appearing in or related to F, like Newton polygons, Newton polynomials, exponential parts and formal solutions, i.e., linear combinations of columns of F. Formal solutions are asymptotic expansions of functional solutions near x0. Their general expression is
y(t)=z(t)t
l
 
eq(t)=z(t)exp �
�		 �
� w(t) dt		 �
�		     (3)
with z(t)=z0(t)+z1(t)ln t+...+zs(t)lnst for vectors zi�K[[t]]n, a constant l�K and a Laurent polynomial q� t-1K[t-1]. Here again w=l t-1+q'. Isolde implements algorithms to compute several formal invariants, in particular exponential parts and formal solutions.

Specialized algorithms have been developed to solve for solutions in several elementary classes of closed form expressions. Their principle is that local data contains essential information on potential closed form solutions. In particular, Isolde implements algorithms to solve for:

2   Exponential Solutions

We now use the formal invariants previously introduced to compute exponential solutions of a linear differential system. Throughout this section, we assume that we know the set of all possible exponential parts of exponential solutions, whose determination is the topic of the next section.

An exponential solution is a vector solution y obtained when the vector z in the general expression (3) reduces to a polynomial vector:
y(x)=p(x)exp �
�		 �
� u(x) dx		 �
�		     for u�K(x) and p�K[x]n.
As an example, an exponential solution of the system described by the matrix
A=
�
�
�
�
�
�
�
�
�
�
�
�
�
-7
x2
-x2
x2
1-x
6-2x+x2
x6(1-x)
1
x2
-4+9x
x(1-x)2
8-8x+5x2
x6
1-x
x2
-4+8x
x(1-x)
�
�
�
�
�
�
�
�
�
�
�
�
�
is exp �
�		 �
� u dx		 �
�
�
�
�
�
0
1
1-x
�
�
�
�
, where u=
1
x2
-
4
x
+
5
1-x
.     (4)
Note the more explicit form exp(� u dx)=x-4(1-x)-5exp(-x-1).

Methods to determine exponential solutions of a system first evaluate or give constraints on the exponential part before computing the polynomial part p. This bases on the local analysis of u. One distinguishes between the singular and regular parts Sx0(u) and Rx0(u) of u at a point x0�K�{�}, defined by:
u=S
 
x0
(u)+R
 
x0
(u),     with Sx0(u)� t-1K[t-1] and Rx0(u)�K[[t]].
For the example above, one gets
S0(u)=
1
t2
-
4
t
    and     R0(u)=5+5t+5t2+���

Basic known algorithms allow us:
  1. for each singularity x0 of a given rational function u, to compute the singular parts Sx0(u) and a finite number of terms of the regular parts Rx0(u) in an efficient way;
  2. to reconstruct a rational function u from its singular parts at all its singularities.
Two different methods based on these algorithms allow the calculation of exponential solutions:
  1. Beke's method:
    1. compute candidates for u by combining the singular parts allowed at all singularities;
    2. set Y=Zexp(� u dx) and search for polynomial solutions Z.
  2. Alternative approach based on Pad� approximants:
    1. bound numerator and denominator of u;
    2. determine the singular part Sx0(u) of u at a (single) singularity x0 so as to be able to compute a Pad� approximant for Rx0(u) next.
The methods contrast to one another inasmuch as Beke's method requires splitting fields for its completeness but avoids the use of Gr�bner bases, while the alternative Pad�-based approach does not appeal to splitting fields but generally requires Gr�bner bases. Besides, note the combinatorial exponential complexity of step (a) in Beke's algorithm.

Solving the example (4) by Beke's method, one first obtains the following sets of exponential parts at 0 and 1:
E0= �
�
�
1
x2
-
4
x
,
a
x2
+
20a-44
57x
�
�
�		 ,     E1= �
�
�		 0,
5
1-x
�
�
�		 ,     where a2+7a-2=0.
With this simple example, this yields only two candidates:
u1=
1
x2
-
4
x
+
5
1-x
    and     u2=
a
x2
+
20a-44
57x
5
1-x
.
With the first candidate, the system (1) reduces to Z'=(A-u1)Z. One then verifies that it admits the polynomial solution already mentioned in (4).

For a formal solution y=zexp(� w dt) such that z has valuation 0 at x0, w is called a generalized exponent. Noticing that each Sx0(u) is a generalized exponent, the idea of the second method is to compute Rx0(u) from a formal series solution, more specifically from the series z after setting w=Sx0(u). More explicitly, for an exponential solution y=pexp(� u dt) with logarithmic derivative U, we have Sx0(u)=Sx0(U), so that
R
 
x0
(U)=R
 
x0
(u)+(ln p)'=(ln z)'
is a vector of rational functions, which can be computed from z as a Pad� approximant. By integration, one obtains Rx0(u) as the rational part and p as the logarithmic part, next an exponential solution if there exists any for this u.

Following up the example (4), choosing the singularity x0=0 and the exponential part candidate w=S0(u)=x-2-4x-1, one first computes the following formal solution (at 0):
y=exp �
�		 �
� S0(u) dx		 �
�
�
�
�
�
0
1+5x+15x2+���
1+4x+14x2+���
�
�
�
�
.
Taking the logarithmic derivative and computing Pad� approximations yields
R0(U)=
�
�
�
�
0
5+5x+���
4+4x+���
�
�
�
�
=
�
�
�
�
�
�
�
�
�
�
0
5
1-x
4
1-x
�
�
�
�
�
�
�
�
�
�
= -
�
�
�
�
0
ln(1-x)5
ln(1-x)4
�
�
�
�
',
from which one gets the exponential solution already mentioned in (4).

3   Exponential Parts

To complete the description of the previous algorithms, we finally describe how to compute all possible exponential parts of solutions of a linear first-order ordinary differential system.

As opposed to the case of a (single) scalar higher order ODE, exponential parts cannot be obtained immediately in the case of the system (1). An obvious indirect method is to compute the Newton polygon by transforming it into an nth order scalar equation, which becomes very costly with the increase of n. This is why several other methods have been developed to transform the initial matrix A into a form from which exponential parts can be read off. All such methods use two special transformations:
  1. the change of unknown Y=TZ, where T is a polynomial matrix with non-zero determinant transforms the system (1) into the equivalent system Z'=(T-1AT-T-1T')Z;
  2. the change of exponential part Y=Zexp(� a dt) leads to the new system Z'=(A-a)Z.
Using the above, algorithms have been proposed to put a differential system into:
  1. companion form, as obtained by the method of cyclic vectors;
  2. Turrittin's canonical form [5], however obtained by a not so constructive process;
  3. Moser's irreducible form [3];
  4. Hilali's and Wazner's super irreducible form, a refined version of Moser's form [2].
We now comment on applications of the last two forms. If a system
xq+1Y'=AY,    for a series A=A0+A1x+���     (5)
admits a solution of the form
zexp �
�
�
�		 �
�
a
xq+1
 dx �
�
�
�		 ,     (6)
then the matrix A0-a necessarily has a zero determinant. Elaborating on this fact, Moser [3] proved that when A is an irreducible system, if A0 is not nilpotent and a is a non-zero root of det(A0-l) with multiplicity m, then there exist m solutions of A of the form (6). Based on this, Barkatou used diagonalization by blocks to devise an algorithm to compute exponential parts [1].

Similarly, a necessary condition for a system like (5) to admit a solution of the form
zexp �
�
�
�		 �
�
a
xk+1
 dx �
�
�
�		 ,     for 0� k� q,
takes the form det(N0-aD0)=0 for two matrices N0 and D0 computed from A0,...,Aq-k. Consider the non-zero polynomial
qk(l)=xsdet(x-kA+l)|x=0
obtained for the appropriate exponent s. Hilali and Wazner [2] proved that when A is a super irreducible system, if a is a non-zero root of multiplicity m of the polynomial qk, then there exist precisely m generalized exponents equal to -ax-(q-k) up to higher valuation terms. Using this fact, Pfl�gel [4] obtained a recursive algorithm to compute all exponential parts of ramification 1, i.e., for the case r=1 in (2). The case of higher ramifications r is work in progress.

References

[1]
Barkatou (M. A.). -- An algorithm to compute the exponential part of a formal fundamental matrix solution of a linear differential system. Applicable Algebra in Engineering, Communication and Computing, vol. 8, n°1, 1997, pp. 1--23.

[2]
Hilali (A.) and Wazner (A.). -- Formes super-irr�ductibles des syst�mes diff�rentiels lin�aires. Numerische Mathematik, vol. 50, n°4, 1987, pp. 429--449.

[3]
Moser (J�rgen). -- The order of a singularity in Fuchs' theory. Mathematische Zeitschrift, vol. 72, 1959/1960, pp. 379--398.

[4]
Pfl�gel (Eckhard). -- An algorithm for computing exponential solutions of first order linear differential equations. In K�chlin (Wolfgang W.) (editor), ISSAC'97 (July 21--23, 1997. Maui, Hawaii, USA). pp. 164--171. -- ACM Press, 1997. Proceedings of ISSAC'97. ISBN 0-89791-875-4.

[5]
Turrittin (H. L.). -- Reduction of ordinary differential equations to the Birkhoff canonical form. Transactions of the American Mathematical Society, vol. 107, 1963, pp. 485--507.

[6]
Wasow (Wolfgang). -- Asymptotic expansions for ordinary differential equations. -- Dover Publications Inc., New York, 1987, x+374p. Reprint of the John Wiley 1976 edition.
This document was translated from LATEX by HEVEA.