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

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:

local problems: compute formal invariants at a point; compute formal solutions at a point;
global problems: compute closed form solutions in a certain class, like that of polynomial, rational, or exponential functions.

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Î M_n(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 x₀ÎKÈ{¥} is a matrix solution of rank n of the form [6]

F(t)=H(t)t

e^Q(t) where t^r=x-x₀ for a positive integer r, (2)

for a matrix of formal power series HÎ M_n(K)[[t]], a constant matrix LÎ M_n(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 M_n(K)[t^-1].

The expression exp(ò W(t) dt)=t^Le^Q(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 x₀. Their general expression is

y(t)=z(t)t

e^q(t)=z(t)exp

æ
è

ó
õ w(t) dt

ö
ø

(3)

with z(t)=z₀(t)+z₁(t)ln t+...+z_s(t)ln^st for vectors z_iÎK[[t]]ⁿ, 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:

polynomial solutions, that are related to formal solutions at ¥ with trivial exponential parts and with no logarithmic component;
rational solutions, whose denominators are bounded by computing the indicial equation (i.e., the Newton polynomial of slope 0) at all finite singularities;
exponential solutions, for which candidates can be computed from the exponential parts at all singularities, as detailed in the next two sections.

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]ⁿ.

As an example, an exponential solution of the system described by the matrix

é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë

-7

x²

-x²

x²

1-x

6-2x+x²

x⁶(1-x)

x²

-4+9x

x(1-x)²

8-8x+5x²

x⁶

1-x

x²

-4+8x

x(1-x)

ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û

is exp

æ
è

ó
õ u dx

ö
ø

é
ê
ê
ë

1-x

ù
ú
ú
û

, where u=

x²

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 S_x_₀(u) and R_x_₀(u) of u at a point x₀ÎKÈ{¥}, defined by:

u=S

x₀

(u)+R

x₀

(u), with S_x_₀(u)Î t^-1K[t^-1] and R_x_₀(u)ÎK[[t]].

For the example above, one gets

S₀(u)=

t²

and R₀(u)=5+5t+5t²+···

Basic known algorithms allow us:

for each singularity x₀ of a given rational function u, to compute the singular parts S_x_₀(u) and a finite number of terms of the regular parts R_x_₀(u) in an efficient way;
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:

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.
Alternative approach based on Padé approximants:
1. bound numerator and denominator of u;
2. determine the singular part S_x_₀(u) of u at a (single) singularity x₀ so as to be able to compute a Padé approximant for R_x_₀(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:

E₀=

ì
í
î

x²

20a-44

57x

ü
ý
þ

, E₁=

ì
í
î

1-x

ü
ý
þ

, where a²+7a-2=0.

With this simple example, this yields only two candidates:

u₁=

x²

1-x

and u₂=

x²

20a-44

57x

1-x

With the first candidate, the system (1) reduces to Z'=(A-u₁)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 x₀, w is called a generalized exponent. Noticing that each S_x_₀(u) is a generalized exponent, the idea of the second method is to compute R_x_₀(u) from a formal series solution, more specifically from the series z after setting w=S_x_₀(u). More explicitly, for an exponential solution y=pexp(ò u dt) with logarithmic derivative U, we have S_x_₀(u)=S_x_₀(U), so that

x₀

(U)=R

x₀

(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 R_x_₀(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 x₀=0 and the exponential part candidate w=S₀(u)=x^-2-4x^-1, one first computes the following formal solution (at 0):

y=exp

æ
è

ó
õ S₀(u) dx

ö
ø

é
ê
ê
ë

1+5x+15x²+···

1+4x+14x²+···

ù
ú
ú
û

Taking the logarithmic derivative and computing Padé approximations yields

R₀(U)=

é
ê
ê
ë

5+5x+···

4+4x+···

ù
ú
ú
û

é
ê
ê
ê
ê
ê
ê
ê
ê
ë

1-x

ù
ú
ú
ú
ú
ú
ú
ú
ú
û

= -

é
ê
ê
ë

ln(1-x)⁵

ln(1-x)⁴

ù
ú
ú
û

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:

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;
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:

companion form, as obtained by the method of cyclic vectors;
Turrittin's canonical form [5], however obtained by a not so constructive process;
Moser's irreducible form [3];
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

x^q+1Y'=AY, for a series A=A₀+A₁x+··· (5)

admits a solution of the form

zexp

æ
ç
ç
è

ó
õ

x^q+1

ö
÷
÷
ø

, (6)

then the matrix A₀-a necessarily has a zero determinant. Elaborating on this fact, Moser [3] proved that when A is an irreducible system, if A₀ is not nilpotent and a is a non-zero root of det(A₀-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

æ
ç
ç
è

ó
õ

x^k+1

ö
÷
÷
ø

, for 0£ k£ q,

takes the form det(N₀-aD₀)=0 for two matrices N₀ and D₀ computed from A₀,...,A_q-k. Consider the non-zero polynomial

q_k(l)=x^sdet(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 q_k, 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 L^AT_EX by H^EV^EA.