6.1 Kovacic algorithm implementation

This paper gives a detailed overview and a number of worked out examples illustrating the Kovacic [1] algorithm for solving second order linear differential equation \(A(x) y''+ B(x) y' + C(x) y=0\) where \(A,B,C\) are rational functions with complex coefficients in the independent variable \(x\). All three cases of the algorithm were implemented in a software package based on an object oriented design and complete source code listing given in the appendix with usage examples. Implementation used the Maple computer algebra language.¹. This package was then used to analyze the distribution of Kovacic algorithm cases on \(3000\) differential equations.

Kovacic [1] gave an algorithm for finding a closed form Liouvillian² solution to any linear second order differential equation \(A y''+ B y' + C y=0\) if such a solution exists. Smith [2] gave an implementation based on a modified version of Kovacic algorithm by Saunders [3].

The current implementation is based on the original paper by Kovacic and uses the new object oriented features in Maple. The accompanied software package have been tested on \(3000\) differential equations with each solution verified using Maple’s odetest. The test suite is included as a separate module. The Appendix describes how to use the software.

The Kovacic algorithm finds one (basis) solution of \(A y''+ B y' + C y=0\). The second basis solution is found using reduction of order. The general solution is a linear combination of the two basis solutions found.

The algorithm starts by writing the input ode \(A y''+ B y' + C y=0\) as

\begin{align*} y''+ a y' + b y=0 \tag {1} \end{align*}

Where \(a=\frac {B}{A},b=\frac {C}{A}\). The substitution

\begin{align*} z &= y e^{\frac {1}{2}\int a\,dx} \tag {2} \end{align*}

is then applied to (1) which transforms it to a second order ode in the new dependent variable \(z(x)\) without the first derivative

\begin{align*} z''&=r z \tag {3} \end{align*}

\(r\) in the above is given by

\begin{align*} r &= \frac {1}{4} a^2 + \frac {1}{2} a' - b \tag {4} \end{align*}

It is ode (3) which is solved by the algorithm and not (1). Equation (3) will be called the DE from now on.

If a solution \(z(x)\) to the DE is found, then the first basis solution to the original ode is obtained using the transformation (2) in reverse

\begin{align*} y &= z e^{-\frac {1}{2}\int a\,dx} \end{align*}

The second solution is found using reduction of order.

These are the four possible cases to consider.

1. DE has solution \(z=e^{\int \omega dx}\) with \(\omega \in \mathbb {C}(x)\).
2. DE has solution \(z=e^{\int \omega dx}\) with \(\omega \) polynomial over \(\mathbb {C}(x)\) of degree \(2\).
3. Solutions of DE are algebraic over \(\mathbb {C}(x)\).
4. DE has no Liouvillian solution.

Before describing how the algorithm works, there are necessary (but not sufficient) conditions that determine which case the DE satisfies. Only those cases that meet the necessary conditions will be attempted.

The following are the necessary conditions for each case. To check each case, let \(r=\frac {s}{t}\) where \(\gcd (s,t) =1\). The order of \(r\) at \(\infty \) (from now on referred to as \(\mathcal {O}(\infty )\)) is defined as \(\deg (t)-\deg (s)\). The poles of \(r\) and the order of each pole need to be determined.

Knowing the order of the poles of \(r\) and \(\mathcal {O}(\infty )\) is all what is needed to determine the necessary conditions for each case. These conditions are the following

1. Case \(1\). Either no pole exists, or if a pole exists, the order must be either one or even. If \(\mathcal {O}(\infty )\) is less \(3\), then it must be even otherwise it can be even or odd.
2. Case \(2\). \(r\) must have at least one pole either of order 2 or odd order greater than 2. There are no conditions on \(\mathcal {O}(\infty )\).
3. Case \(3\). \(r\) must have a pole either of order \(1\) or \(2\). No other order is allowed. \(\mathcal {O}(\infty )\) must be at least \(2\).

If the conditions of a case are not satisfied then the case will be attempted as the algorithm guarantees that there will be no Liouvillian solution. However if the conditions are satisfied, this does not necessarily mean a solution exists. As an example \(y''=1/x^6 y\) satisfies only case one, but running the algorithm on case one shows that there is no Liouvillian solution.

The following table summarizes the above conditions for each case.

Case	Allowed pole order for \(r\)	Allowed value for \(\mathcal {O}(\infty )\)
1	\(\left \{ 0,1,2,4,6,8,\cdots \right \} \)	\(\left \{ \cdots ,-6,-4,-2,0,2,3,4,5,6,\cdots \right \} \)
2	Need to have at least one pole of order \(2\) or pole of odd order greater than \(2\). Any other pole order is allowed as long as the above condition is satisfied. The following are examples of pole orders which are allowed. \(\{1,2\}\),\(\{1,3\}\),\(\{2\}\),\(\{3\}\),\(\{3,4\}\),\(\{1,2,5\}\).	no conditions
3	\(\left \{ 1,2\right \} \)	\(\left \{ 2,3,4,5,6,7,\cdots \right \} \)

Table 6.1: Necessary conditions for each Kovacic case

Some observations: In case one, no odd order pole is allowed except for order 1. Case one is the only case that could have no pole in \(r\), which is the same as a pole of order zero. Case two and three require at least one pole. For case three, only poles of order \(1\) or \(2\) are allowed. If \(\mathcal {O}(\infty )\) is zero, then only possibility is either case one or two. For case one, if \(\mathcal {O}(\infty )\) is negative, then it must be even.

The above table also shows that when \(r\) has only one pole of order \(2\) and \(\mathcal {O}(\infty )\) equals \(2\) or higher then all three cases are possible. Also, if \(r\) has two poles one of order \(1\) and the other of order \(2\) and \(\mathcal {O}(\infty )\) equals \(2\) or higher then all three cases are possible.

These are the only two possibilities where all three cases have the same necessary conditions.

¹The complete Kovacic package in one mpl file accompany the arXiv version of this paper

²Wikipedia defines Liouvillian function as function of one variable which is the composition of a finite number of arithmetic operations \((+,-, \times , \div )\), exponentials, constants, solutions of algebraic equations (a generalization of nth roots), and antiderivatives. Kovacic in his original paper says “Such a solution may involve exponentials, indefinite integrals and solutions of polynomial equations. (As we are considering functions of a complex variable, we need not explicitly mention trigonometric functions, they can be written in terms of exponentials. Note that logarithms are indefinite integrals and hence are allowed.”