Solving Higher Order ODEs using Mathematica 5.2

Devendra Kapadia, Wolfram Research, Inc.

1. Introduction

The Mathematica function DSolve has been equipped with several modern algorithms for solving higher order linear ordinary differential equations (ODEs) in Version 5.2. The aim of this notebook is to explain the motivation for these developments and to provide some information and examples which illustrate the new functionality.

In Mathematica 5.1, we had focussed on adding methods for solving first order and second order ODEs such as Abel equations, hypergeometric-type equations and  equations with non-rational coefficients using DSolve. As explained in the Advanced Documentation for DSolve, the code structure for this function is hierarchical, so that the problem of solving  ODEs of order greater than 2 is often reduced to that of solving a first order or second order ODE. Within the last few years, a deeper understanding of several aspects of higher order ODEs (such as factorization techniques) has emerged which makes it possible to carrry out this reduction in a systematic way. Also, higher order ODEs (particularly orders 3 and 4) are increasingly  being seen in scientific models. Thus, we were interested in widening the application of the methods implemented in Version 5.1 to higher order ODEs.

The plan for the notebook is as follows.

In Section 2, we will review the methods for solving higher order ODEs which were already available in V 5.1. This will serve dual purposes.  First, it will provide a catalog of examples which can be solved by these methods. It will also allow us to discuss the important class of rational-exponential solutions which are used heavily in the new algorithms.

In Section 3, we will discuss the implementation of the Bronstein-Mulders-Weil-van Hoeij algorithm for solving linear ODES of arbitrary order that are symmetric powers of second order ODEs. Thus, given any linear second order ODE (say, Legendre's equation) whose general solution has a basis {u, v}, one can solve a higher order ODE of order n+1 whose solutions are {(u^n),  (u^(n - 1)) v,..., (v^n)} ('symmetric powers" of solutions of the second order ODE). The speed of these evaluations is impressive and most examples finish in less than 10 seconds.

Section 4 deals with a generalization of the notion of symmetric power in which we start with a pair of second order ODEs with bases {(u^1), (u^2)} and {(v^1), (v^2)} for their general solutions. An algorithm due to van Hoeij allows us to solve a fourth order ODE with {(u^1) (v^1),  (u^1) (v^2), (u^2) 
          (v^1) ,  (u^2) (v^2)} as a basis for its general solution. This fourth order ODE is called the symmetric product of the given second order ODEs. Once again, we present examples involving well-known second order ODEs such as Bessel's equation and Airy's equation to illustrate the algorithm.

In Section 5, we deal with the important notion of factorization for a differential operator. This simply means the process of writing the operator as a composition of lower order operators (for instance, a third order operator can sometimes be factored in terms of a second order operator and a first order operator respectively). While an algorithm for factorization has existed for more than a century, it is only during the last decade that efficient implementations, suitable for a computer algebra system such as Mathematica, have become available. We have implemented the Schwarz-Bronstein algorithm for factorization and will give examples which show that the implementation works well for ODEs of order up to 5.

We end this introduction by noting that it will be convenient to switch back and forth between (homogeneous) differential equations and the  corresponding differential operators since the algorithms really refer to the differential operators and the final step of integrating lower order ODEs to find the solutions  is straightforward in all cases.

2. Review of functionality available in Version 5.1

In Mathematica 5.1, DSolve  can solve the following types of higher order linear ODEs.

Equations with constant coefficients

This equation has constant coefficients and the solution is a combination of exponentials, as in the second-order case.

$Version

DSolve[y'''[x] - 4*y''[x] - 25*y'[x] + 28*y[x] == 0, 
      y[x], x]

Cauchy - Euler equations

A Cauchy-Euler equation is an ODE of the form

(x^n)  (y^((n) ))(x) +(Subscript[a, 1 ]) (x^(n - 1)) (y^((n - 1)))(x) +(Subscript[a, 2 ]) (x^(n - 2)) (y^((n - 2)))(x) ... +(Subscript[a, n]) y(x)==0.

Here is an equation of Cauchy-Euler type. It is solved by transformation to an ODE with  constant coefficients.

DSolve[
      x^3*y'''[x]  - 2*x^2*y''[x]  - x*y'[x]  + y[x] == 
        0,  y[x], x]

Legendre Linear equations

The following fourth order ODE is an example of a Legendre linear equation, which is a generalization of Cauchy-Euler equations (we allow (((a x + b))^n) instead of (x^n) as in the Cauchy - Euler case).

DSolve[(3  x + 5)^4*y''''[x]  + (3  x + 5)^3*y'''[x] - 
          2*(3  x + 5)^2*y''[x]  + 7 (3  x + 5)*y'[x]  + 
          y[x] == 0,  y[x], x]

Exact ODEs

Here is an example of an exact  ODE. The condition for exactness is shown to be satisfied below.

a0  =  1; a1 = -1; a2 = 5  x; a3 = 5;

ExactODE = 
      a0*y'''[x]  + a1*y''[x] + a2*y'[x] + a3*y[x]

conditionforexactness =  -D[a0,  {x, 3}] + D[a1, {x, 2}] - 
        D[a2, x] + a3

sol = DSolve[ExactODE == 0,  y, x] // Timing

Special Product Solutions involving Airy and Bessel functions

The solution of this equation has  a nice symmetrical structure . In fact, this is a special case of symmetric powers which has worked since V 4.1.

Timing[
      sol = DSolve[
          x^3*y'''[x]  + (4  x^3 - 3  x)*y'[x]  + 
              3  y[x] == 0,  y, x]]

This plot shows the oscillatory behavior of the solution for a particular choice of the arbitrary constants.

Plot[y[x] /. sol[[1]] /. {C[1] -> 2, C[2] -> 3, 
            C[3] -> 1} // Evaluate, {x,  0.5, 20}, 
      PlotStyle -> {Red}]

Rational-Exponential Solutions with reduction of order

The following ODE is solved using an algorithm for finding  “rational-exponential” solutions, which are combinations of rational functions and exponentials of the integrals of rational functions. Once the rational-exponential solution (e^x) has been computed, reduction of order (from three to two)  is used to  produce a complete solution for the given ODE in terms of exponential and Legendre functions.

DSolve[-((30 y[x])/7)+((30 y'[x])/7)+2 x y'[x]-y''[x]-2 x y''[x]+(x^2) y''[x]+Derivative[3][y][x]-(x^2) Derivative[3][y][x]==0, y,x]//Timing

Hypergeometric Type Solutions

A variety of transformations are used in DSolve to try and reduce a given ODE to one of hypergeometric type.

The solution to this linear ODE is expressed in terms of HypergeometricPFQ.

DSolve[
      x^2*y'''[x] + 
          3  x*y''[x] + (2 - x)*y'[x] + (4/5) y[x] == 
        0, y, x]

x^2*y'''[x] +  
              3  x*y''[x]  + (2 - x)*y'[x] + (4/5) 
                y[x] /. %[[1]] /.  {x ->  Random[Complex]} // 
        Simplify // Chop

3. Symmetric Powers of Second Order Differential Operators (new in 5.2)

As noted in the introduction, if {u, v} is a basis for the solution space  of a homogeneous linear second order ODE, then the second symmetric power of this equation is a linear third order ODE with  solution basis given by {(u^2), u v, (v^2)}.  More generally, the nth symmetric power of the second order ODE is a linear (n+1)th order equation with solution basis given by {(u^n),  (u^(n - 1)) v,..., (v^n)}.

Thus, if a higher order linear ODE is a symmetric power of a linear second order ODE, then its general solution has a very simple structure. Conversely, starting from any classical example of a linear second order ODE (Bessel, Airy, Legendre, etc.), we can construct equations of arbitrary high order whose solution depends in a simple way on the solutions of the second order equation.

Given an nth order linear ODE, the Bronstein-Mulders-Weil-vanHoeij algorithm for solving symmetric powers proceeds as follows:

1. Assuming that this is indeed a symmetric power, the coefficients of the corresponding second order ODE are calculated in terms of the coefficients of the input ODE, using a criterion due to van Hoeij [8].

2. The  (n-1)th symmetric power of the second order ODE is calculated using an iterative method due to Bronstein, Mulders and Weil [2]. Computationally, this is the most intensive step but due to caching of  derivatives, subsequent evaluations involving the same second order ODE finish quite rapidly. We then check that the resulting symmetric power is identical with the input ODE.

3. We calculate the general solution of the second order ODE and use it to write down the general solution of the input ODE.

Here are some examples which work with the new implementation.

Order 3  - Hypergeometric2F1

Timing[DSolve[(((2 a b)/((1 - x) (x^2))) - ((2 a b)/((((1 - x))^2) x)) - ((4 a b (c - (1 + a + b) x))/((((1 - x))^2) (x^2)))) y[x]+(((-1 - a - b)/((1 - x) x)) - ((4 a b)/((1 - x) x)) - ((c - (1 + a + b) x)/((1 - x) (x^2))) + ((c - (1 + a + b) x)/((((1 - x))^2) x)) + ((2 (((c - (1 + a + b) x))^2))/((((1 - x))^2) (x^2)))) y'[x]+((3 (c - (1 + a + b) x) y''[x])/((1 - x) x))+Derivative[3][y][x]==0,y,x]]

Order 4 - HypergeometricU and LaguerreL

Timing[DSolve[(((6 n)/(x^3)) - ((21 n (1 + a - x))/(x^3)) + ((18 n (((1 + a - x))^2))/(x^3)) - ((6 n)/(x^2)) + ((9 (n^2))/(x^2))) y[x]+(((2 (1 + a - x))/(x^3)) - ((7 (((1 + a - x))^2))/(x^3)) + ((6 (((1 + a - x))^3))/(x^3)) + 
                    (2/(x^2)) - ((10 n)/(x^2)) - ((7 (1 + a - x))/(x^2)) + ((30 n (1 + a - x))/(x^2))) y'[x]+(-((4 (1 + a - 
                                x))/(x^2)) + ((11 (((1 + a - x))^2))/(x^2)) - (4/x) + ((10 n)/x)) y''[x]+((6 (1 + a - x) Derivative[3][y][x])/x)+Derivative[4][y][x]==0, y,x]]

Order 7 - LegendreP and LegendreQ

Timing[sol=DSolve[(-((192 n (1 + 
                                  n) (x^5))/(((1 - (x^2)))^6)) + ((2496 n (1 + n) (x^3))/(((1 - (x^2)))^5)) - ((8160 (n^2) (((1 + n))^2) (x^3))/(((1 - (x^2)))^5)) - ((1632 n (1 + n) x)/(((1 - (x^2)))^4)) + ((7824 (n^2) (((1 + n))^2) x)/(((1 - (x^2)))^4)) - ((6912 (n^3) (((1 + n))^3) x)/(((1 - (x^2)))^4))) y[x]+(((64 (x^6))/(((1 - (x^2)))^6)) - ((1824 (x^4))/(((1 - (x^2)))^5)) + ((8256 n (1 + n) (x^4))/(((1 - (x^2)))^5)) + ((2880 (x^2))/(((1 - (x^2)))^4)) - ((19200 n (1 + n) (x^2))/(((1 - (x^2)))^4)) + ((22896 (n^2) (((1 + n))^2) (x^2))/(((1 - (x^2)))^4)) - 
                      (272/(((1 - (x^2)))^3)) + ((2208 n (1 + n))/(((1 - (x^2)))^3)) - ((4384 (n^2) (((1 + n))^2))/(((1 - (x^2)))^3)) + ((2304 (n^3) (((1 + n))^3))/(((1 - (x^2)))^3))) y'[x]+(-((2016 (x^5))/(((1 - (x^2)))^5)) + ((9408 (x^3))/(((1 - (x^2)))^4)) - ((19488 n (1 + n) (x^3))/(((1 - (x^2)))^4)) - ((3696 x)/(((1 - (x^2)))^3)) + ((12768 n (1 + n) x)/(((1 - (x^2)))^3)) - ((9408 (n^2) (((1 + n))^2) x)/(((1 - (x^2)))^3))) y''[x]+(((4816 (x^4))/(((1 - (x^2)))^4)) - ((7168 (x^2))/(((1 - (x^2)))^3)) + ((9632 n (1 + n) (x^2))/(((1 - (x^2)))^3)) + 
                      (616/(((1 - (x^2)))^2)) - ((1456 n (1 + n))/(((1 - (x^2)))^2)) + ((784 (n^2) (((1 + n))^2))/(((1 - (x^2)))^2))) Derivative[3][y][x]+(-((2800 (x^3))/(((1 - (x^2)))^3)) + ((1400 x)/(((1 - (x^2)))^2)) - ((1400 n (1 + n) x)/(((1 - (x^2)))^2))) Derivative[4][y][x]+(((560 (x^2))/(((1 - (x^2)))^2)) - 
                      (70/(1 - (x^2))) + ((56 n (1 + n))/(1 - (x^2)))) Derivative[5][y][x]-((42 x Derivative[6][y][x])/(1 - 
                    (x^2)))+Derivative[7][y][x]==0, y,x]]

Order 10 - HermiteH and Hypergeometric1F1

Timing[DSolve[(5552640 n - 36006912 (n^2) + 82705536 (n^3) - 
                    80797824 (n^4) + 28576800 (n^5) - 305593344 n (x^2) + 
                    1573323264 (n^2) (x^2) - 2536095744 (n^3) (x^2) + 
                    1298716416 (n^4) (x^2) + 2061656064 n (x^4) - 
                    7521942528 (n^2) (x^4) + 6401576448 (n^3) (x^4) - 
                    3630292992 n (x^6) + 7073796096 (n^2) (x^6) + 
                    1672151040 n (x^8)) y[x]+(-7796768 x + 130077728 n x - 
                    488359872 (n^2) x + 667149120 (n^3) x - 
                    304479648 (n^4) x + 133161216 (x^3) - 
                    1794152448 n (x^3) + 4758842880 (n^2) (x^3) - 
                    3422352384 (n^3) (x^3) - 497981952 (x^5) + 
                    4852818432 n (x^5) - 6849722880 (n^2) (x^5) + 
                    573272064 (x^7) - 3058311168 n (x^7) - 
                    185794560 (x^9)) y'[x]+(1041664 - 10993664 n + 33934912 (n^2) - 
                    40822848 (n^3) + 16915536 (n^4) - 55211488 (x^2) + 
                    465523168 n (x^2) - 1010674368 (n^2) (x^2) + 
                    638054208 (n^3) (x^2) + 357351168 (x^4) - 
                    2151553536 n (x^4) + 2473092864 (n^2) (x^4) - 
                    600897024 (x^6) + 1950423552 n (x^6) + 
                    262803456 (x^8)) y''[x]+(7089360 x - 46738560 n x + 87695520 (n^2) x - 
                    49769280 (n^3) x - 94726560 (x^3) + 442200000 n (x^3) - 
                    437279040 (n^2) (x^3) + 247476480 (x^5) - 
                    616619520 n (x^5) - 150105600 (x^7)) Derivative[3][y][x]+(-285440 + 1594560 n - 2666400 (n^2) + 
                    1382480 (n^3) + 11759280 (x^2) - 46231680 n (x^2) + 
                    40587360 (n^2) (x^2) - 52762560 (x^4) + 
                    109934880 n (x^4) + 46315520 (x^6)) Derivative[4][y][x]+(-690480 x + 2383920 n x - 1896048 (n^2) x + 
                    6360480 (x^3) - 11575872 n (x^3) - 8618400 (x^5)) Derivative[5][y][x]+(15456 - 48048 n + 35112 (n^2) - 437136 (x^2) + 
                    713328 n (x^2) + 1012368 (x^4)) Derivative[6][y][x]+(15960 x - 23760 n x - 75600 (x^3)) Derivative[7][y][x]+(-240 + 330 n + 3480 (x^2)) Derivative[8][y][x]-90 x Derivative[9][y][x]+Derivative[10][y][x]==0,y,x]]

Kamke 3.20

DSolve[y'''[x]  - 6  x*y''[x]  + 
          2 (4  x^2 + 2 a - 1)*y'[x]  - 8 a x y[x] == 0,  
      y, x]

Kamke 3.44

DSolve[
      x^2*y'''[x]  - 
          3  x*(x - m)*y''[
              x]  + (2  x^2 + 4 (n - m) x + m (2  m - 1))*y'[
              x]  - 2  n (2  x - 2  m + 1)* y[x] == 0,  y, 
      x]

Polyanin and Zaitsev, No. 48  (page 456)

DSolve[y'''[x]  + 
          3 a x*y''[x]  + (2  a^2 x^2 +  a + b)*y'[x]  + 
          a b x y[x] == 0,  y, x]

4. Symmetric Products of Second Order Differential Operators (new in 5.2)

A higher order ODE is a symmetric product of two or more second order ODEs if it has a solution basis which is composed of products of solutions for the second order ODEs. For instance, a fourth order ODE is a symmetric product if there exist two second order ODEs with solutions bases  {(u^1), (u^2)} and {(v^1), (v^2)} such that  {(u^1) (v^1),  (u^1) (v^2), (u^2) 
          (v^1) ,  (u^2) (v^2)} is a basis for the general solution of the fourth order ODE.

In [7], Mark van Hoeij described an algorithm for solving such a fourth order ODE. His algorithm breaks up into two cases.

Case 1. Here, the  two components in the symmetric product have missing terms in y'[x] i.e. they are second order ODEs of the type

y''[x] + a[x]*y[x] == 0.

In this case, we simply find the symmetric product of the second order ODEs with general coefficients and compare coefficients with the 
input equation

y''''x] + b[x]*y'''[x] + c[x]*y''[x] + d[x]*y'[x] + e[x]*y[x] == 0

to arrive at the result.

Case 2. In this case, the second order ODEs are of the type

y''[x] + a1[x]*y'[x] + a2[x]*y[x] ==  0

and the algorithm reduces the problem to finding rational-exponential  solutions of an associated third order linear ODE. The relevant solutions c(x) of this third order ODE have the property that ((c(x))^4) is a rational function of x. Finding these rational-exponential solutions is computationally expensive particularly if there are many singularities in the equation. If this step is successful, then a modified fourth order operator to which Case 1 applies can be obtained.

In both cases, once the second order operators have been determined, we use the available methods for second order ODEs to return a general solution for the given fourth order ODE.

We will now give some examples which work with the new implementation.

Example 1 (from van Hoeij's paper)

Timing[sol=DSolve[4 (x^4) y[x]-((8 y'[x])/x)+4 y''[x]-((2 Derivative[3][y][x])/x)+Derivative[4][y][x]==0,y,x]]

4*x^4*y[x] - (8*Derivative[1][y][x])/x + 
              4*Derivative[2][y][
                  x] - (2*Derivative[3][y][x])/
                x + Derivative[4][y][x] /. 
            sol[[1]] /. {x ->  Random[]} // Simplify // 
      Chop

Example 2 (Airy times Trigonometric)

DSolve[(((2 + 3 x + 3 (x^2) + (x^3)) y[x])/(1 + x))+(((-5 - x) y'[x])/(1 + x))+(2 - 2 x) y''[x]-((Derivative[3][y][x])/(1 + x))+Derivative[4][y][x]==0,y,x]

Simplify[(((2 + 3 x + 3 (x^2) + (x^3)) y[x])/(1 + x))+(((-5 - x) y'[x])/(1 + x))+(2 - 2 x) y''[x]-((Derivative[3][y][x])/(1 + x))+Derivative[4][y][x]/.%[[1]]]

Example 3 (Airy times Liouvillian)

DSolve[Derivative[4][y][x]+((2 (-8 (x^2) + (x^3) + 4 + x) Derivative[3][y][x])/(-8 x + 2 + (x^2)))+(((-8 (x^3) + (x^4) - 4 x + (x^2) + 4) y''[x])/(-8 x + 2 + (x^2)))-((2 ((x^2) - 2) y'[x])/(-8 x + 2 + (x^2)))-((x (-12 (x^3) + (x^4) - 24 x + 35 (x^2) + 6) y[
                    x])/(-8 x + 2 + (x^2)))==0,y,x]

Simplify[Derivative[4][y][x]+((2 (-8 (x^2) + (x^3) + 4 + x) Derivative[3][y][x])/(-8 x + 
              2 + (x^2)))+(((-8 (x^3) + (x^4) - 4 x + (x^2) + 4) y''[x])/(-8 x + 
              2 + (x^2)))-((2 ((x^2) - 2) y'[x])/(-8 x + 
              2 + (x^2)))-((x (-12 (x^3) + (x^4) - 24 x + 35 (x^2) + 6) y[
                  x])/(-8 x + 2 + (x^2)))/.%[[1]]]

Example 4 (Bessel times Exponential)

Timing[DSolve[Derivative[4][y][x]+((4 (x - 1) Derivative[3][y][x])/(x (-1 + 2 x)))-((2 (2 (x^3) + (x^2) + 1) y''[x])/((x^2) (-1 + 2 x)))-((4 ((x^2) - 1 - 2 x) y'[x])/((x^2) (-1 + 2 x)))+(((2 (x^3) - 5 (x^2) + 6 x - 1) y[
                      x])/((x^2) (-1 + 2 x)))==0,y,x]]

FullSimplify[Derivative[4][y][x]+((4 (x - 1) Derivative[3][y][x])/(x (-1 + 2 x)))-((2 (2 (x^3) + (x^2) + 1) y''[x])/((x^2) (-1 + 2 x)))-((4 ((x^2) - 1 - 2 x) y'[x])/((x^2) (-1 + 2 x)))+(((2 (x^3) - 5 (x^2) + 6 x - 1) y[
                  x])/((x^2) (-1 + 2 x)))/.%[[2]]]

Example 5 (Airy times HermiteH)

DSolve[(((4 + 5 x + 4 (x^2) + 7 (x^3) + 6 (x^4) + 3 (x^5) + 
                        (x^6)) y[x])/(1 + x + (x^2)))+(((-5 + x + 5 (x^2) + 2 (x^3)) y'[x])/(1 + x + 
                  (x^2)))+(2 - 2 x + 2 (x^2)) y''[x]+(((-1 - 2 x) Derivative[3][y][x])/(1 + x + 
                  (x^2)))+Derivative[4][y][x]==0, y,x]//Timing

(((4 + 5 x + 4 (x^2) + 7 (x^3) + 6 (x^4) + 3 (x^5) + 
                        (x^6)) y[x])/(1 + x + (x^2)))+(((-5 + x + 5 (x^2) + 2 (x^3)) y'[x])/(1 + x + 
                  (x^2)))+(2 - 2 x + 2 (x^2)) y''[x]+(((-1 - 2 x) Derivative[3][y][x])/(1 + x + 
                  (x^2)))+Derivative[4][y][x]/.%[[2]]/. {x ->  Random[]}//Simplify//Chop

5. Solution of  Higher Order ODEs using Factorization (new in 5.2)

We will begin by recalling the basic definitions and ideas underlying the factorization of differential operators and solution of the corresponding differential equations.

A linear differential operator P of order 'n' is reducible if it can be written as

P(y) = R(Q(y))

where Q is a linear differential operator of order 'm' with m <  n. Q is said to be a right factor of P while R is called a left factor of P.

Clearly, every solution of Q (i.e.of Q(y) = 0) is also a solution of P. Thus, any basis {y1, y2,..., ym} for the solution space of Q is a subbasis of the solution space of P.

Assume that a right factor Q has been found and that a basis for its solution space has been determined.Then, this can be completed to give a basis for P by solving the (n-m) non-homogeneous ODEs.

Q(y) =  z1, Q(y) =  z2, ...

where {z1, z2, ...} is a basis for the solution space of  R.

Since DSolve does not return partial solutions, an algorithm for factorization will succeed on a given example if:

1. We can find a right factor Q for it.

2. We can determine the corresponding left factor R (this is an easy step).

3. We can find a basis for the solution spaces of Q and R (not always possible).

The problem of finding first order right factors is equivalent to finding rational-exponential solutions of the ODE. This functionality has been available since Mathematica 4.1 and is illustrated by the following example.

leftfactor = (D[#, {x, 2}] - x*#) &;

rightfactor = y'[x] - 
          x*y[x];

eqn = 
      Collect[leftfactor[
          rightfactor], {y[x], y'[x], y''[x], y'''[
            x]}]

sol = DSolve[eqn == 0, y, x] // Timing

Simplify[eqn /. sol[[2]]]

For this reason, we will focus on the problem of finding right factors of order 2 or greater in this section.

We will now sketch the basic ideas underlying the Schwarz-Bronstein algorithm for factorization of differential operators.

Let P be a linear differential operator. Then the subdeterminants of the Wronskian of P satisfy certain linear ODEs called the associated equations for P. The coefficients of the associated equations are polynomial functions of  the coefficients of P.

The coefficients of P can be recovered using the rational-exponential solutions of the associated equations.

If Q is a candidate right factor of P, then the subdeterminants of its Wronskian satisfy a subset of the associated equations for P (which are known to us). Hence we can determine the coefficients of Q using the associated equations for P.

In general, it is difficult even to set up the associated equations (this requires some form of Gaussian elimination). For third order ODEs, the expressions are fairly simple and are given in the paper by Schwarz [3].  We have implemented an efficient algorithm due to Bronstein [1] which works in the general case and allows us to search for factors of different orders in succession (starting from order  two).

Once the factors of the ODE have been found, it is straightforward to find the general solution as explained in the Introduction.

Here are some examples which work with our implementation of the Schwarz-Bronstein algorithm in V 5.2.

Example 1  -  Airy functions - order 3

leftfactor = (D[#, x] - #) &;

\(\(\(rightfactor = \(y''\)[x] - x*y[x];\)\(
\)
    \)\)

eqn = 
      Collect[leftfactor[
          rightfactor], {y[x], y'[x], y''[x], y'''[
            x]}]

Timing[sol = DSolve[eqn == 0, y, x]]

Simplify[eqn /. sol[[1]]]

Example 2  -  Bessel functions  - order 3

leftfactor = (D[#, x] - x*#) &;

rightfactor = 
        x^2*y''[x] + x*y'[x] + (x^2 - 1)*y[x];

eqn = 
      Collect[Expand[
          leftfactor[rightfactor]], {y[x], y'[x], y''[x], y'''[
            x]}]

Timing[sol = DSolve[eqn == 0, y, x]]

FullSimplify[eqn /. sol[[1]]]

Example 3  -  Legendre functions  - order 3

\(\(\(leftfactor = \((D[#, x] + x*#)\) &;\)\(
\)
    \)\)

rightfactor = (1 - x^2)*y''[x] - 
          2*x*y'[x] + (3/5) ((3/5) + 1)*y[x];

\(\(\(eqn = 
      Collect[leftfactor[
          rightfactor], {y[x], \(y'\)[x], \(y''\)[x], \(y'''\)[x]}]\)\(
\)
    \)\)

Timing[sol = DSolve[eqn == 0, y, x]] // InputForm

eqn /. sol[[1]] /. {x -> Random[]} // N // 
        Simplify // Chop

Example 4  -  Bronstein's original example  - order 4

eqn  =  y''''[t]  - 2 t y''[t]  - 2 y'[t]  + 
          t^2 y[t];

(sol = DSolve[eqn == 0,  y, t]) // Timing // 
      InputForm

(y1[
            t_] = (y[t] /. sol[[1]] // FunctionExpand // 
              Simplify)) // InputForm

eqn /. y ->  y1 /.  t ->  Random[] // Simplify // 
      Chop

Example 5  - Bessel functions - order 4

leftfactor  =  (D[#,  {x, 2}] - #) &;

rightfactor  = 
        x^2*y''[x]  +  x*y'[x] + x^2*y[x];

eqn  = Together[leftfactor[rightfactor]]

(sol = DSolve[eqn == 0,  y, x]) // Timing

eqn /. sol[[1]] /. {x ->  Random[]} // Simplify // 
      Chop

Example 6 - Airy functions - order 5

leftfactor1 =  (x*D[#,  {x,  2}] + D[#, x] + 
              x*#) &;

leftfactor2 = (D[#, x] - (1/x)*#) &;

rightfactor  = (y''[x] - x*y[x]);

(eqn = Together[leftfactor1[leftfactor2[rightfactor]]]) // 
      InputForm

DSolve[rightfactor == 0,  y, x]

(sol = DSolve[eqn == 0,  y, x]) // Timing // 
      InputForm

eqn /. sol[[1]] /.  {x ->  Random[]} // Simplify // 
      Chop

6. Outlook

As we have seen in the earlier sections, the use of symmetric powers, symmetric products and factorization enlarges the scope of DSolve for solving higher ODEs in a very substantial way.

For future versions, we would like to implement other schemes for solving higher order ODEs, including ideas based on Mark van Hoeij's  Ph. D. thesis.

Another possibility is to implement new Mathematica functions for dealing with the differential operators themselves (without carrying out the final step of solving differential equations). This will open up new applications and is also worthwhile from a computational point of view because, in the examples considered by us, most of the time is spent in solving the lower order ODEs.

7. References

[1] M. Bronstein, 1994. An improved algorithm for factoring linear ordinary differential operators, Proceedings of ISSAC' 94, ACM Press, 336-340.

[2] M. Bronstein, T. Mulders and J-A. Weil, 1997. On symmetric powers of differential operators, Proceedings of ISSAC'97, ACM Press, 156-163

[3] F. Schwarz, 1989. A factorization algorithm for linear ordinary differential equations, Proceedings of ISSAC' 89, ACM Press, 17-25.

[4] M. F. Singer, 1981. Liouvillian solutions of of nth order homogeneous linear differential equations, Amer. J. Math., 103, 661-682.

[5] M. F. Singer, 1985. Solving homogeneous linear differential equations in terms of second order linear differential equations,  American Journal of Mathematics, 107, p. 663-696.

[6]. M. van Hoeij, 1996. Factorization of Differential Operators, Ph. D. thesis.

[7] M. van Hoeij, 2000. Decomposing a fourth order linear differential equation as a symmetric product, Banach Center Publications, Warsaw.

[8] M. van Hoeij. The two highest coefficients of a symmmetric power of a second order operator (preprint).