Solving Third Order Linear Differential Equations in Terms of Second Order Equations

Mark van Hoeij^{∗}

Dept. of Mathematics, Florida State University

Tallahassee, FL 32306-4510, USA

hoeij@math.fsu.edu

ABSTRACT

This paper presents a simplified version of a method by Michael Singer for reducing a third order linear ode to a second order linear ode whenever possible. An implementation is available as well.

Categories and Subject Descriptors

G.4 [Mathematical Software]: Algorithm design and analysis

General Terms

Algorithms

Keywords

Linear Differential Equations, Reduction of Order

1. INTRODUCTION

Let k be a differential field of characteristic 0, and let k[∂] be the ring of differential operators with coefficients in k (see Section 1.1 for more details on the notations). Let L ∈ k[∂] and let L(y) = 0 be the corresponding linear ode (ordinary differential equation). In [19] Michael Singer described in which situations L(y) = 0 has so-called eulerian solutions (defined in [19, Sect. 2]), which by [19, Thm. 4.3] is equivalent to having a non-zero solution that can be written in terms of solutions of linear ode’s of second order (with coefficients in \overline{k}, the algebraic closure of k). If L has order 3, we consider these three (not mutually exclusive) cases:

(1). L is the symmetric square of an operator L_{2} of order 2.

(2). L is reducible in \overline{k}[∂].

(3). L is gauge equivalent to a symmetric square in \overline{k}[∂].

∗Supported by NSF grant 0511544 Cases (1),(3) resp. (2) are part (d) resp. (a),(b),(c) in [19, Thm. 4.3]. If none of (1),(2),(3) holds then L is not solvable in terms of second order linear ode’s in the sense of [19].

A similar problem is to decide when L allows a non-zero solution that can be written in terms of solutions of second order linear ode’s whose coefficients are now restricted to k. For this restricted situation, replace cases (2),(3) by

(k2). L is reducible in k[∂].

(k3). L is gauge equivalent to a symmetric square in k[∂].

Detecting case (1) and finding the corresponding operator L_{2} is easy (see [19, Lemma 3.4] or Section 1.1). Case (1) remains the same if we restrict to coefficients in k, because if L is the symmetric square of a monic operator L_{2} ∈ \overline{k}[∂], then [19, Lemma 3.4] shows that L_{2} is in k[∂]. Formulas for testing if an n’th order operator equals the (n − 1)’th symmetric power of a second order operator are known as well, see for example [3], [15, Prop. 3.2], or [6, Capitel 5].

For case (k2) one has to compute a first order right or left factor of L in k[∂] if it exists. If k = ℂ(x) then there are several algorithms to do this (e.g. [17, Sect. 4.1] or [4] or the references therein) and several computer algebra systems have implementations to find such factors. If L is reducible in \overline{k}[∂] (case (2)) but irreducible in k[∂] then there are Liouvillian solutions, for which there exist algorithms, see for example [17, Sect. 4.3], [21], [23], or [13]. The latter algorithm, with a number of unpublished improvements, was implemented by Manuel Bronstein.

The symmetric square of a third order operator L has order 5 in case (1), and order 6 otherwise (Section 1.1 contains a proof, but it also follows from the proof of [23, Lemma 3.1] or from classical results [6, Capitel 4]). Assume now that cases (1),(2) do not hold and denote L_{6} as the symmetric square of L. Singer showed that then case (3) holds iff L_{6} has a first-order right-hand factor ∂ − r, where r is either in k, or algebraic of degree 3 over k. Moreover, case (k3) holds iff: (i) L_{6} has a first order right-hand factor in k[∂], and (ii) a certain conic (equation (4.2.1) in [19]) has a non-zero solution over k.

If k = ℂ(x) (or more generally, a C_{1}-field) then condition (ii) can be omitted because for such k, every conic over k is solvable over k. For other fields k, such as ℂ(x, e^{x}) or ℚ(x), to show that condition (ii) is actually necessary, one needs to show that conics without solutions over k can occur in this context (i.e. in equation [19, (4.2.1)] for some L ∈ k[∂] that satisfies condition (i)). Such examples were first constructed in [11], where it was shown explicitly that up to birational equivalence every conic over k does occur in this context. This also implies that to implement the reduction from order 3 to order 2, one can not avoid solving equation [19, (4.2.1)] or some conic equivalent to it. For this reason the author implemented an algorithm [10] for solving conics over k = ℚ(x), and more generally, any field F (t_{1}, …, t_{r}) for which a conic-solver over F is available. A similar algorithm for k = ℝ(x) was given earlier by Josef Schicho in [18].

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee.

ISSAC’07, July 29–August 1, 2007, Waterloo, Ontario, Canada.

Copyright 2007 ACM 978-1-59593-743-8/07/0007 ...$5.00.

An alternative approach to perform the reduction from order 3 to order 2 in case (3) was given in the proof of [11, Prop. 3.1], and in [14]. Here the conic is constructed in a different way, but like in [19], the reduction from order 3 to order 2 over k is only possible when this conic has a point defined over k. So this conic must be birationally equivalent to equation [19, (4.2.1)]. However, this does not imply that the two approaches are computationally equivalent, and indeed, examples (see [7]) show that the two conics do not have the same coefficients. In most examples tried so far, the coefficients in the [11] approach were larger. Therefore we will follow the method given in [19], with some algorithmic improvements explained below.

The expression I_{0} in [19, Prop. 4.2] is an exponential solution of L_{6} which means that it corresponds to a first order right-hand factor ∂ − r of L_{6} (here r = I′_{0}/I_{0} is in \overline{k} in case (3) and in k in case (k3)). It is one of the coefficients of the conic [19, (4.2.1)]. Like [19] we will compute this I_{0}. The main difference is how the other coefficients of the conic are computed. In [19] two of those coefficients (denoted I_{1} and I_{2}) are computed by solving two other 6’th order operators (denoted L_{a}^{◯ s 2} and L_{b}^{◯ s 2} in [19, Sect. 6]). Instead, we will give a formula for the conic that requires only the computation of I_{0}.

Thus, we only need to compute an exponential solution of one 6’th order ode instead of three, and then get the remaining coefficients for free^{1}. This means less work for the algorithm, and makes the algorithm much easier to implement because a number of technical problems (the second half of section 6 in [19]) are avoided this way (for instance, solving those two other 6’th order operators only determines I_{1} and I_{2} up to constant factors. Those constants are important but are not known a priori, causing complications we avoid.)

Comparison to prior work: The main problem in this paper is to reduce 3’rd order equations to 2’nd order equations whenever possible. This problem has already been solved theoretically in [19]. The main application is to solve 3’rd order equations, but it can also be useful for 4’th order equations; it was shown in [5, 16] that order 3 → 2 reduction can be used to solve similar problems for order 4 (this will be illustrated in Section 3). Despite these applications, the algorithm in [19] has not been implemented for two reasons. First, as already mentioned, there were a number of technical problems that made the algorithm difficult to implement. Second, since those problems were treated by solving systems of polynomial equations, it is likely that the resulting algorithm would have been too slow to be practical. Both issues are addressed in this paper, resulting in an algorithm that is both practical and easy to implement (see [7]). So the progress made here is not mathematical but algorithmic in nature.

1Similar progress was made by Bronstein [2] for the problem of factoring differential operators.

1.1 Preliminaries and Notations

This section will list notations and some facts about differential operators that will be needed in this paper. Proofs and more details can be found in [17] or [19].

Let k be a differential field of characteristic 0. A differential operator over k is an operator of the form L = a_{n}∂^{n} + ⋯ + a_{1}∂ + a_{0} that acts as L(y) = a_{n}y^{(n)} + ⋯ + a_{1}y′ + a_{0}y. If a_{n} ≠ 0 then the order of L is ord(L) = n. The set of differential operators over k forms a non-commutative ring k[∂], where multiplication is composition of operators.

The field of constants C_{k} is the set of elements of k with derivative 0. If C_{k} is algebraically closed, then there exists a differential field Ω, whose field of constants is again C_{k}, with the following property: For every L ∈ k[∂]−{0}, the set V (L) := {y ∈ Ω|L(y) = 0} is a C_{k}-vector space of dimension ord(L). Such Ω is called a universal extension, and V (L) is called the solution space of L.

An exponential solution of L over k is a non-zero u ∈ V (L) for which u′/u ∈ k. This corresponds to a first order righthand ∂ − r of L where r = u′/u. We will denote a non-zero solution u of ∂ − r as exp(∫ r).

Given L_{1}, L_{2} ∈ k[∂], the smallest order monic operator L ∈ k[∂] for which y_{1}y_{2} ∈ V (L) for every y_{1} ∈ V (L_{1}) and y_{2} ∈ V (L_{2}) will be denoted as L = L_{1} ◯ s L_{2}. This operator L is called the symmetric product in [19]. The operator L^{◯ s m} is called the m’th symmetric power of L, it is the symmetric product of m copies of L.

For any r ∈ k, the map S_{r} given by ∑ a_{i}∂^{i} ↦ ∑ a_{i}(∂−r)^{i} is an automorphism of the ring k[∂]. If L is monic then S_{r}(L) = L ◯ s (∂ − r).

If

L = ∂^{3} + a_{2}∂^{2} + a_{1}∂ + a_{0} ∈ k[∂]

then we can test if L is a symmetric square (case (1)) by following the proof of [19, Lemma 3.4], and one finds that L = L_{2}^{◯ s 2} for some second order operator L_{2} if and only if

a′′_{2} + 6a_{0} − 3a′_{1} + 2a_{2}(a′_{2} − a_{1} + 2a_{2}^{2}/9) = 0  (1)

in which case L_{2} must be in k[∂] because of the explicit formula (same as [19, (3.4.2)]):

L_{2} = ∂^{2} + \frac{a_{2}}{3} ∂ − \frac{1}{4} ( \frac{a′_{2}}{3} − a_{1} + \frac{2}{9} a_{2}^{2}) .  (2)

In general, the coefficients of L^{◯ s m} can be found by solving linear equations over k, see [19, Sect. 3] or [1]. For L^{◯ s 2} = A_{n}∂^{n} +…+A_{0} we get a homogeneous system of 6 equations in n + 1 unknowns A_{0}, …, A_{n}. If n = 6 then such system always has a non-zero solution. If n = 5 this system will only have a non-zero solution if its determinant vanishes, a condition that a computation shows to be equivalent with equation (1) and hence with case (1). So if we are not in case (1) then L^{◯ s 2} has order 6.

If we are in case (1) then L = L_{2}^{◯ s 2} with L_{2} ∈ k[∂] given above. Then L^{◯ s 2} = L_{2}^{◯ s 4} has order 5 by [19, Lemma 3.2(b)]. Thus, case (1) is very easy to detect (just check equation (1)) in which case the reduction from order 3 to order 2 is easy. From now on we will assume that L is not a symmetric square, or equivalently, that L^{◯ s 2} has order 6. The goal is now to find some gauge transformation R that sends L to another third order operator L_{R} with ord(L_{R}^{◯ s 2}) = 5. Then L_{R} must be the symmetric square of some L_{2}, so we can solve L_{R} (and hence L by computing the inverse gauge transformation) in terms of solutions of L_{2}.

1.2 Gauge transformations

Let L, R ∈ k[∂]−{0}, and apply operator R to the solution space of L. The result R(V (L)) is the solution space of another operator, that we denote as L_{R}, whose coefficients can be found by solving linear equations, see [19, Lemma 6.4].

Two operators L, L̃ ∈ k[∂] are called gauge equivalent if there exists an R ∈ k[∂] that bijectively maps the solutions of L to the solutions of L̃. In this case L̃ = L_{R}, and R is called a gauge transformation from L to L̃.

If R ∈ k[∂] is a gauge transformation from L to L̃, i.e., if R : V (L) → V (L̃) is bijection, then R maps no non-zero solutions of L to 0, which means that the greatest common right divisor (see [17, Sect. 2.1] for details) of R and L in k[∂] is 1. Then we can find R̃, T ∈ k[∂] with R̃R + T L = 1 with the extended Euclidean algorithm. The operator R̃R acts as the identity on V (L) since it is congruent to 1 modulo L. Hence R̃ is the inverse gauge transformation, i.e., the map R̃ : V (L̃) → V (L) is the inverse of R : V (L) → V (L̃).

As an example, let L = ∂^{2} − x (the Airy equation), and L̃ = ∂^{2} − \frac{1}{x} ∂ − x. These operators are gauge equivalent because L̃ = L_{R} where R = ∂. This means that by applying R to the solutions of L we get the solutions of L̃. The inverse gauge transformation is \frac{1}{x} ∂, it maps V (L̃) back to V (L).

2. REDUCING ORDER 3 TO ORDER 2

Assumptions: k is a differential field of characteristic 0, L_{3} ∈ k[∂] has order 3, and its solutions can be written in terms of solutions of second order linear ode’s in the sense of [19]. We will also assume that we are not in case (1) or (2) since there already are algorithms for those cases.

2.1 Computing a first order factor

The assumptions imply that we are in case (3), and that L_{3}^{◯ s 2} has a first order right-hand factor ∂−r where r is either in k, or k(r) is an algebraic extension of k of degree 3, see [19, Thm. 4.3(d)]. However, if [k(r) : k] = 3 then L_{3} must have Liouvillian solutions by [22, Thm. 4.7(i)(a)] and this implies case (2), contrary to our assumptions above, so only the case r ∈ k remains.

The same reasoning applies if we replace k by the coefficient field k′ ⊆ k of L_{3} (the smallest differential field for which L_{3} ∈ k′[∂]) and thus we only need to search for righthand factors ∂ − r of L_{3}^{◯ s 2} with r ∈ k′. This saves computation time, because if for example k′ = C(x) for some field of constants C, then to find r we only need the algorithm in [4, Sect. 5.4] and not the more complicated algorithm in [4, Sect. 6 and 7].

Another algorithmic improvement is the following. If the coefficient of ∂^{2} in L_{3} is 0 (applying S_{a_{2}/3} will accomplish this if this coefficient was a_{2}) then it follows from [22] that the cube of the exponential solution exp(∫ r) will be in k (the cubes of the characters are trivial for each of the groups in [22, Thm. 4.7(i)]). This means that if we use for example [4] to compute r, then we only need the exponents in \frac{1}{3} ℤ, and no other (generalized) exponents (see also [4, Sect. 7.2]).

Let L := L_{3} ◯ s (∂ +r/2). Then L^{◯ s 2} has a first order righthand factor (∂ −r) ◯ s (∂ +r/2)^{◯ s 2} = ∂, in other words, 1 is a solution of L^{◯ s 2}. Solving L is equivalent to solving L_{3} since their solution spaces differ only a factor exp(∫ r/2) from each other. We will proceed with L instead of L_{3} because that will lead to shorter formulas in Section 2.2.

2.2 Computing the conic

Write L = ∂^{3} + a_{2}∂^{2} + a_{1}∂ + a_{0} ∈ k[∂]. Given a non-zero R = b_{0} + b_{1}∂ + b_{2}∂^{2} ∈ k[∂] the linear map

R : V (L) → V (L_{R})

induces a linear map

R_{2} : V (L^{◯ s 2}) → V (L_{R}^{◯ s 2}).

Now V (L^{◯ s 2}) has dimension 6 because L is not a symmetric square. We want L_{R} to be a symmetric square, so the dimension of V (L_{R}^{◯ s 2}) should be 5. That is equivalent to R_{2} having a 1-dimensional kernel, which in turn corresponds to an exponential solution (see [17, Lemma 4.8 part 3]). In Section 2.1 we replaced L_{3} by L to make this exponential solution have value 1, so R_{2}(1) = 0. Now R depends linearly on the unknowns b_{0}, b_{1}, b_{2} and hence R_{2} depends quadratically on b_{0}, b_{1}, b_{2}. One can verify (see below) that R_{2}(1) = 0 reduces to

b_{0}^{2} + cb_{1}^{2} − 2cb_{0}b_{2} + c′b_{1}b_{2} + (c′′/2 + a_{1}c + a_{2}c′/2)b_{2}^{2} = 0 (3)

where

c := \frac{a′′_{0} + 7a′_{0}a_{2}/3 + a_{0}(4a_{2}^{2}/9 + 4a_{1} − a′_{2}/3)}{a′′_{2} + 6a_{0} − 3a′_{1} + 2a_{2}(a′_{2} − a_{1} + 2a_{2}^{2}/9)}  (4)

There are multiple equivalent ways to represent formulas (3) and (4) because there exists a differential relation between a_{0}, a_{1}, a_{2} (namely: L^{◯ s 2}(1) = 0). For instance, c′ can be replaced by 2(a_{0} − a_{2}c)/3 because they are equivalent modulo the differential relation for a_{0}, a_{1}, a_{2}. Note that the denominator of (4) is not 0, otherwise we would be in case (1), see equation (1) in Section 1.1.

That R_{2}(1) = 0 is equivalent (modulo the a_{0}, a_{1}, a_{2} relation) to (3) is easy to verify with a computer computation, but formulas (3),(4) can also be found by hand, as follows. Write

R(y)^{2} = (b_{0}y + b_{1}y′ + b_{2}y′′)^{2} =  (5)
Y_{1}b_{0}^{2} + 2Y_{2}b_{0}b_{1} + Y_{3}b_{1}^{2} + 2Y_{4}b_{0}b_{2} + 2Y_{5}b_{1}b_{2} + Y_{6}b_{2}^{2}

where

(Y_{1}, …, Y_{6}) := (y^{2}, yy′, (y′)^{2}, yy′′, y′y′′, (y′′)^{2}).

With L(y) = 0 we have y′′′ = −a_{2}y′′ − a_{1}y′ − a_{0}y which implies

Y′ = ( \begin{matrix}0 & 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 & 0 \\ −a_{0} & −a_{1} & 0 & −a_{2} & 1 & 0 \\ 0 & −a_{0} & −a_{1} & 0 & −a_{2} & 1 \\ 0 & 0 & 0 & −2a_{0} & −2a_{1} & −2a_{2} \end{matrix} ) Y

where Y is the transpose of (Y_{1}, …, Y_{6}). Setting Y_{1} = y^{2} to 1 (the exponential solution of L^{◯ s 2}), and denoting Y_{3} by c, then from rows 1, 2, 3, 5 in the above matrix we quickly find Y_{2} = yy′ = Y′_{1}/2 = 0, Y_{3} + Y_{4} = Y′_{2} = 0 so Y_{4} = −c, Y_{5} = Y′_{3}/2 = c′/2, and Y′_{5} = −a_{0}Y_{2} − a_{1}Y_{3} − a_{2}Y_{5} + Y_{6} so Y_{6} = c′′/2 + a_{1}c + a_{2}c′/2. Substituting these into (5) gives (3).

From row 4 we get c′ = −Y′_{4} = a_{0}Y_{1} + a_{1}Y_{2} + a_{2}Y_{4} − Y_{5} = a_{0} − a_{2}c − c′/2 which gives the relation c′ = 2(a_{0} − a_{2}c)/3 mentioned before. Row 6 expresses Y′_{6} in terms of Y_{1}, …, Y_{6}, spelling this out gives a linear differential equation for c of order 3. Since c′ can be written in terms of c, the same must also be true for c′′ and c′′′ (write c′′ = (c′)′ = (2(a_{0} − a_{2}c)/3)′ = 2(a′_{0} − a′_{2}c − a_{2}c′)/3, and replace c′ by 2(a_{0} − a_{2}c)/3. Then repeat for c′′′ = (c′′)′). This way we can eliminate the derivatives of c from the row 6 equation. What remains is a linear equation in c that gives us (4).

2.3 An example

Let k = ℂ(x) and

L_{3} = ∂^{3} − \frac{1}{x(x − 1)} ∂^{2} − \frac{4(4x^{2} − 2x + 1)}{x^{2}} ∂ + \frac{12}{x(x − 1)}

This example was constructed to have Bessel type solutions (see [7] how to make such examples) so it is solvable in terms of solutions of second order ode’s. Then at least one of cases (1),(2),(3) holds by [19].

L is irreducible in k[∂] and is not a symmetric square. There are no Liouvillian solutions (Bessel functions are not Liouvillian) so L is irreducible in \overline{k}[∂] as well. So only case (3) remains.

Let L_{6} = L_{3}^{◯ s 2} and compute its exponential solutions. We find the solution u = 1, so r = u′/u = 0 and L = L_{3} in this example. Let a_{i} be the coefficient of ∂^{i} in L. Then equations (4) and (3) become c = −12 and

b_{0}^{2} − 12b_{1}^{2} + 24b_{0}b_{2} + 0b_{1}b_{2} + \frac{48(4x^{2} − 2x + 1)}{x^{2}} b_{2}^{2} = 0.

To solve this equation, we can use [24], or complete the square by substituting b_{0} ↦ b_{0} +cb_{2} and b_{1} ↦ b_{1} −c′b_{2}/(2c) in (3) and then use [10]. Completing the square is simplified by the absence of a b_{0}b_{1} term in equation (3) (the coefficient 2Y_{2} of this term vanishes because in Section 2.1 we rescaled the exponential solution to Y_{1} = 1, and Y′_{1} = 2Y_{2}). In the example we get

b_{0}^{2} − 12b_{1}^{2} + \frac{48(x − 1)^{2}}{x^{2}} b_{2}^{2} = 0.

This equation can easily be solved by hand in this particular example, but in general one needs an algorithm to solve such equations. Our implementation [7] uses [10] and finds

(b_{0}, b_{1}, b_{2}) = (0, −2(x − 1), x).

Reversing the substitution, we get

(b_{0}, b_{1}, b_{2}) = (−12x, −2(x − 1), x).

Now let R = b_{2}∂^{2} + b_{1}∂ + b_{0}, make R monic, and we get

R = ∂^{2} − \frac{2(x − 1)}{x} ∂ − 12.

This maps V (L) to V (L_{R}) for some operator

L_{R} = ∂^{3} + A_{2}∂^{2} + A_{1}∂ + A_{0}

where

A_{2} = \frac{−3}{(x − 1)x}

and A_{0}, A_{1} are large expressions we will omit. If s ∈ k with s ≠ 0 then we can replace R by sR, and multiply all solutions of L_{R} by s, by replacing L_{R} with

S_{\frac{s′}{s}} (L_{R}) = ∂^{3} + (A_{2} − 3 \frac{s′}{s})∂^{2} + ⋯

We aim to minimize the coefficient of ∂^{2} by finding some s ∈ k for which A_{2}−3 \frac{s′}{s} is as small as possible. For k = ℂ(x) this can be done as follows. Apply the integration algorithm on \frac{1}{3} A_{2}, then take the logarithmic terms c_{i}log(f_{i}), and then let s be the product of f_{i}^{[c_{i}]} taken over those i for which c_{i} ∈ ℚ. Here [c_{i}] denotes c_{i} rounded to the nearest integer. In our example this leads to s = x/(x − 1), so we get

R = \frac{x}{x − 1} ∂^{2} − 2∂ − 12 \frac{x}{x − 1} .

Now recompute L_{R} for this new R and we find

L_{R} = ∂^{3} − \frac{8(2x − 1)}{x} ∂ − \frac{4}{x^{2}} .

We compute the inverse gauge transformation (Section 1.2) and find

R̃ = − \frac{1}{24} ∂ − \frac{1}{12} \frac{x − 1}{x}

which maps V (L_{R}) to V (L). Applying formula (2) to L_{R} yields

L_{2} = ∂^{2} − \frac{4x − 2}{x} .

Our implementation [7] follows the above steps to compute L_{2} and R̃. In the example, L_{2} has the following basis of solutions

y_{1} = x (I_{0}(2x) − I_{1}(2x)) , y_{2} = x (K_{0}(2x) + K_{1}(2x))

where I_{α} and K_{α} are the Bessel I and K functions. Then y_{1}^{2}, y_{1}y_{2}, y_{2}^{2} is a basis of V (L_{R}), so

R̃(y_{1}^{2}), R̃(y_{1}y_{2}), R̃(y_{2}^{2})

is a basis of V (L). We find

R̃(y_{1}^{2}) = \frac{x}{12} · (I_{0}(2x) − I_{1}(2x)) · (xI_{0}(2x) − (x + 1)I_{1}(2x))

and similar expressions for the other two solutions of L.

2.4 The algorithm

Algorithm: ReduceOrder 3 → 2

Input: The field k, a monic third order L_{3} ∈ k[∂] which is assumed to be irreducible in \overline{k}[∂], and a boolean B.

Output: L_{2}, R̃ in k[∂] (or in \overline{k}[∂] if B = true) where R̃ bijectively maps V (L_{2}^{◯ s 2}) to V (L) if such L_{2}, R̃ exist.

1. L_{6} := L_{3}^{◯ s 2}

2. If ord(L_{6}) = 5 then let R̃ := 1 and L_{2} as in equation (2). Then return L_{2}, R̃ and stop.

3. Find r ∈ k for which ∂ − r is a right-hand factor of L_{6}. If such r does not exist, return “Order reduction not possible” and stop.

4. L := S_{−r/2}(L_{3}), see Section 1.1 for notations.

5. Find, if it exists, a non-zero solution (b_{0}, b_{1}, b_{2}) over k of (3) where a_{0}, a_{1}, a_{2} are the coefficients of L.

6. If no such solution exists, then do the following: If B = true (field extensions are allowed) then set (b_{0}, b_{1}, b_{2}) := (\sqrt{−c}, 1, 0) where c is as in (4), otherwise return “Order reduction not possible” and stop.

7. Let R := S_{r/2}(b_{0} + b_{1}∂ + b_{2}∂^{2}).

8. Compute (L_{3})_{R}, see Section 1.2. If possible, try as in Section 2.3 to find some s ∈ k such that replacing R by sR and updating (L_{3})_{R} accordingly leads to a smaller expression (L_{3})_{R}.

9. Use (2) to find L_{2} for which L_{2}^{◯ s 2} = (L_{3})_{R}.

10. Compute R̃, the inverse transformation of R, by solving (e.g. with the extended Euclidean algorithm) the equation R̃R + T L_{3} = 1.

11. Return L_{2}, R̃ and stop.

An implementation of this algorithm can be found on [7]. In step 6, the expression (\sqrt{−c}, 1, 0) is always a solution of (3), however, if at all possible we want to perform the reduction without introducing field extensions that could complicate the problem of solving L_{2}. After all, to get explicit solutions for L_{3} we still need to solve L_{2}.

The gauge transformation b_{0} + b_{1}∂ + b_{2}∂^{2} in step 7 sends L to a symmetric square. Then the gauge transformation S_{r/2}(b_{0} + b_{1}∂ + b_{2}∂^{2}) must send L_{3} to a symmetric square because L_{3} = S_{r/2}(L) and S_{r/2} is an automorphism of the ring of differential operators.

Selecting a different point on conic (3) will lead to a different output. Two different L_{2}’s, say L_{2a} and L_{2b}, in the output of the same L_{3} must be projectively equivalent (see Theorem 4.7 in [11]) which means that L_{2a} is gauge equivalent to S_{t}(L_{2b}) for some t. That implies that the solutions of each one can be expressed in terms of solutions of the other. Even so, it is still possible that a computer algebra system like Maple or Mathematica would solve one of L_{2a}, L_{2b} but not the other. In fact, in our experiments that turned out to be very common. Because of this, there is a strong incentive to try to find the “easiest” output L_{2}. For this reason our implementation computes not one point in step 5 but several points (see also Sections 2.2 and 6 in [10]). It then chooses the smallest one because that significantly increases the likelihood that L_{2} can be solved with current computer algebra systems.

Potential improvements: Before calling the above algorithm, it makes sense to compute local data (the p-curvature and the generalized exponents, see [4]) and check if this data is compatible with case (3). This way one can avoid unnecessary calls to this algorithm. The generalized exponents can be re-used to speed up step 3 using [12].

3. EXAMPLES OF ORDER 4

It was shown in [5, 16] that order 3 → 2 reduction can be used to solve similar problems for order 4. We will illustrate that with two examples. Given some irreducible L_{4} ∈ k[∂] of order 4, consider the following questions:

q1 Does there exist some second order L_{2} ∈ k[∂] such that L_{4} = L_{2}^{◯ s 3}?

q2 Do there exist second order operators L_{a}, L_{b} ∈ k[∂] such that L_{4} = L_{a} ◯ s L_{b}?

q3 Does there exist some second order L_{2} ∈ k[∂] such that L_{4} is gauge equivalent to L_{2}^{◯ s 3}? q4 Do there exist second order operators L_{a}, L_{b} ∈ k[∂] such that L_{4} is gauge equivalent to L_{a} ◯ s L_{b}?

Questions q1 resp. q2 are solved in [3] resp. [9], while q3 and q4 are illustrated below.

3.1 An example for q3

Let k = ℂ(x) and

L_{4} = ∂^{4} − \frac{2}{x} ∂^{3} − 10x∂^{2} − 24∂ + \frac{68 + 9x^{3}}{x} ∈ k[∂].

This L_{4} is irreducible in k[∂] (use a differential factorization algorithm, e.g. [5, Appendix A], [17, Sect. 4.2] or [8]). With [3] resp. [9] one can show quickly that we are not in case q1 resp. q2.

The method in [16, Chap. 4, Prop. 7] is to compute L^{◯ s 2}, and then to find a third order right-hand factor, if one exists. One finds

L_{3} = ∂^{3} − \frac{3x^{2}}{x^{3} − 7} ∂^{2} − \frac{x(4x^{3} − 49)}{x^{3} − 7} ∂ + \frac{2(5x^{3} − 56)}{x^{3} − 7} .

Note that this step is computationally expensive because L_{4}^{◯ s 2} is large (it has order 10) and because computing factors of order 3 takes more time than computing factors of order 1. Next is to apply the algorithm ReduceOrder 3 → 2. We find

L_{2} = ∂^{2} − x

(the Airy equation) so L_{2} and L_{3} (and L_{4} as we shall see) can be solved in terms of Airy functions. Then compute a gauge transformation from L_{2}^{◯ s 3} to L_{4} (in Maple this can be done with DEtools[Homomorphisms]). We find

R = x∂^{2} + ∂ − x^{2}.

So R(y) = xy′′ + y′ − x^{2}y is a solution of L_{4} for any solution y of L_{2}^{◯ s 3}. Hence

R(y_{1}^{i}y_{2}^{3−i}), i ∈ {0, 1, 2, 3}

is a basis of solutions of L_{4}, where y_{1}, y_{2} are the Airy functions (y_{1}, y_{2} is a basis of V (L_{2})).

3.2 An example for q4

Let k = ℂ(x) and

L_{4} = ∂^{4} + \frac{3 − 2x^{2}}{8x^{2}} ∂^{2} + \frac{3}{4x^{3}} ∂ − \frac{567}{256x^{4}} .

L_{4} is irreducible, and cases q1, q2, q3 do not hold.

The method in [16, Chap. 4, Prop. 10] or [14, Prop. 6.1] means computing a first order factor of L_{4}^{◯ s 2} (which is less expensive than computing a third order factor as was done in Section 3.1). After this, a quadratic equation in 4 unknowns needs to be determined and solved. The latter is not implemented, so at the moment it is easier to follow a strategy by Compoint and Weil which goes as follows. Let L_{6} be the second exterior power of L_{4} and determine its two right-hand factors of order 3 (DFactorLCLM in Maple). We find

L_{3,a} = ∂^{3} − \frac{2}{x} ∂^{2} − \frac{x^{2} − 7}{4x^{2}} ∂ − \frac{7}{4x^{3}}

and

L_{3,b} = ∂^{3} − \frac{2}{x} ∂^{2} − \frac{x^{2} − 12}{4x^{2}} ∂ − \frac{3}{x^{3}} .

Next is to apply algorithm ReduceOrder 3 → 2 on these two inputs. We find

L_{2,a} = ∂^{2} − \frac{x^{2} + 21}{16x^{2}} , L_{2,b} = ∂^{2} − \frac{x^{2} + 12}{16x^{2}} .

Now let L_{ab} := L_{2,a} ◯ s L_{2,b}. Then L_{ab} ◯ s (∂ − r) should be gauge equivalent to L_{4} for some r ∈ k. In this example L_{ab} itself is already gauge equivalent to L_{4}; Maple finds a gauge transformation

R = 4x^{2}∂^{3} + 9x∂^{2} − (x^{2} + \frac{21}{4} )∂ + \frac{75}{16x} − \frac{9x}{4} .

Then, for any solution y_{a} for L_{2,a} and y_{b} for L_{2,b} the expression R(y_{a}y_{b}) is a solution of L_{4}, and a basis of solutions of L_{4} is obtained this way. For example, if y_{a} = x^{1/2}I_{5/4}(x/4) and y_{b} = x^{1/2}I_{1}(x/4), where I_{α} is the Bessel I function, one finds this solution R(y_{a}y_{b}) =

\frac{9x^{2}}{32} (I_{0}(x/4)I_{1/4}(x/4) − I_{1}(x/4)I_{5/4}(x/4)) ∈ V (L_{4}).

4. FUTURE WORK

Methods for order higher than 3 can be found in [19, 20, 16, 14, 5]. Section 3 illustrates two of these methods on explicit examples of order 4, but improvements might be possible. For instance, in Section 3.1 one computes a 3’rd order factor of an operator of order 10, a time consuming step. However, the second exterior power of L_{4} has order 6 and a factor of order 1 that would take much less time to compute. This raises the question if this could be used to give a more efficient method. Whether Section 3.2 is optimal or not is not yet clear either.

The paper [14] gives numerous explicit constructions that can be implemented (software for solving systems, such as the ISOLDE package by Barkatou and Pflügel, is useful here). Some of these cases need a generalization of the conic solver.

It makes sense to treat the lower order cases first, because for high order it becomes harder to obtain an efficient implementation, and high order might also be less common in practice. So for future work, the next natural step would be to see if the order 4 cases can be improved, and to implement some of the higher order cases found in the above mentioned references.

5. REFERENCES

[1] M. Bronstein, T. Mulders, J-A. Weil, On symmetric powers of differential operators, ISSAC’1997 Proceedings, ACM, New York, 156-163, (1997).

[2] M. Bronstein, An improved algorithm for factoring linear ordinary differential operators, ISSAC’1994 Proceedings, ACM, New York, 336-340, (1994).

[3] R. Chalkley, Relative invariants for homogeneous linear differential equations, J. Differential Equations, 80, 107-153, (1989).

[4] T. Cluzeau, M. van Hoeij, A Modular Algorithm to Compute the Exponential Solutions of a Linear Differential Operator, J. Symb. Comput., 38, 1043-1076, (2004).

[5] E. Compoint, J-A. Weil, Absolute reducibility of differential operators and Galois groups, J. Algebra, 275, 77-105, (2004).

[6] G. Fano, Ueber lineare homogene Differentialgleichungen mit algebraischen Relationen zwischen den Fundamentallosungen, Math. Ann., 53, 493-590, (1900).

[7] M. van Hoeij, An implementation and examples for algorithm ReduceOrder

www.math.fsu.edu/~hoeij/files/ReduceOrder

[8] M. van Hoeij, Factorization of Differential Operators with Rational Functions Coefficients, J. Symb. Comput., 24, 537-561, (1997).

[9] M. van Hoeij, Decomposing a 4’th order linear differential equation as a symmetric product, Banach Center Publications, 58, 89-96, (2002).

[10] M. van Hoeij, J. Cremona, Solving conics over function fields, accepted for publication in JNTB.

www.math.fsu.edu/~hoeij/files/ConicProgram

[11] M. van Hoeij, M. van der Put, Descent for differential modules and skew fields. J. of Algebra, 296, 18-55 (2006).

[12] M. van Hoeij, J-A. Weil, An algorithm for computing invariants of differential Galois groups, J. Pure Appl. Algebra, 117-118, 353-379 (1997).

[13] M. van Hoeij, J-F. Ragot, F. Ulmer, J-A. Weil, Liouvillian solutions of linear differential equations of order three and higher, J. Symb. Comput., 28, 589-609 (1999).

[14] K. A. Nguyen, M. van der Put, Solving linear differential equations, preprint, (2006).

[15] W. R. Oudshoorn, M. van der Put, Lie symmetries and differential Galoisgroups of linear equations, Math. Comput., 71, 349-361, (2002).

[16] A. C. Person, Solving Homogeneous Linear Differential Equations of Order 4 in Terms of Equations of Smaller Order, PhD thesis, http://www.lib.ncsu.edu /theses/available/etd-08062002-104315/ (2002).

[17] M. van der Put, M. F. Singer, Galois Theory of linear Differential Equations, Grundlehren der mathematischen Wissenschaften, 328, Springer (2003).

[18] J. Schicho, Rational parametrization of real algebraic surfaces, ISSAC’1998 Proceedings, ACM, New York, 302-308, (1998).

[19] M. F. Singer, Solving Homogeneous Linear Differential Equations in Terms of Second Order Linear Differential Equations, Am. J. of Math., 107, 663-696, (1985).

[20] M. F. Singer, Algebraic Relations Among Solutions of Linear Differential Equations: Fano’s Theorem, Am. J. of Math., 110, 115-143, (1988).

[21] M. F. Singer, F. Ulmer, Linear Differential Equations and Products of Linear Forms, J. of Pure and Appl. Algebra, 117-118, 549-563, (1997).

[22] M. F. Singer, F. Ulmer, Galois Groups of Second and Third Order Linear Differential Equations J. Symb. Comput., 16, 9-36, (1993).

[23] F. Ulmer, Liouvillian solutions of third order differential equations, J. Symb. Comput., 36, 855-889, (2003).

[24] C. van de Woestijne, Surface Parametrisation without Diagonalisation, ISSAC’2006 Proceedings, ACM, New York, 340-344, (2006).