Linear Ordinary Differential Equations: breaking through the order 2 barrier

Manuel Bronstein

Institut für Wissenschaftliches Rechnen

ETH - Zentrum CH-8092 Zürich, Switzerland

bronstein@inf.ethz.ch

Abstract

A major subproblem for algorithms that either factor ordinary linear differential equations or compute their closed form solutions is to find their solutions y which satisfy y′/y ∈ \overline{K}(x) where K is the constant field for the coefficients of the equation. While a decision procedure for this subproblem was known in the 19^{th} century, it requires factoring polynomials over \overline{K} and has not been implemented in full generality. We present here an efficient algorithm for this subproblem, which has been implemented in the AXIOM computer algebra system for equations of arbitrary order over arbitrary fields of characteristic 0. This algorithm never needs to compute with the individual complex singularities of the equation, and algebraic numbers are added only when they appear in the potential solutions. Implementation of the complete Singer algorithm for n = 2, 3 based on this building block is in progress.

1 Introduction

Consider the linear ordinary differential equation

A_{n}(x)y^{(n)} + … + A_{1}(x)y′ + A_{0}(x)y = 0  (1)

where the A_{i}’s are polynomials with coefficients in some constant field K. The problem of finding whether (1) has a “closed-form” solution has been of great interest to mathematicians since the second half of the 19^{th} century. While the problem of finding all the Liouvillian solutions of (1) is decidable in theory for any order ([8]), the published decision procedure is not considered a practical algorithm and has not been implemented. Currently, even the most advanced computer algebra systems are limited to the second order case (n = 2), for which they use Kovacic’s algorithm ([6]). Even for that case, implementations are usually incomplete and authors have reported difficulties in dealing with arbitrary complex singularities ([11]). These difficulties usually stem from the need of those algorithms to compute with the singularities of (1) in \overline{K}.

We describe here an efficient algorithm, fully implemented in the AXIOM computer algebra system, that returns a basis over \overline{K} for the space of the solutions y of (1) such that y′/y ∈ \overline{K}(x) (we call them the exponential solutions of (1) over \overline{K} in the rest of this paper). The algorithm is “rational” in the sense that algebraic elements over K are introduced only if they appear in the potential solutions, and not in the singularities of the equation as was previously required. Thus, we do not need to factor A_{n}(x), even over K.

While the exponential solutions form only a subspace of the Liouvillian solutions, the main subalgorithm of the algorithms that find the Liouvillian solutions ([6, 8]) is to find the exponential solutions of higher order equations. So, as applications of the algorithm described here we get:

• a “rational” version of Kovacic’s algorithm for n = 2,

• with the bounds of [10], an efficient version of Singer’s algorithm for n = 3,

• a version of the linear differential operator factoring algorithm of [7] which factors over \overline{K}(x) (the existing implementation of that algorithm factors over Q(x) only).

All of those applications are currently being implemented ^{1}. In addition, the search for the exponential solutions alone is being used as an incomplete n^{th}-order solver in the AXIOM computer algebra system, and is successful on a large number of higher order equations from [5] and other sources.

1The last one in conjunction with F.Schwarz at the GMD, Bonn

2 Rational solutions

A solution y of (1) is called rational if y ∈ \overline{K}(x). The rational solutions form a finite dimensional vector space \overline{V} over \overline{K}. Let V be the vector space generated over K by the solutions of (1) which are in K(x). By Theorem 9.1 of [4], dim_{K}(V) = dim\overline{K}(\overline{V}) and a basis for V over K is also a basis for \overline{V} over \overline{K}, so it is sufficient to look for the rational solutions in K(x). A fast rational algorithm that finds a basis for V over K is described in [1, 2, 3] and has been implemented in the computer algebra systems AXIOM, Maple and REDUCE, so we do not repeat its description here. Its main innovation is the use of balanced factorisation which we outline here since it will be used in the next part of the algorithm.

Recall that for P, Q ∈ K[x] \ {0}, deg(P) > 0, the order of Q at P, denoted ν_{P} (Q) is the unique integer m such that P^{m} | Q and P^{m+1} |̸ Q.

Definition 1 Let A, B ∈ K[x]. We say that A is balanced w.r.t. B if ν_{P} (B) = ν_{Q}(B) for any two irreducible factors P and Q of A in K[x]. Let S ⊆ K[x]. We say that A = A_{1}^{e_{1}} ⋯ A_{m}^{e_{m}} is a balanced factorisation of A w.r.t S if each A_{i} is squarefree and balanced w.r.t B for any B ∈ S, and (A_{i}, A_{j}) = (1) for i ≠ j.

In other words, a polynomial A is balanced w.r.t. to B if all the irreducible factors of A appear with the same power in B. A balanced factorisation of A w.r.t. a set of polynomials S is a factorisation of A into squarefree coprime factors, each of which is balanced w.r.t. each element of S.

An algorithm for computing balanced factorisations appeared in [1] and slightly modified in [3]. The interesting point computationally is that it only requires gcd computations in K[x], and is at most as expensive as computing a square-free-gcd-free basis for {A} ∪ S.

3 The associated Riccati equation

In order to look for the other exponential solutions of (1) over \overline{K}, we consider the associated Ricatti equation defined as follow: let u be a new differential indeterminate and define the sequence of differential polynomials (P_{i})_{i≥0} by

\{ \begin{matrix}P_{0} & = & 1 \\ P_{i} & = & P_{i−1}′ + uP_{i−1} for i > 0. \end{matrix}  (2)

i.e. P_{0} = 1, P_{1} = u, P_{2} = u′ + u^{2}, P_{3} = u′′ + 3uu′ + u^{3}, etc...

Definition 2 The associated Riccati equation of (1) is the non-linear ordinary differential equation:

A_{n}(x)P_{n} + … + A_{1}(x)P_{1} + A_{0}(x)P_{0} = 0. (3)

Example:The associated Ricatti equation of y′′ = ry, r ∈ K(x) (the equation solved by Kovacic’s algorithm) is u′ + u^{2} = r, which is a classical Ricatti equation.

The key property of the associated Ricatti equation, which is proven by substituting e^{∫ u} for y in equation (1) is: if y ≠ 0 is any solution of (1), then u = y′/y is a solution of (3); conversely, if u is any solution of (3),

then y = e^{∫ u} is a solution of (1).

Since for an exponential solution y of (1), u = y′/y ∈ \overline{K}(x), finding the exponential solutions of (1) is equivalent to finding the rational solutions of (3).

4 Rational solutions of the Ricatti equation

We now look for the solutions u of (3) which are in \overline{K}(x). In opposition to the linear case, it can happen that there are solutions in \overline{K}(x) and none in K(x). For example the associated Riccati equation of y′′ + y = 0 is u′ + u^{2} = −1 which has the 2 solutions ±\sqrt{−1} but no solution in Q(x). So it may be necessary to extend K at some point in the algorithm.

In order to bound the order of the singularities of a solution of (3), we need the following construction which is somewhat similar to the Newton polygon of equation (1) at a singularity:

Definition 3 Let S = {Q_{1}, …, Q_{m}} ⊂ K[x]. We define the Riccati exponents of S at ∞ to be:

Γ_{∞}(S) = \{\frac{deg(Q_{i}) − deg(Q_{j})}{i − j} for 1 ≤ i ≠ j ≤ n

such that Q_{i} ≠ 0 ≠ Q_{j}} ∩ N \ {0}

where N is the set of the nonnegative integers, and the Riccati bound of S at ∞ to be:

δ_{∞}(S) = \{ \begin{matrix}0 & if Γ_{∞}(S) = ∅ \\ max(Γ_{∞}(S)) & otherwise. \end{matrix}

Similarly, for P ∈ K[x]\{0}, deg(P) > 0, we define the Riccati exponents of S at P to be:

Γ_{P} (S) = \{\frac{ν_{P} (Q_{i}) − ν_{P} (Q_{j})}{i − j} for 1 ≤ i ≠ j ≤ n

such that Q_{i} ≠ 0 ≠ Q_{j}} ∩ N \ {0}

and the Riccati bound of S at P to be:

δ_{P} (S) = \{ \begin{matrix}1 & if Γ_{P} (S) = ∅ \\ max(Γ_{P} (S)) & otherwise. \end{matrix}  (4)

Our algorithm is based on the algorithm of [9], but uses the following Theorem to avoid factoring A_{n}. This Theorem is essentially the base case (′ = d/dx) of Theorem 8.4 of [4]:

Theorem 1 Let A_{0}, …, A_{n} ∈ K[x] with A_{n} ≠ 0 and monic. Let C_{1}^{e_{1}} … C_{q}^{e_{q}} be a balanced factorisation of A_{n} w.r.t. S = \{A_{0}, …, A_{n}\}. Then, any solution u ∈ \overline{K}(x) of equation (3) can be written as

u = \frac{B_{1}}{C_{1}^{δ_{C_{1}} (S)}} + … + \frac{B_{q}}{C_{q}^{δ_{Cq} (S)}} + \frac{Q′}{Q} + P  (5)

where δ_{C_{i}}(S) is given by (4), Q, P, B_{1}, …, B_{q} ∈ \overline{K}[x], (Q, A_{n}) = (1) and deg(B_{i}) < δ_{C_{i}}(S) deg(C_{i}) for each i.

The algorithm proceeds by computing a finite set of possibilities for each summand of (5). Everytime an additive component α is found, we do the change of variable u = α + \overline{u} in equation (3) and obtain a new Ricatti equation for \overline{u}. We show first how such a change of variable can be done by computing with operators over K(x, α).

4.1 Additive change of variable in the Ricatti equation

Suppose that we have found a component α ∈ \overline{K}(x) of the sum (5). We set u = α + \overline{u} and want a new Ricatti equation for \overline{u}. Let E = K(x, α)[∂_{x}] be the ring of the linear ordinary differential operators over K(x, α) where ∂_{x} is differentiation with respect to x (such a ring is a provided as a type in AXIOM, and is easy to implement in the other systems). The operator corresponding to

equation (1) is L = ∑_{i=0}^{n} A_{i}∂_{x}^{i}. Let y = e^{∫ u} and \overline{y} =

e^{∫ \overline{u}}. Since u = α + \overline{u}, we have y = e^{∫ α}\overline{y}. Also, if u is a solution of (3), then Ly = 0, so L(\overline{y}e^{∫ α}) = 0. But,

from the arithmetic in E, L(\overline{y}e^{∫ α}) = e^{∫ α}\overline{L}\overline{y} where \overline{L} = ∑_{i=0}^{n} A_{i}(∂_{x} + α)^{i}. So the Ricatti equation satisfied by \overline{u} is just the associated Ricatti equation of \overline{L}, and the change of variable procedure is:

change var(A_{0}, …, A_{n}, α) OUTPUT: C_{0}, …, C_{n} such that when u = α + \overline{u} is a solution of (3), \overline{u} is a solution of ∑_{i=0}^{n} C_{i}P_{i} = 0.

• \overline{L} ← ∑_{i=0}^{n} A_{i}(∂_{x} + α)^{i} ∈ K(x, α)[∂_{x}]

• return [coefficient of ∂_{x}^{i} in \overline{L} for i = 0, …, n]

For the case where arithmetic in E is not readily available, there is a closed-form formula for the operator \overline{L} of the first step, namely

(∂_{x} + α)^{i} = ∑_{j=0}^{i} \binom{i}{j}P_{i−j}(α)∂_{x}^{j}

where the P_{k}(α)’s are given by the recursion (2) with u replaced by α.

4.2 The polynomial part

We introduce the notation lc(P) for the leading coefficient of a polynomial P, and we use the notation of Theorem 1 in the rest of this paper. From the proof of Proposition 2.3 of [9] we have deg(P) ≤ δ_{∞}(S), and if deg(P) = m, then H_{m}(lc(P)) = 0 where H_{m} ∈ K[Y] is given by:

d_{m} ← max ; \begin{matrix}0 ≤ i ≤ n \\ A_{i} ≠ 0 \end{matrix} (im + deg(A_{i}))

H_{m} = ∑ ; \begin{matrix}0 ≤ i ≤ n \\ A_{i} ≠ 0 \\ im + deg(A_{i}) = d_{m} \end{matrix} lc(A_{i})Y^{i}.  (6)

So the potential leading coefficients of P are the roots of H_{m}, and the algorithm that returns all the potential polynomial parts in (5) is:

poly part(A_{0}, …, A_{n})

OUTPUT: A set {P_{1}, …, P_{t}} of polynomials such that any solution u ∈ \overline{K}(x) of (3) can be written in the form (5) with P = P_{i} for some i.

• return

poly part bounded(A_{0}, …, A_{n},

δ_{∞}({A_{0}, …, A_{n}}))

poly part bounded(A_{0}, …, A_{n}, m)

• ans ← {0}

• if m < 0 then return ans

• H_{m} ← polynomial given by (6)

• for a ∈ \overline{K} | H_{m}(a) = 0 do α ← ax^{m} {D_{0}, …, D_{n}} ← change var(A_{0}, …, A_{n}, α) {P_{1}, …, P_{t}} ←

poly part bounded(D_{0}, …, D_{n}, m − 1) ans ← ans ∪ {α + P_{1}, …, α + P_{t}}

• return ans

4.3 The P-adic part

Notation: for A, B ∈ K[x], B ≠ 0, we write ϕ_{B}(A) for the remainder of A modulo B. Let C be one of the C_{i}’s of Theorem 1, δ be δ_{C}(S) and rewrite the term R/C^{δ} of the sum (5) as

\frac{R}{C^{δ}} = \frac{R_{δ}}{C^{δ}} + … + \frac{R_{1}}{C}  (7)

where deg(R_{i}) < deg(C). We call (7) the P-adic part of u at C. From the proofs of Theorem 8.4 of [4] and Proposition 2.3 of [9], we have

H_{δ}(R_{δ}) ≡ 0 (mod C)  (8)

where H_{δ} ∈ K[x][Y] is given by

n_{δ} ← max ; \begin{matrix}0 ≤ i ≤ n \\ A_{i} ≠ 0 \end{matrix} (iδ − ν_{C}(A_{i}))

H_{δ} = ∑ ; \begin{matrix}0 ≤ i ≤ n \\ A_{i} ≠ 0 \\ iδ − ν_{C}(A_{i}) = n_{δ} \end{matrix} ϕ_{C}(\frac{A_{i}}{C^{ν_{C}(A_{i})}} )S_{δ,i}(Y). (9)

where

S_{δ,i}(Y) = \{ \begin{matrix}Y^{i} & if δ > 1 \\ ∏_{j=0}^{i−1}(Y − jC′) & if δ = 1. \end{matrix}

So the potential leading coefficients of (7) are the solutions of (8) in \overline{K}[x] of degree at most deg(C)−1. There are two alternatives for solving equation (8):

1. Since deg(R_{δ}) < deg(C) we write R_{δ} = a_{0} +a_{1}x+ … a_{m−1}x^{m−1} where m = deg(C) and the a_{i}’s are indeterminates. Plugging into equation (8), reducing modulo C and equating to 0, we get a zero-dimensional non-linear algebraic system for a_{0}, …, a_{m−1}, and we use the existing non-linear algebraic solvers to find all the solutions in \overline{K}.

2. The algorithm of [9] factors C into irreducibles over K, then for each factor Q_{i}, it factors H_{δ}(γ)(Y) over \overline{K} where Q_{i}(γ) = 0. For each linear factor Y − S(γ) of H_{δ}(γ)(Y), R = S(x) is a potential solution of (8).

Our implementation uses the first method (non-linear algebraic system), since Gröbner basis implementations are available over arbitrary fields and solvers are now reasonably efficient for zero-dimensional systems with a small number of solutions (which are the systems arising when the initial differential equation has a small order). Depending on the actual constant field in each case, it can payoff to factor C further thereby reducing the number of variables, but testing with examples from [5] has not shown the need for factoring C.

With a solver for equation (8) available, the algorithm for computing all the potential P-adic parts at C in (5) parallels exactly the algorithm for the polynomial part (it is in fact the same if one considers the degree function to be the negative of the order at infinity):

padic part(A_{0}, …, A_{n}, C) OUTPUT: A set {f_{1}, …, f_{t}} of fractions such that any solution u ∈ \overline{K}(x) of (3) can be written in the form (5) with B/C^{δ_{C}(S)} = f_{i} for some i.

• return

padic part bounded(A_{0}, …, A_{n}, C,

δ_{C}({A_{0}, …, A_{n}}))

padic part bounded(A_{0}, …, A_{n}, C, δ)

• ans ← {0}

• if δ < 1 then return ans

• H_{δ} ← polynomial given by (9)

• for R ∈ \overline{K}[x] | deg(R) < deg(C)&H_{δ}(R) ≡ 0(C) do α ← R/C^{m} {D_{0}, …, D_{n}} ← change var(A_{0}, …, A_{n}, α) {f_{1}, …, f_{t}} ←

padic part bounded(D_{0}, …, D_{n}, C, δ − 1) ans ← ans ∪ {α + f_{1}, …, α + f_{t}}

• return ans

Once the individual candidates for each P-adic part can computed, finding the entire singular part of (5) (except for the Q′/Q term) is done by computing all the potential P-adic parts at each balanced factor C_{i} and grouping them:

singular part(A_{0}, …, A_{n}) OUTPUT: A set {f_{1}, …, f_{t}} of fractions such that any solution u ∈ \overline{K}(x) of (3) can be written in the form f_{i} + Q′/Q + P with Q and P as in (5).

• C_{1}^{e_{1}} … C_{q}^{e_{q}} ← balanced factorisation of A_{n} w.r.t. {A_{0}, …, A_{n}}

• return singular part1(A_{0}, …, A_{n}, C_{1}, …, C_{q}, q)

singular part1(A_{0}, …, A_{n}, C_{1}, …, C_{m}, m)

• ans ← {0}

• if m < 1 then return ans • {f_{1}, …, f_{t}} ← padic part(A_{0}, …, A_{n}, C_{m})

• for i = 1 … t do {D_{0}, …, D_{n}} ← change var(A_{0}, …, A_{n}, f_{i}) {g_{1}, …, g_{s}} ←

singular part1(D_{0}, …, D_{n},

C_{1}, …, C_{m−1}, m − 1)

ans ← ans ∪ {f_{i} + g_{1}, …, f_{i} + g_{s}}

• return ans

4.4 The logarithmic derivative part

We can now compute a finite set of possibilities for all the quantities appearing in (5) except for the Q′/Q term. For each possibility α we do the change of variable u = α + Q′/Q and obtain a Ricatti equation R = ∑_{i=0}^{n} D_{i}P_{i} such that R(Q′/Q) = 0. But, since

e^{∫ Q′/Q} = Q, we get L(Q) = 0 where L = ∑_{i=0}^{n} D_{i}∂_{x}^{i}. We can find a basis {f_{1}, …, f_{z}} for the rational solutions of Ly = 0 by the algorithm of section 2. If this basis is empty, then equation (3) has no solution with this particular α, otherwise {α + f_{1}, …, α + f_{z}} are all solutions of (3). After grouping all such solutions for the various possible α’s we get all the solutions of (3), hence a basis for the exponential solutions of (1).

solve ricatti(A_{0}, …, A_{n}) OUTPUT: The set of all the solutions of equation (3) in \overline{K}(x).

• ans ← {0}

• {P_{1}, …, P_{t}} ← poly part(A_{0}, …, A_{n})

• for i = 1 … t do {D_{0}, …, D_{n}} ← change var(A_{0}, …, A_{n}, P_{i}) {g_{1}, …, g_{s}} ← singular part(D_{0}, …, D_{n}) for j = 1 … s do

{F_{0}, …, F_{n}} ← change var(D_{0}, …, D_{n}, g_{j}) {h_{0}, …, h_{z}} ←

basis for the rational solutions of ∑_{i=0}^{n} F_{i}y^{(i)} = 0

ans ← ans ∪ {f_{i} + g_{j} + h_{1}′/h_{1}, …, f_{i} + g_{j} + h_{z}′/h_{z}}

• return ans

5 Example

Consider the equation

(x^{2} − x^{4})y′′ − (8x^{3} − 5x + 1)y′ + (3 − 12x^{2})y = 0 (10)

which has no rational solution. Its associated Ricatti equation is

(x^{2} −x^{4})(u′ +u^{2})−(8x^{3} −5x+1)u+3−12x^{2} = 0. (11)

A balanced factorisation of A_{2} = x^{2} − x^{4} w.r.t. S = {A_{0}, A_{1}, A_{2}} is A_{2} = C_{1}C_{2}^{2} = (1 − x^{2})x^{2}, and the order information of the coefficients of (10) is in the following table:

\begin{matrix}i & A_{i} & ν_{x}(A_{i}) & ν_{1−x^{2}}(A_{i}) & deg(A_{i}) \\ 0 & 3 − 12x^{2} & 0 & 0 & 2 \\ 1 & −8x^{3} + 5x − 1 & 0 & 0 & 3 \\ 2 & x^{2} − x^{4} & 2 & 1 & 4 \end{matrix} .

(12)

The polynomial part: From (12), Γ_{∞}(S) = {1, 2} so δ_{∞}(S) = 2. By formula (6), d_{2} = max_{i=0,1,2}(2i + deg(A_{i})) = 8, and H_{2} = lc(A_{2})Y^{2} = −Y^{2}. The only root of H_{2} is 0, so the degree of the polynomial part of a rational solution of (11) is at most 1. Again, by formula (6), d_{1} = max_{i=0,1,2}(i + deg(A_{i})) = 6, and H_{1} = lc(A_{2})Y^{2} = −Y^{2}, so the polynomial part must have degree 0. By formula (6) once more, d_{0} = max_{i=0,1,2}(deg(A_{i})) = 4, and H_{0} = lc(A_{2})Y^{2} = −Y^{2}, so any rational solution of (11) must be of the form (5) with P = 0. The P-adic part: From (12), Γ_{x}(S) = {1, 2}, so δ_{x}(S) = 2. By formula (9), n_{2} = max_{i=0,1,2}(2i − ν_{x}(A_{i})) = 2, and

H_{2} = ϕ_{x}(\frac{A_{2}}{x^{2}} )Y^{2} + ϕ_{x}(A_{1})Y = Y^{2} − Y

so we look for the solutions R of degree 0 in \overline{K}[x] of R^{2} − R ≡ 0 (mod x). Since deg(R) = 0, we set R = a ∈ \overline{K}, and solve a^{2} − a = 0 in \overline{K}. The 2 solutions are a = 0 and a = 1, so the 2 possible terms of order −2 at x = 0 are 0 and 1/x^{2}. We must look for the term of order −1 in both of those cases. Case a = 0: Setting δ = 1, we have by formula (9): n_{1} = max_{i=0,1,2}(i − ν_{x}(A_{i})) = 1, and H_{1} = ϕ_{x}(A_{1})Y = −Y, which has only 0 for solution, so the only possible P-adic part at x in this case is 0. Case a = 1: We must first do the change of variable u = 1/x^{2} + \overline{u}, and the Ricatti equation for \overline{u} is the one associated with

\overline{L} = A_{2}(∂_{x} + \frac{1}{x^{2}} )^{2} + A_{1}(∂_{x} + \frac{1}{x^{2}} ) + A_{0}

= −\frac{B_{2}∂_{x}^{2} + B_{1}∂_{x} + B_{0}}{x}

where B_{2} = x^{5} −x^{3}, B_{1} = 8x^{4} +2x^{3} −5x^{2} −x and B_{0} = 12x^{3} + 6x^{2} − 2x − 3. We first look for the term of order −1 at x of \overline{u}. From formula (9), n_{1} = max_{i=0,1,2}(i − ν_{x}(B_{i})) = 0, and H_{1} = ϕ_{x}(B_{1}/x)Y = −Y, so the the only possible P-adic part for \overline{u} at x is 0. The last remaining P-adic part for \overline{u} is at C_{1} = 1 − x^{2}, and we find that Γ_{1−x^{2}}({B_{0}, B_{1}, B_{2}}) = {1}, so δ_{1−x^{2}}({B_{0}, B_{1}, B_{2}}) = 1. By formula (9), n_{1} = max_{i=0,1,2}(i − ν_{1−x^{2}}(B_{i})) = 1, and

H_{1} = ϕ_{1−x^{2}}(\frac{B_{2}}{1 − x^{2}} )Y (Y + 2x) + ϕ_{1−x^{2}}(B_{1})Y

= −xY (Y + 2x) + (x + 3)Y
= −xY^{2} − (2x^{2} − x − 3)Y

so we look for the solutions R of degree 1 in \overline{K}[x] of −xR^{2} − (2x^{2} − x − 3)R ≡ 0 (mod 1 − x^{2}), which reduces to −xR^{2} + (x + 1)R ≡ 0 (mod 1 − x^{2}). Since deg(R) = 1, we set R = a + bx, a, b ∈ \overline{K}, plug in the above equation, reduce modulo 1 − x^{2} and we get the following system for a and b:

a + b − 2ab = −a^{2} − b^{2} + a + b = 0

which has the 2 solutions a = b = 0 and a = b = 1. So the P-adic expansion of \overline{u} at 1 − x^{2} is either 0 or (1+x)/(1−x^{2}) = 1/(1−x). So the singular part of \overline{u} is either 0 or 1/(1 − x), so the singular part of u is either 0, 1/x^{2} or 1/x^{2} + 1/(1 − x). Since u has no polynomial part, we are left with computing the Q′/Q part for each of those 3 cases. The logarithmic derivative part: We have 3 cases, depending on the potential singular part of u. Case 0: Then, since (10) has no nonzero rational solution, (11) has no rational solution. Case 1/x^{2}: Then, we check whether B_{2}y′′+B_{1}y′+B_{0} = 0, has nonzero rational solutions. This equation has a one-dimensional space of rational solutions generated by h = 1/(x^{4} − x^{3}), so the corresponding solution of (11) is

u_{1} = \frac{1}{x^{2}} + \frac{h′}{h} = −\frac{4x^{2} − 4x + 1}{x^{3} − x^{2}}.

Case 1/x^{2} + 1/(x − 1): Then we must do the change of variable \overline{u} = 1/(1 − x) + ũ, and the Ricatti equation for ũ is the one associated with

\overline{L} = B_{2}(∂_{x} + \frac{1}{1 − x})^{2} + B_{1}(∂_{x} + \frac{1}{1 − x}) + B_{0}

= (x^{5} − x^{3})∂_{x}^{2} + (6x^{4} − 5x^{2} − x)∂_{x} + 6x^{3} − 3x − 3.

This equation has a one-dimensional space of rational solutions generated by h = 1/x^{3}, so the corresponding solution of (11) is

u_{2} = \frac{1}{x^{2}} +\frac{1}{1 − x}+\frac{h′}{h} = \frac{1}{x^{2}} +\frac{1}{1 − x}−\frac{3}{x} = −\frac{4x^{2} − 4x + 1}{x^{3} − x^{2}}.

Since u_{1} = u_{2}, it is the only rational solution. The corresponding exponential solution of (10) is

y = e^{∫ u_{1}} = \frac{e^{−1/x}}{x^{3}(1 − x)}

and there is only a 1-dimensional vector space of exponential solutions. The second solution is found by using this one to reduce the order and solving the resulting first order equation.

We now give a few examples solved with the AXIOM implementation of this algorithm. Since this implementation is also used as a top-level differential equation solver, the following modification is made: whenever m independent rational solutions y_{1}, …, y_{m} are found, we use reduction of order to reduce equation (1) to an equation of order n − m. This new equation can still have nonzero rational solutions (which are not rational solutions of the initial equation), so we can repeat the search for rational solutions until we get to an equation with no rational solutions. The solutions found in this way are always Liouvillian, but some may be non-exponential over \overline{K} (because of the integration steps). The first two examples are from [5] and have complex singularitites which do not appear in the exponential solutions:

                                                  (x**4+x)*d(y x,x,2)+(x**3-1)*d(y x,x)-x**2*y x
                                                         4      ,,        3      ,       2
                                                  (5)  (x  + x)y  (x) + (x  - 1)y (x) - x y(x)
                                                  -> solve(%,y,x).basis
                                                     (6)
                                                         +------+
                                                        3| 3
                                                       [\|x  + 1 ,
                                                         +------+   x
                                                        3| 3      ++          %F
                                                        \|x  + 1  |   ------------------- d%F ]
                                                                 ++             +-------+
                                                                         3     3|  3
                                                                      (%F  + 1)\|%F  + 1
                                                  -> 4*(x**2+x+1)**2*d(y x,x,2) -3*y x
                                                          4     3      2           ,,
                                                  (11) (4x  + 8x  + 12x  + 8x + 4)y  (x) - 3y(x)
                                                  -> solve(%,y,x).basis
                                                         +----------+              +----------+
                                                         | 2               2x + 1  | 2
                                                  (12) [\|x  + x + 1 ,atan(------)\|x  + x + 1 ]
                                                                             +-+
                                                                            \|3

The highest-order example tried to date is from E.Tournier’s thesis ([12]):

                                                  -> eq

(2)
6 (vi)         5 (v)          4 (iv)
x y    (x) + 20x y   (x) + 122x y    (x)
+
3 ,,,          2 ,,
258x y   (x) + 156x y  (x)
+
,
8xy (x) - 8y(x)
-> solve(%,y,x).basis
+---+
1    log(x)   - \|- 1 log(x)
(3)  [--,x,------,%e              ,
2      2
x      x
+---+             2
\|- 1 log(x) log(x)  + log(x)
%e            ,----------------]
2
x

References

[1] S.A. Abramov, Rational Solutions of Linear Differential and Difference Equations with Polynomial Coefficients (in russian), Journal of Computational Mathematics and Mathematical Physics 29, No.11, 1611–1620 (1989).

[2] S.A. Abramov & K.Yu. Kvashenko, Fast Algorithms for the Search of the Rational Solutions of Linear Differential Equations with Polynomial Coefficients, in “Proceedings of ISSAC ’91”, Bonn, ACM Press, 267–270 (1991).

[3] M. Bronstein, The Risch Differential Equation on an Algebraic Curve, in “Proceedings of ISSAC ’91”, Bonn, ACM Press, 241–246 (1991).

[4] M. Bronstein, On Solutions of Linear Ordinary Differential Equations in their Coefficient Field, Journal of Symbolic Computation 13, No.4, 413– 439 (1992).

[5] E. Kamke, Differentialgleichungen, Lösungmethoden und Lösungen, I. Gewöhnliche Differentialgleichungen, Akad. Verlag., Leipzig, 1959.

[6] Kovacic J.J.: An Algorithm for Solving Second Order Linear Homogeneous Differential Equations, Journal of Symbolic Computation 2, 3–43 (1986).

[7] F. Schwarz, A Factorisation Algorithm for Linear Ordinary Differential Equations, in “Proceedings of ISSAC ’89”, Portland, ACM Press, 17–25 (1989).

[8] M.F. Singer, Liouvillian Solutions of n^{th} Order Homogeneous Linear Differential Equations, American Journal of Mathematics 103, No. 4, 661–682 (1981).

[9] M.F. Singer, Liouvillian Solutions of Linear Differential Equations with Liouvillian Coefficients, Journal of Symbolic Computation 11, 251–273 (1991).

[10] M.F. Singer & F. Ulmer, Bounds and Necessary Conditions for Liouvillian Solutions of (Third Order) Linear Differential Equations, manuscript 1992, abstract in these proceedings.

[11] C. Smith, A Discussion and Implementation of Kovacic’s Algorithm for Ordinary Differential Equations, Res. Report CS 84-35, Dpt. of Computer Science, Univ. of Waterloo, Ontario, 1984.

[12] E. Tournier, Solutions Formelles d’Equations Différentielles, Ths̀e d’état, IMAG, Grenoble, 1987.