数理解析研究所講究録
第 848 巻 1993 年 100-108

100

Algorithm of solving second orderlinear ordinary differential equations
and its implementation in REDUCE
L. M. Berkovich
Dept. Algebra and Geometry, Samara State University, Russia

F. L. Berkovich
Dept. Applied Mathematics, Samara State Avia-Space University, Russia

1. Introduction

The construction of algorithms for finding formal solutions for some classes of equations
is a main purpose of any constructive theory of ordinary differential equations theory.
Explicit formulas are the most important and conclude in itself all available information.
It is necessary to have them also to convince in our physics intuition and to compare
different theories including the bounds of their applicability.
Euler, Liouville, Kummer, Jacoby and other mathematicians discovered that the basic method of integration and investigation of differential equations is an introducing theconvenient changes of variables which reduce the original equations to simpler form. However, they didn’t propose the algorithms for finding this transformations and, as a result,search for suitable substitutions had an euristic character.
From the other hand, deep relations between linear ODE and algebraic equations havebeen known for a long time and that lead to the method of factorization of differentialoperators. Here many results had a non-constructive character, too.
One of the authors developed in [1] the effective method of transformation of differential equations in which there were used both change of variables and factorization ofdifferential operators.
There exist the different approaches for receiving the Liouvillian solutions of ODE, i.e.,
those ones that may be represented in a finite form using elementary, algebraic functions
and quadratures.
An interesting though not complete, review of different methods for receiving theLiouvillian and other formal solutions of ODE may be found in [2].
Modem computer algebra systems are becoming the powerful means of PC implementation of the exact methods of ODE investigating and integrating. As a result, manyusers are able now to use practically those methods available earlier for specialists only.

101
In this paper we propose our algorithm for finding the general solutions of nonhomogeneous second order LODE
$Ly\equiv a_{2}(x)y’’+a_{1}(x)y’+a_{0}(x)y=f(x)$where coefficients
belong to some differential field $K$ and they are arbitrary
differentiable functions, possibly, containing parameters. The heart of an algorithm is asearch for variable change that reduces the correspondent homogeneous equation $Ly=0$
to one with constant coefficients.
Algorithm is implemented in computer algebra system REDUCE.
$a_{2},$ $a_{1},$$a_{0}$2. Method
2.1. Search for the Kummer-Liouville transformation
Let the equation(1)$y”+a_{1}(x)y’+a_{0}(x)y=0,$ $a_{1}\in C^{1}(I),$ $a_{0}\in C(I)$where $I=\{x|a<x<b\}$ be given. Here for simplicity we assumed $a_{2}=1$ . Let us apply
here the Kummer-Liouville $(KL)$ transformation , i.e. the variable change:$y=v(x)z(t),$ $dt=u(x)dx$ ,$u,$ $v\in C^{2}(I),$ $uv\neq 0,\forall x\in I$.

(2)

Due to the St\"ackel-Lie theorem (2) is a most general transformation that keeps the orderand linearity of an equation (1). it is a basic part of the Kummer problem (see [3]) of
reduction of(l)to the equation ofaform$z”(t)+b_{1}(t)z’(t)+b_{0}(t)z(t)=0$ ,$b_{1}\in C^{1}(J),$ $b_{0}\in C(J)$(3)where $J=\{t|c<t<d\}$ .
The Kummer problem is always solvable [4] and, therefore, there always exists theKL-transformation that reduces (1) to (3). However, the problem of reduction of (1) tothe equation with constant coefficients$z”(t)+b_{1}z(t)+b_{0}z(t)=0,$$b_{1},$$b_{0}=const$(4)is of the most essential interest.
The KL-transformation may be found usingLemma 1. The equation (1) may be reduced to (4) by the KL-transformation, for which
the kernel $u(x)$ satisfies to the second order Kumner-Schwarz equation$\frac{1}{2}\frac{u’’}{u}-\frac{3}{4}(\frac{u’}{u})^{2}-\frac{1}{4}\delta u^{2}=A_{0}(x)$

,

(5)

where
$A_{0}(x)=a_{0}- \frac{1}{4}a_{1}^{2}-\frac{1}{2}a_{1}’$

(6)

102
is a semi-invariant of (1), $\delta=b_{1}^{2}-4k$ is a discriminant of a characteristic equation$r^{2}+b_{1}r+b_{0}=0$.

(7)

Then, the multiplier $v(x)$ of KL-transformation satisfies to the equation resulted(8)$v(x)=|u|^{-1/2} \exp(-\frac{1}{2}\int a_{1}dx)\exp(\frac{1}{2}b_{1}\int udx)$

and, moreover, $v(x)$ and $u(x)$ are related by a differential equation$v”+a_{1}(x)v’+a_{0}(x)v-b_{0}u^{2}v=0$.

(9)Lemma 1 allows to find constructive advices on finding $u(x)$ .
Let us consider the equation$y_{1}’’+a_{1}y_{1}’+a_{01}y_{1}=0,$(10)$a_{01}=a_{0}(x)-b_{0}u^{2}(x)$. This equation and (1) have the same kernel
resulted from (9) after the change$u(x)$ of the KL-transformation. Since the coefficientdue to the formula$varrow y_{1}$$a_{0}(x)$(11)$a_{0}(x)=a_{01}(x)+b_{0}u^{2}(x)$as an additive part (up to the multiplicative constantcontains), then it isor $A_{0}(x)$ or theirreasonable to choose candidates for $u(x)$ among the expressions
additive parts.
Then the function $w(x)=u^{2}(x)$ must satisfy to the equation$b_{0}\neq 0$$u^{2}(x)$$a_{0}(x)$$\frac{1}{4}\frac{w’’}{w}-\frac{5}{16}(\frac{w’}{w})^{2}-\frac{1}{4}\delta w=A_{0}$.

(12)Remark: in order to apply this statements, at least one of expressions
be not zero.

$\delta$or$b_{0}$

should

2.2. Factorization
Let us say that the equation (1) admits factorization if its differential operator$L=D^{2}+a_{1}D+a_{0},$(13)$D=d/dx$

may be represented as a product of first order operators$L=(D-\alpha_{2})(D-\alpha_{1}),$$\alpha_{1}=\alpha_{1}(x),$ $\alpha_{2}=\alpha_{2}(x)$.

(14)

Here the differential analog of Viet formulas$a_{1}=-(\alpha_{1}+\alpha_{2}),$

remains solid.

$a_{0}=\alpha_{1}\alpha_{2}-a_{1}’$

(15)

103
Due to the G.Mammana’s theorem, (1) always admits fact rization.
$0$

Lemma 2. The equation (1) may be reduced to (4) by the transformation (2) and admits
factorization$L=(D- \frac{v’}{v}-\frac{u’}{u}-r_{2}u)(D-\frac{v’}{v}-r_{1}u)y=0$(16)

where$r_{1},$$r_{2}$

are roots of the characteristic equation (7).

Lemma 9. Zeroes (roots) of factorization may be represented in the form$\alpha_{1}=-\frac{1}{2}\frac{u’}{u}-\frac{1}{2}a_{1}+\frac{\sqrt{\delta}}{2}u$,$\alpha_{2}=\frac{1}{2}\frac{u’}{u}-\frac{1}{2}a_{1}-\frac{\sqrt{\delta}}{2}u$.

(17)This formula connects “zeroes” of factorization with Kummer-Liouville transformation.
Let us now find the feedback of KL-transformation with “zeroes” of factorization.
, we may require
Using the arbitrariness of determining the characteristic roots$r_{1}=r_{2}=0$ . Then abovermula takes a form
$r_{1},$ $r_{2}$$fo$$\alpha_{1}=\frac{v’}{v}$

,

(18)

$\alpha_{2}=\frac{v’}{v}+\frac{u’}{u}$from where$v=e^{\int a_{1}dx}$,$u=e^{\int(\alpha_{2}-\alpha_{1})dx}$.

(19)

So, we may formulate
Lemma 4. The equatiori (1) by transformation$y=e^{\int\alpha_{1}dx}$,$dt=e^{\int(\alpha_{2}-\alpha_{1})dx}dx$(20)

may be reduced to the equation with constant coefficients$z”(t)=0$ .

(21)Since in this case the KL-transformation is represented using the “roots” of factorization, the elementary procedure for finding factorization may be pointed.
Let us consider the Kummer-Schwarz equation (5) that due to (15) will take a form
$A_{\{}=- \frac{1}{4}(\alpha_{2}-\alpha_{1})^{2}-\frac{1}{2}(\alpha_{2}-\alpha_{1})’$(22)

from where$\alpha=\alpha_{2}-\alpha_{1}$So, we nay find$\alpha$as $\alpha=-2\sqrt{w}$ , where$w$.

is an a additive part of semi-invariant$A_{0}$.

1042.3. Fundamental system of solutions of LODE
Case 1: equation (1) admits the factorization (19).
Lemma 5. Equation (1) has the following fundamental system of solutions (FSS):$y_{1,2}=u^{-1/2} \exp(-\frac{1}{2}\int a_{1}dx)\exp(\pm\frac{\sqrt{\delta}}{2}\int udx)$$y_{1}=u^{-1/2} \exp(-\frac{1}{2}\int a_{1}dx),$,$y_{2}=y_{1} \int udx$$\delta\neq 0$,,$\delta=0$

(23)

.

(24)

Case 2: $b_{1}=b_{0}=0$ .
Lemma 6. If the factorization (14) is known, then FSS as follows:$y_{1}=e^{\int\alpha_{1}dx}$,$y_{2}=y_{1} \int e^{\int(\alpha_{2}-\alpha_{1})dx}dx$.

(25)2.4. Partial solution
Let the nonhomogeneous equation(26)$y”+a_{1}y’+a+Oy=f(x)$of correspondent equation (1) is known, then the partial
be given. If the FSS
for Lemma 6 has the form
solution$\{y_{1}, y_{2}\}$$y^{*}$(27)$y^{*}=-y_{1} \int e^{\int a_{1}dx}y_{2}fdx+y_{2}\int e^{\int a_{1}dx}y_{1}fdx$

while for Lemma 5$y^{*}= \frac{1}{2\sqrt{-b_{0}}}(y_{1}\int e^{\int a_{1}dx}y_{2}fdx-y_{2}\int e^{\int a_{1}dx}y_{1}fdx)$,$b_{0}\neq 0$.

(28)

2.5. Semi-invariants and special cases of a Kummer-Liouville transformation(invariant by the transformation of dependent vari2.5.1. Semi-invariantable).
As it was already mentioned above, has a form
$J_{0}$$J_{0}$$J_{0}=a_{0}- \frac{1}{4}a_{1}^{2}-\frac{1}{2}a_{1}’$.

$IfJ_{0}=const,$ $thentheKL$ -transformation will takea form
$y= \exp(-\frac{1}{2}\int a_{1}dx)z,$

Factorization of operator$L$$dt=dx$ .

(29)

in this case will become commutative:$L=(D+ \frac{1}{2}a_{1}+\sqrt{b_{0}})(D+\frac{1}{2}a_{1}-\sqrt{b_{0}})$

.

(30)

105
2.5.2. Semi-invariant$J_{1}$(invariant by the transformation of independent vari-able).
$J_{1}=a_{0}e^{2\int a_{1}dx}(b_{1} \int e^{-\int a_{1}dx}dx+c)^{2}$.

(31)

.

(32)

If $J_{1}=const$ , then the KL-transformation will take a form
$y=z,$

We may determine whether$dt=- \frac{e^{-\int a_{1}dx}}{b_{1}\int e^{-\int a_{1}dx}dx+c}dx$$J_{1}$

is a constant or not, using the formula$\frac{a_{1}}{\sqrt{a_{0}}}+\frac{1}{2}\frac{a_{0}’}{a_{0}\sqrt{a_{0}}}=b_{1}=const$.

(33)

2.6. Equations solvable algebraically
2.6.1. Exponential solutions
Let the equation $Ly=0$ have an exponential solution
the characteristic equation$y=e^{\lambda x}$, where $\lambda=const$ . Then(34)

$r^{2}+a_{1}(x)r+a_{0}(x)=0$

has among its roots

not a function but a number :$\lambda$$r_{1},$ $r_{2}$$r_{1,2}=- \frac{a_{1}}{2}\pm\sqrt{\frac{a_{1}^{2}}{4}-a_{0}}$Factorization

$L$.

(35)

takes a form:$L=(D+a_{1}+\lambda)(D-\lambda)$

.

(36)

2.6.2. Adjoint equations
By definition, the adjoint for $Ly=0$ equation $L^{*}y=0$ is an equation$L^{*}y\equiv y’’-a_{1}y’+(a_{0}-a_{1}’)y=0$

It admits factorization$L^{*}=(D+\alpha_{1})(D+\alpha_{2})$

.

.

(37)
(38)

Let us consider the characteristic equation for adjoint one:$r^{2}-a_{1}r+a_{0}-a_{1}’=0$

If one of its roots is a number , then the factorization$\lambda$.

(39)$L^{*}$$L^{*}=(D-a_{1}+\lambda)(D-\lambda)$

Simultaneously, the factorization$L$.

takes a form:(40)

is of the form

$L=(D+\lambda)(D+a_{1}-\lambda)$

.

(41)

106
2.6.3. Exact equation, then we have an exact equation, for which the factorization
If$\lambda=0$$L=D(D+a_{1})$$L$

has a form.

(42)

3. Implementation

The REDUCE program SOLDE which implements the algorithm described above, has
the following main characteristics:
INPUT: a2 al $a0,f$
OUTPUT: u.v,
bl $bO$ ,$/{}_{0}Coefficients$$\iota/.variable$of given equation

change$’/*coefficients$ of reduced equationof factorization

alfal alfa2,yl, $y2$ ,$|/.coefficients$

yp$’/.partial$ solution for nonhomogeneous ODE$’/.FSS$

The method is, really, an algorithm itself, so the computer program SOLDE written incomputer algebra system REDUCE, follows to the second part of this paper in detail, so
there is no need to describe it here. It should be pointed only, that for make the program
more convenient for users and faster the second author decided to include the coefficientunder consideration. It may be explained by the fact that sometimes, that was proved
theoretically, the candidate for $u(x)$ may be taken from $u(x)=1/q(x)$ , where $q(x)$ is aor its multiplicative part.
Program SOLDE was tested on hundreds of equations and it had success in 75% cases.
Failure in the left 25% are caused by two reasons: 1) the other important transformation,of Euler-Imshenetsky-Darboux, was not included to the program yet, though there alreadyexist all the necessary formulas for its application here, and 2) the algorithm doesn’t workwhen $b_{1}=b_{0}=0$ .
For illustrate the current possibilities of this program, let us name the equations fromthe demonstration file:$a_{2}$$a_{2}(x)$$***$I.Equations with constant coefficients $***$$Y’ l+A*Y’+B*Y=F(X)$$***$

II. Euler’$Y’ l+A*Y)+A*A/4*Y=F(X)$$Y^{l}$$s$

equations $***$$+A/X*Y’+B/(X*X)*Y=F(X)/(X*X)$$***$
III.Equations of a form :$Y’)+A1(x)*Y^{l}=F(X)$$***$ IY.Equations with exponential coefficients $***$$Y’+A*Y)+B*E**(2*A*X)*Y=F(X)$$***Y.Equations$ with trigonometrical functions $***$$Y’$ $‘+2*A*C0T(A*X)*Y)+(B*B-A*A)*Y=0$’$‘+(M*M+A/SIN(2*M*X)**2)*Y=0$$***YI.Equations$ with hyperbolic functions $***$***$Y$

107$Y$‘$’+2*TANH(X)*Y$ ’ $+B*Y=0$$‘+(-M*M+A/(SINH(M*X)**4))*Y=0$$***YII.Equations$ with algebraic coefficients $***$Y)$Y^{\prime l}+8*A/(5*(A*X+B))*Y^{l}+C*(A*X+B)**(1/5)/(5*(A*X+B))*Y=0$$***YIII.Equations$ with mixed coefficients $***$$Y”+2*A/X*Y’+((B*B*E**(2*C*X)-1/4)*C*C+A*(A-1)/(X*X))*Y=0$$Y))+(1/(4*X*X)+1/(X*X)/(P*LOG(X)+Q)**4))*Y=0$$***$ IX.Equations with arbitrary functions $***$$Y^{l}$$+2*F(X)*Y^{l}+(F(X)*F(X)+F(X)+G)(X)/2/G(X)-3/4*G^{l}(X)*G(X)/$$G(X)*G(X)-A*G(X)*G(X))*Y=0$***X.Equations with rational coefficients $***$$Y^{ll}+D/(A*X*X+B*X+C)**2*Y=0$$Y))+(-M*(M+1)/(X*X)+(1/(-P/(2*M+1)*X**(-M)+Q*X**(M+1))**4)*Y=0$

High efficiency of this implementation is caused by the fact that entire package wasdeveloped around the procedure SOLDE which allows to the user to investigate his ODE
from different sides, and even the program doesn’t give an answer itself, it often helps
the user in finding the approaches to his equation making more effective and faster hisoperations and testing his propositions. Thus, each ODE may be investigated during thewhole session by the procedures shown in Table 1.
Table 1.

Procedures in the package.

Example of this program in action is given below:

108
2:

case of this equation with $a=1$ was solved by Kovacic [5]another way$|/_{0}Specia1$$|/.in$solde $(x^{\sim}2,0,-a*x+3/16 , 0)$ ;$*Summary$ of the $operations*$$***************************************************************$$*The$

equation was:$*The$

semi-invariant by dependent variable:$(X^{2})*y^{l}$$‘+$(0) $*y)+$$( - ——-)*y\underline{1}\underline{6}*\underline{A}*\underline{X}-\underline{3}=0$16$——-\underline{1}\underline{6}*\underline{A}*\underline{X}-\underline{3}_{-}$$J0=$$16*X^{2}$

1/4$*The$

transformation:$y=$(X$z’$$(t)+(0)*z^{J}(t)+$leads to$*The$ factorization:)$*$) $*z$ ,dt

1$=$(———) dx
SQRT (X)$( -A)*z(t)=0$$4*SQRT(A)*X+$ SQRT(X)$*$$(D-($$L=(D-(————————))$$4*SQRT(X)*X$$4*SQRT(A)*X+$ SQRT(X)$———————–))$$4*SQRT(X)*X$$*Fundamental$
1/4$Y1=$Xsystem of solutions of $Ly=0$ :$2*SQRT(X)*SQRT(A)$$*E$

1/4$x$$Y2=——————–$$2*SQRT$ (X)*SQRT (A)$E$$**************************************************************$The previous version of the program is described in [6].
References[1] Berkovich L.M., Factorization and Tmnsformations of Ordinary Differential Equations, Saratov Univ. Publishers, Saratov, 1989 (in Russian)[2] Singer M. F., Formal Solutions of Differential Equations, J. Symb. Comp., 1990, 10, 59-94[3] Boruvka 0., Lineare Differential Tmnsformationen 2. Ordnung, Berlin, 1967

[4] Berkovich L. M., Transformations of Linear Ordinary Differential Equations, Kuibyshev
State Univ, 1978 (in Russian)[5] Kovacic J., An Algorithm for Solving Second Order Linear Homogeneous Differential equations, J. Symb. Comp., 1986, 2, 1
[6] Berkovich L. M., Gerdt V. P., Kostova Z. T., Nechaevsky M. L., Second Order Reducible
Linear Differential Equations, Preprint JINR E5-89-141, Dubna, 1989.