Constructive factorization of LPDO in two variables

Richard Beals^{+} , Elena Kartashova^{∗}

+ Yale University, New Haven, CT, USA ∗ RISC, J.Kepler University, Linz, Austria

e-mails: richard.beals@yale.edu, lena@risc.uni-linz.ac.at

Abstract

We study conditions under which a partial differential operator of arbitrary order n in two variables or ordinary linear differential operator admits a factorization with a first-order factor on the left. The factorization process consists of solving, recursively, systems of linear equations, subject to certain differential compatibility conditions. In the generic case of partial differential operators one does not have to solve a differential equation. In special degenerate cases, such as ordinary differential, the problem is finally reduced to the solution of some Riccati equation(s). The conditions of factorization are given explicitly for second- and, and an outline is given for the higher-order case.

1 Introduction

Factorization of linear differential operators (LDO) is a very well-studied problem and a lot of pure existence theorems has been proved (e.g. [1], [8], [9]). The case of linear ordinary differential operators (LODO) has been studied in more detail; various algorithms are known for factoring LODOs over different differential fields, e.g. [3], [10], [11] and others. Probably the first algorithm for the case of the simplest possible differential field (rational functions) is described in [2]. In this case factorization into irreducible factors is known to be unique [6], [7] and in general, the problem finally reduces to the solution of a Riccati equation.

Much less is known about the factorization of linear partial differential operators (LPDOs) which would appear to be much more complicated than LODO. In several variables the naive definition of factorization as representing a given n-th order operator as a composition of lower order operators does not have a nice uniqueness property; this is illustrated by an example of E. Landau (formulated for left composition) [5]:

(∂_{x} + 1)(∂_{x} + 1)(∂_{x} + x∂_{y}) = [∂_{x}^{2} + x∂_{x}∂_{y} + ∂_{x} + (2 + x)∂_{y}](∂_{x} + 1).

Nevertheless, some interesting results have been obtained also in this direction; recently Grigoriev and Schwarz [4] have given an algorithm for factoring an LPDO with separable symbol. Some generalized definitions of factorization gave rise to another factorization algorithms for factoring LPDO [12] and for factoring of systems of LPDEs with rational function coefficients [13], [14].

It seems to be a general opinion that the naive factorization – even for secondand third order differential operators in two variables – requires “the solution of a partial Riccati equation, which in turn requires to solve a general first-order ODE and possibly an ordinary Riccati equation. The bottleneck for designing a factorization algorithm for LPDO is the general first-order ODE which makes the full problem intractable at present because in general there are no solution algorithms available” [4].

We consider it possible to construct an explicit algorithm which could be used for absolute factorization of LPDO of arbitrary order n in two variables into linear factors. The word ”absolute” means that we do not fix the coefficient field from the very beginning and that our only demand on the coefficients is that they be smooth, i.e. they belong to an appropriate differential field. The procedure proposed here is to find a first order left factor (when possible) in contrast to the use of right factorization, which is common in the papers of last few decades. Of course the existence of a certain right factor of a LPDO is equivalent to the existence of a corresponding left factor of the transpose of that operator, so in principle nothing is lost by considering left factorization. Moreover taking transposes is trivial algebraically, so there is also nothing lost from the point of view of algorithmic computation.

The use of left factorization makes all the algebraic calculations needed much easier. It also yields explicit necessary and sufficient conditions on the coefficients of the operator for the existence of a factor. The simple form of these conditions allows us one to see the principal difference between factorization of LPDOs and LODOs - in the first case the factorization problem can be solved by pure algebraic methods in the ”generic” case; the LODO is one example of the various degenerate cases in which it is necessary to solve a Riccati equation in order to factor in two variables.

In Sec. 2 the general idea of factorization (in the generic case) by purely algebraic methods is presented and explicit conditions of factorization for an LPDO of order 2 are written out. It is also shown how the possibility of pure algebraic factorization for LPDO does not contradict the necessity to solve a Riccati equation in order to factor an LOPO. Similar results for LPDO and LODO of order three are demonstrated in Sec. 3, while in Sec. 4 the general procedure for factorization for LPDO of arbitrary order n is presented. Some interesting examples are constructed in Sec. 5 and Sec. 6 contains a brief concluding discussion.

2 LPDO of order 2 in two variables

We consider an operator

A_{2} = ∑_{j+k≤2} a_{jk}∂_{x}^{j}∂_{y}^{k} = a_{20}∂_{x}^{2} + a_{11}∂_{x}∂_{y} + a_{02}∂_{y}^{2} + a_{10}∂_{x} + a_{01}∂_{y} + a_{00}.  (1)

with smooth coefficients and look for a factorization

A_{2} = (p_{1}∂_{x} + p_{2}∂_{y} + p_{3})(p_{4}∂_{x} + p_{5}∂_{y} + p_{6}).

Let us write down the equations on p_{i} explicitly, keeping in mind the rule of left composition, i.e. that ∂_{x}(α∂_{y}) = ∂_{x}(α)∂_{y} + α∂_{xy}.

Then in all cases

\begin{cases}a_{20} & = p_{1}p_{4} \\ a_{11} & = p_{2}p_{4} + p_{1}p_{5} \\ a_{02} & = p_{2}p_{5} \\ a_{10} & = L(p_{4}) + p_{3}p_{4} + p_{1}p_{6} \\ a_{01} & = L(p_{5}) + p_{3}p_{5} + p_{2}p_{6} \\ a_{00} & = L(p_{6}) + p_{3}p_{6} \end{cases}  (2SysP)

where we use the notation L = p_{1}∂_{x} + p_{2}∂_{y}.

We begin with some preliminary remarks. To some extent the factorization problem is to be treated as a local problem; for example if some coefficient does not vanish identically, we may wish to restrict to a region where it has no zeros. Of course insofar as the formulas obtained are analytic functions of the data, if the data are also analytic then the formulas may continue globally. With this in mind, we note that if the operator is not (globally) of order one, then we may, after a linear change of variables if necessary, assume that

a_{20} ≠ 0.

Then necessarily p_{1} ≠ 0, and we may assume without loss of generality that

p_{1} = 1.

Then the first three equations of 2SysP, describing the highest order terms are equations in the variables p_{2}, p_{4}, p_{5} and to find them we have to find roots of a quadratic polynomial. Having computed p_{2}, p_{4}, p_{5} one can plug them into two next equations of 2SysP and get a linear system of equations in two variables p_{3}, p_{6} which can easily be solved if ω is a simple root. Then the last equation of 2SysP will give us the condition of factorization.

In fact, at the first step from

\begin{cases}a_{20} = p_{4} \\ a_{11} = p_{2}p_{4} + p_{5} \\ a_{02} = p_{2}p_{5} \end{cases}

it follows that

P_{2}(−p_{2}) := a_{20}(−p_{2})^{2} + a_{11}(−p_{2}) + a_{02} = 0

The choice of different roots of the characteristic polynomial P_{2} gives us different (potential) factorizations of the initial operator. (As we shall see, there may be no such factorization, or there may be one with one choice of root and not with another choice.) Let ω be a root of P_{2}, and take

p_{2} = −ω, L = ∂_{x} − ω∂_{y}.

This leads to a linear system for p_{4}, p_{5} with ω as parameter:

\begin{cases}p_{4} = a_{20}, \\ −ωp_{4} + p_{5} = a_{11} \end{cases}

or, in matrix form,

[\begin{matrix}1 & 0 \\ −ω & 1 \end{matrix}] [\begin{matrix}p_{4} \\ p_{5} \end{matrix}] = [\begin{matrix}a_{20} \\ a_{11} \end{matrix}]; [\begin{matrix}p_{4} \\ p_{5} \end{matrix}] = [\begin{matrix}1 & 0 \\ ω & 1 \end{matrix}] [\begin{matrix}a_{20} \\ a_{11} \end{matrix}] .

Thus

\begin{cases}p_{1} = 1 \\ p_{2} = −ω \\ p_{4} = a_{20} \\ p_{5} = a_{20}ω + a_{11} \end{cases}  (2Pol)

At the second step, from

\begin{cases}a_{10} & = L(p_{4}) + p_{3}p_{4} + p_{1}p_{6} \\ a_{01} & = L(p_{5}) + p_{3}p_{5} + p_{2}p_{6} \end{cases}

and from (2Pol) we get

\begin{cases}a_{10} & = La_{20} + p_{3}a_{20} + p_{6} \\ a_{01} & = L(a_{11} + a_{20}ω) + p_{3}(a_{11} + a_{20}ω) − ωp_{6}. \end{cases}  (2Lin)

If the root ω is simple, i.e. P′_{2}(ω) = 2a_{20}ω + a_{11} ≠ 0, then these equations have the unique solution

p_{3} = \frac{ωa_{10} + a_{01} − ωLa_{20} − L(a_{20}ω + a_{11})}{2a_{20}ω + a_{11}}

p_{6} = \frac{(a_{20}ω + a_{11})(a_{10} − L)a_{20} − a_{20}(a_{01} − L(a_{20}ω + a_{11})}{2a_{20}ω + a_{11}}.

At this point all coefficients p_{1}, p_{2}, ..., p_{6} have been computed.

At the third step from

a_{00} = L(p_{6}) + p_{3}p_{6}

the corresponding condition of factorization can be written out explicitly:

a_{00} = L \{\frac{ωa_{10} + a_{01} − L(2a_{20}ω + a_{11})}{2a_{20}ω + a_{11}}\} + \frac{ωa_{10} + a_{01} − L(2a_{20}ω + a_{11})}{2a_{20}ω + a_{11}} ×

×\frac{a_{20}(a_{01} − L)(a_{20}ω + a_{11}) + (a_{20}ω + a_{11})(a_{10} − La_{20})}{2a_{20}ω + a_{11}},

and we have just to check it. Thus factorization in two variables, when possible, can (generally) be accomplished by purely algebraic means.

Remark. In the preceding derivation we assumed that ω is a simple root so that (a_{11} + 2a_{20}ω) ≠ 0. Suppose now that the condition fails identically in some region. Then by eliminating p_{6} from (2Lin) we obtain the necessary condition

a_{10}ω + a_{01} − L(a_{20}ω + a_{11}) − ωLa_{20} = 0.  (2Degen)

If this condition is satisfied, then we may use the first equation of (2Lin) to determine p_{6} from p_{3}. Then the last equation of (2SysP) gives a Riccati equation for p_{3}.

Let us look explicitly at the case of an ordinary differential operator of second order as a particular case of our initial LPDO. Here ω is a double root of the characteristic equation and the necessary condtion (2Degen) is necessarily satisfied. In fact 2SysP takes the form

\begin{cases}a_{20} & = p_{4} \\ a_{10} & = L(p_{4}) + p_{3}p_{4} + p_{6} \\ a_{00} & = L(p_{6}) + p_{3}p_{6} \end{cases}  (2SysO)

i.e. we have only three non-trivial equations. The first equation shows again that p_{4} = a_{20}, Lp_{4} = ȧ_{20} is known, and p_{6} can written out explicitly as p_{6} = (a_{10} − ȧ_{20} − a_{20})p_{3}. Substitution of this expression for p_{4} into the equation for a_{00} gives us immediately

a_{00} = ȧ_{10} − ä_{20} − a_{20}ṗ_{3} − ȧ_{20}p_{3} + (a_{10} − ȧ_{20} − a_{20}p_{3})p_{3}  (2)

which is the Riccati equation

ψ̇ + ψ^{2} + \frac{2ȧ_{20} − a_{10}}{a_{20}} ψ + \frac{a_{00} + ä_{20} − ȧ_{10}}{a_{20}} = 0, for p_{3} = ψ.

In order to be able to find the general solution explicitly one has to know at least one of its solutions. Of course the solution will depend on the explicit form of a_{ij} = a_{ij}(x).

3 LPDO of order 3

Now we consider an operator

A_{3} = ∑_{j+k≤3} a_{jk}∂_{x}^{j}∂_{y}^{k} = a_{30}∂_{x}^{3} + a_{21}∂_{x}^{2}∂_{y} + a_{12}∂_{x}∂_{y}^{2} + a_{03}∂y^{3}  (3)
+a_{20}∂_{x}^{2} + a_{11}∂_{x}∂_{y} + a_{02}∂_{y}^{2} + a_{10}∂_{x} + a_{01}∂_{y} + a_{00}.

with smooth coefficients and look for a factorization

A_{3} = (p_{1}∂_{x} + p_{2}∂_{y} + p_{3})(p_{4}∂_{x}^{2} + p_{5}∂_{x}∂_{y} + p_{6}∂_{y}^{2} + p_{7}∂_{x} + p_{8}∂_{y} + p_{9}).

The conditions of factorization are described by the following system:

\begin{cases}a_{30} & = p_{1}p_{4} \\ a_{21} & = p_{2}p_{4} + p_{1}p_{5} \\ a_{12} & = p_{2}p_{5} + p_{1}p_{6} \\ a_{03} & = p_{2}p_{6} \\ a_{20} & = L(p_{4}) + p_{3}p_{4} + p_{1}p_{7} \\ a_{11} & = L(p_{5}) + p_{3}p_{5} + p_{2}p_{7} + p_{1}p_{8} \\ a_{02} & = L(p_{6}) + p_{3}p_{6} + p_{2}p_{8} \\ a_{10} & = L(p_{7}) + p_{3}p_{7} + p_{1}p_{9} \\ a_{01} & = L(p_{8}) + p_{3}p_{8} + p_{2}p_{9} \\ a_{00} & = L(p_{9}) + p_{3}p_{9} \end{cases}  (3SysP)

with L = p_{1}∂_{x} + p_{2}∂_{y}.

Once again we may assume without loss of generality that the coefficient of the term of highest order in ∂_{x} does not vanish, and that the linear factor is normalized:

a_{30} ≠ 0, p_{1} = 1.

The first four equations of 3SysVar describing the highest order terms are equations in the variables p_{2}, p_{4}, p_{5}, p_{6}. Solving these equations requires the choice of a root −p_{2} of a certain polynomial of third degree. Once this choice has been made, the remaining top order coefficients p_{4}, p_{5}, p_{6} are easily found. The top order coefficients can now be plugged into the next four equations of 3SysP. The first three of these four equations will now be a linear system of equations in the variables p_{3}, p_{7}, p_{8} which is easily solved. The next equation is now a linear equation on variable p_{9} which means that all variables p_{i}, i = 1, ..., 9 have been found. The last two equations of 3SysP will give us then the conditions of factorization.

Namely, at the first step from

\begin{cases}a_{30} = p_{4} \\ a_{21} = p_{2}p_{4} + p_{5} \\ a_{12} = p_{2}p_{5} + p_{6} \\ a_{03} = p_{2}p_{6} \end{cases}

it follows that

P_{3}(−p_{2}) := a_{30}(−p_{2})^{3} + a_{21}(−p_{2})^{2} + a_{12}(−p_{2}) + a_{03} = 0.

As for the case of second order, taking p_{2} = −ω, where ω is a root of the characteristic polynomial P_{3} we get a linear system in p_{4}, p_{5}, p_{6} with ω as parameter. Then again

p_{2} = −ω, L = ∂_{x} − ω∂_{y},

which leads to

\begin{cases}a_{30} = p_{4} \\ a_{21} = −ωp_{4} + p_{5} \\ a_{12} = −ωp_{5} + p_{6} \end{cases}

i.e.

[\begin{matrix}1 & 0 & 0 \\ −ω & 1 & 0 \\ 0 & −ω & 1 \end{matrix}] [\begin{matrix}p_{4} \\ p_{5} \\ p_{6} \end{matrix}] = [\begin{matrix}a_{30} \\ a_{21} \\ a_{12} \end{matrix}];  [\begin{matrix}p_{4} \\ p_{5} \\ p_{6} \end{matrix}] = [\begin{matrix}1 & 0 & 0 \\ ω & 1 & 0 \\ ω^{2} & ω & 1 \end{matrix}] [\begin{matrix}a_{30} \\ a_{21} \\ a_{12} \end{matrix}] .

Thus

\begin{cases}p_{1} = 1 \\ p_{2} = −ω \\ p_{4} = a_{30} \\ p_{5} = a_{30}ω + a_{21} \\ p_{6} = a_{30}ω^{2} + a_{21}ω + a_{12}. \end{cases}  (3Pol)

At the second step, from

\begin{cases}a_{20} & = L(p_{4}) + p_{3}p_{4} + p_{1}p_{7} \\ a_{11} & = L(p_{5}) + p_{3}p_{5} + p_{2}p_{7} + p_{1}p_{8} \\ a_{02} & = L(p_{6}) + p_{3}p_{6} + p_{2}p_{8} \end{cases}

and (2Pol) we get

\begin{cases}a_{20} − La_{30} & = p_{3}a_{30} + p_{7} \\ a_{11} − L(a_{30}ω + a_{21}) & = p_{3}(a_{30}ω + a_{21}) − ωp_{7} + p_{8} \\ a_{02} − L(a_{30}ω^{2} + a_{21}ω + a_{12}) & = p_{3}(a_{30}ω^{2} + a_{21}ω + a_{12}) − ωp_{8}. \end{cases}  (3Lin∗)

As a linear system for p_{3}, p_{7}, p_{8} this has determinant

3a_{30}ω^{2} + 2a_{21}ω + a_{12} = P′(ω),

so if ω is a simple root the system has unique solution

p_{3} = \frac{ω^{2}(a_{20} − La_{30}) + ω(a_{11} − L(a_{30}ω + a_{21} + a_{02})) − L(a_{30}ω^{2} + a_{21}ω + a_{12})}{3a_{30}ω^{2} + 2a_{21}ω + a_{12}};

p_{7} = \frac{a_{20} − La_{30}}{3a_{30}ω^{2} + 2a_{21}ω + a_{12}} − \frac{a_{30}}{3a_{30}ω^{2} + 2a_{21}ω + a_{12}} · p_{3};

p_{8} = \frac{ω(a_{20} − La_{30}) + a_{11} − L(a_{30}ω + a_{21})}{3a_{30}ω^{2} + 2a_{21}ω + a_{12}} − \frac{a_{30}ω + a_{21}}{3a_{30}ω^{2} + 2a_{21}ω + a_{12}} · p_{3}.

In order to find the last coefficient p_{9} we use the next equation of (3SysVar), namely:

a_{10} = L(p_{7}) + p_{3}p_{7} + p_{1}p_{9},  (3Lin ∗ ∗).

At this point all coefficients p_{i}, i = 1, ...9 have been computed, under the assumption that ω is a simple root.

At the third step from

\begin{cases}a_{01} = L(p_{8}) + p_{3}p_{8} + p_{2}p_{9} \\ a_{00} = L(p_{9}) + p_{3}p_{9} \end{cases}

all the necessary conditions for factorization can be written out. We do not do so here because the formulas are tedious and do not add anything to understanding the main idea. If the conditions are satisfied, the explicit factorization formulae could be written out as for the second-order operator. The difference is that in this case the polynomial defined by the highest order terms is of degree 3 and we have not one but two conditions of factorization.

Regarding now LODO as a particular case of LPDO let rewrite 3SysP with p_{1} = 1 as

\begin{cases}a_{30} & = p_{4} \\ a_{20} & = p_{3}p_{4} + p_{1}p_{7} + ṗ_{4} \\ a_{10} & = p_{3}p_{7} + p_{9} + ṗ_{7} \\ a_{00} & = p_{3}p_{9} + ṗ_{9}. \end{cases}  (3SysO)

As in the second order case ω = 0 is a multiple root of P_{3}. After eliminating p_{4} = a_{30} and solving the second and third equations for p_{7} and p_{9}, respectively, in terms of p_{3}, we see finally that the factorization problem for a third-order LODO is equivalent to the following system of first order nonlinear ODEs for the single unknown function p_{3}:

\begin{cases}a_{10} = p_{7}(a_{20} − ȧ_{30} − p_{7})/a_{30} + p_{9} + ṗ_{7} \\ a_{00} = p_{9}(a_{20} − ȧ_{30} − p_{7})/a_{30} + ṗ_{9}. \end{cases}

or, equivalently, a certain second-order nonlinear ODE.

4 LPDO of arbitrary order

Now consider the factorization problem for the LPDO of order n in the general form

A_{n} = ∑_{j+k≤n} a_{jk}∂_{x}^{j}∂_{y}^{k} = (p_{1}∂_{x} + p_{2}∂_{y} + p_{3})( ∑_{j+k<n} p_{jk}∂_{x}^{k}∂_{y}^{j})  (4)

and use the same assumptions and notations as above:

• some a_{jk} with j + k = n is not = 0, so

• after a linear change of variables one can assume that a_{n0} ≠ 0,

• set p_{1} = 1 and p_{jk} = 0 if j < 0 or k < 0,

• denote p_{2} as −ω and L = p_{1}∂_{x} + p_{2}∂_{y} = ∂_{x} − ω∂_{y}. Then the equations describing terms of highest order take the form p_{1}p_{j−1,k}+p_{2}p_{j,k−1} = a_{jk}, i.e.

p_{j−1,k} − ωp_{j,k−1} = a_{jk}, j + k = n,  (5)

while the equations of lower order take the form

Lp_{jk} + p_{3}p_{jk} + p_{j−1,k} − ωp_{j,k−1} = a_{jk}, j + k < n  (6)

Keeping in mind that j = n − k for the first n equations (5) let us rewrite them as

p_{n−1−k,k} − ωp_{n−k,k−1} = a_{n−k,k}, k = 0, 1, ..., n

or in matrix form as

[\begin{matrix}1 & 0 & 0 & ... & 0 \\ −ω & 1 & 0 & ... & 0 \\ 0 & −ω & 1 & ... & 0 \\ ... \\ 0 & 0 & ... & −ω & 1 \end{matrix}] [\begin{matrix}p_{n−1,0} \\ p_{n−2,1} \\ p_{n−3,2} \\ ... \\ p_{0,n−1} \end{matrix}] = [\begin{matrix}a_{n,0} \\ a_{n−1,1} \\ a_{n−2,2} \\ ... \\ a_{1,n−1} \end{matrix}] .  (7)

If we knew ω, the unique solution could be written out as

[\begin{matrix}p_{n−1,0} \\ p_{n−2,1} \\ p_{n−3,2} \\ ... \\ p_{0,n−1} \end{matrix}] = [\begin{matrix}1 & 0 & 0 & ... & 0 \\ ω & 1 & 0 & ... & 0 \\ ω^{2} & ω & 1 & ... & 0 \\ ... \\ ω^{n−1} & ω^{n−2} & ... & ω & 1 \end{matrix}] [\begin{matrix}a_{n,0} \\ a_{n−1,1} \\ a_{n−2,2} \\ ... \\ a_{1,n−1} \end{matrix}] ,

so

p_{n−k,k} = a_{n,0}ω^{k} + a_{n−1,1}ω^{k−1} + ⋯ + a_{n−k,k}.  (8)

The equation with k = n, together with ωp_{0,n−1} = a_{0,n} gives us finally

P_{n}(ω) := a_{n,0}ω^{n} + a_{n−1,1}ω^{n−1} + ... + a_{1,n−1}ω + a_{0,n} = 0.

Choosing ω to be any root of this polynomial, we get then a unique determination of p_{j,k} for j + k = n.

Consider now the n equations of (6) corresponding to the terms of order n − 1, i.e. the case j + k = n − 1:

p_{n−2−k,k} − ωp_{n−k−1,k−1} = b_{n−1−k,k} − p_{3}p_{n−1−k,k},  (9)

where the functions p_{n−1−k,k} and

b_{n−1−k,k} = a_{n−1−k,k} − Lp_{n−1−k,k}

are known already. In matrix form these equations are

[\begin{matrix}p_{n−1,0} & 1 & 0 & 0 & ... & 0 \\ p_{n−2,1} & −ω & 1 & 0 & ... & 0 \\ p_{n−3,2} & 0 & −ω & 1 & ... & 0 \\ ... \\ p_{0,n−1} & 0 & 0 & ... & 0 & −ω \end{matrix}] [\begin{matrix}p_{3} \\ p_{n−2,0} \\ p_{n−3,1} \\ ... \\ p_{0,n−2} \end{matrix}] = [\begin{matrix}b_{n−1,0} \\ b_{n−2,1} \\ b_{n−3,2} \\ ... \\ b_{0,n−1} \end{matrix}]

Inductively, the determinant of a matrix of this form is

(−1)^{n−1}[p_{n−1,0}ω^{n−1} + p_{n−2,1}ω^{n−2} + ⋯ + p_{0,n−1}].

Equations (8) imply that the term in brackets is

na_{n,0}ω^{n−1} + (n − 1)a_{n−2,1}ω^{n−2} + ⋯ + a_{1,n−1} = P′_{n}(ω).

Therefore if ω is a simple root, the equations (9) for p_{3} and the p_{n−2−k,k} have a unique solution. Multiplying the j-th equation of the system by ω^{n−j} and adding, we determine

P′_{n}(ω) p_{3} = b_{n−1,0}ω^{n−1} + b_{n−2,1}ω^{n−2} + ⋯ + b_{0,n−1}.  (10)

If ω is a simple root, then this equation determines p_{3}, and we may determine the p_{ni2−k,k} by rewriting the first n − 1 of these equations with the p_{3} terms on the right hand side, so that the solution for the p_{n−2−k,k} in terms of p_{3} and other known functions is

[\begin{matrix}p_{n−2,0} \\ p_{n−3,1} \\ p_{n−4,2} \\ ... \\ p_{0,n−2} \end{matrix}] = [\begin{matrix}1 & 0 & 0 & ... & 0 \\ ω & 1 & 0 & ... & 0 \\ ω^{2} & ω & 1 & ... & 0 \\ ... \\ ω^{n−2} & ω^{n−3} & ... & ω & 1 \end{matrix}] [\begin{matrix}b_{n−1,0} − p_{3}p_{n−1,0} \\ b_{n−2,1} − p_{3}p_{n−2,1} \\ b_{n−3,2} − p_{3}p_{n−3,2} \\ ... \\ b_{1,n−2} − p_{3}p_{1,n−2} \end{matrix}]  (11)

as before.

Continuing to assume that ω is a simple root, let us follow this process one more step. The first n − 2 of the n − 1 equations for the coefficients p_{n−3−k,k} can be written in a matrix form like (7), with known functions on the right hand side, so the solution can be written down as above. Thus at this point we have determined the coefficients p_{j} of the first factor and all coefficients p_{jk}, j +k ≤ n−3 of the second factor. However there is one more equation at this step:

−ωp_{0,n−3} = a_{1,n−3} − (L + p_{3})p_{1,n−3}.

Since both sides have already been determined, this is a necessary condition for the existence of a factorization with this choice of root ω. Continuing, the same situation occurs at each step. Thus if the condition is satisfied at one step we may proceed to the next step, and either the process fails or there are exactly n − 1 polynomial conditions that are satisfied by the coefficients, the (simple) root ω, and their derivatives.

Suppose we drop the assumption that ω is a simple root and assume that P′_{n}(ω) ≡ 0 in some region. Then the equation (10), with zero right hand side, is the necessary and sufficient condition for solvability of the full system of equations for p_{3} and the p_{n−2−k.k}. If this condition is satisfied, then we may make any choice of p_{3} and determine the p_{n−2−k,k} from (11). At subsequent stages we will encounter factorization conditions that are expressible as polynomial equations in p_{3} and its derivatives (with coefficients that are polynomials in the a_{jk} and their derivatives and in the root ω and its derivatives); these can be thought of as generalized Riccati equations for p_{3}.

In summary: starting with a simple root of the polynomial P_{n}(ω) in the generic case a_{n,0} ≠ 0, the factorization process proceeds algebraically (in the coefficients a_{jk}, the root ω, and their derivatives) as far as it can, as determined by n − 2 polynomial obstructions. Starting with a double root, then subject to one such polynomial constraint, the remaining constraints are n − 2 generalized Riccati equations for p_{3}.

We call this eliminating procedure order-reduction for obvious reasons. In trying to factor an operator of order n as a product of a first order operator and an operator of order n − 1, we look at the equations associated to the terms of given degree of the operator of order n, in descending order. The first set of equations gives us the top order terms of both factors (modulo choice of a root of the characteristic equation), the next set determines the next lower order coefficients, and so on.

5 Some examples

5.1 General second-order LPDO

• Example 5.1.1. Hyperbolic operator:

a_{20} = 0, a_{11} = 1, a_{02} = 0, i.e. L_{2} := ∂_{x}∂_{y} + a_{10}∂_{x} + a_{01}∂_{y} + a_{00}.

In this case either of two different factorizations may be possible:

– (5.1.1.1 ) if a_{00} = ∂_{x}a_{10} + a_{10}a_{01}, then

L_{2} = (∂_{x} + a_{01}) (∂_{y} + a_{10})

or – (5.1.1.2 ) if a_{00} = ∂_{y}a_{01} + a_{10}a_{01}, then

L_{2} = (∂_{y} + a_{10}) (∂_{x} + a_{01}) .

Consider the case of such an operator known as a Poisson operator:

L_{2} = ∂_{x}∂_{y} + \frac{1}{x + y} (α∂_{x} + β∂_{y}) + \frac{γ}{(x + y)^{2}}

where α, β, γ are constants.

Due to (5.1.1.1) and (5.1.1.2) factorization is only possible if γ = α(β − 1) or γ = β(α − 1) and its explicit forms are, respectively:

L_{2} = (∂_{x} + \frac{β}{x + y} ) (∂_{y} + \frac{α}{x + y} )

or

L_{2} = (∂_{y} + \frac{α}{x + y} ) (∂_{x} + \frac{β}{x + y} ) .

We note that in the particular case γ = 0, the Poisson operator reduces to the operator considered in Example 4 from [4] and the results there can be obtained immediately from the formulas above.

• Example 5.1.2. Parabolic operator The operator

A_{p} = a_{20}∂_{x}^{2} + a_{11}∂_{x}∂_{y} + a_{02}∂_{y}^{2} + a_{10}∂_{x} + a_{01}∂_{y} + a_{00}

is parabolic only if a_{11}^{2}−4a_{20}a_{02} = 0. This condition means that the characteristic polynomial

P_{2} := a_{20}ω^{2} + a_{11}ω + a_{02}

has a double root ω = −a_{11}/a_{20} and it follows from the discussion in Sec. 2 that factorization requires solution of a Riccati equation for p_{3}.

• Example 5.1.3. Elliptic operator The operator

A_{e} = a_{20}∂_{x}^{2} + a_{11}∂_{x}∂_{y} + a_{02}∂_{y}^{2} + a_{10}∂_{x} + a_{01}∂_{y} + a_{00}

(with real coefficients) is elliptic by definition if a_{11}^{2}−4a_{20}a_{02} < 0, so that the roots of the characteristic polynomial P_{2} are complex. Therefore, the factorization is only possible over a field that contains C.

5.2 Role of zero term; order of factors

Consider the following second order operators in two variables:

A_{a} = ∂_{x}^{2} − ∂_{y}^{2} + x∂_{y} + y∂_{x} + \frac{1}{4}(y^{2} − x^{2}) + a,  (12)

where a is a real constant.

• Example 5.2.1 We look for a factorization in the form

A_{a} = (p_{1}∂_{x} + p_{2}∂_{y} + p_{3})(p_{4}∂_{x} + p_{5}∂_{y} + p_{6}).

Necessarily p_{1}p_{4} = 1, and with no essential loss of generality we take p_{1} = p_{4} = 1. Then it is not difficult to show that the factorization can only have one of two forms:

A_{a} = L_{a}^{+}L_{a}^{−} = (∂_{x} + ∂_{y} + b)(∂_{x} − ∂_{y} + c)  (13)

or

A_{a} = L_{a}^{−}L_{a}^{+} = (∂_{x} − ∂_{y} + b)(∂_{x} + ∂_{y} + c)  (14)

It follows from the results above, and can be checked directly as an exercise, that the following is true:

– The operator

A_{a} = ∂_{x}^{2} − ∂_{y}^{2} + y∂_{x} + x∂_{y} + \frac{1}{4}(y^{2} − x^{2}) + a,

does not have a factorization of either form (13) or (14) unless a = ±1.

– The operator

A_{1} = ∂_{x}^{2} − ∂_{y}^{2} + x∂_{y} + y∂_{x} + \frac{1}{4}(y^{2} − x^{2}) + 1

has a factorization of the form (13) but has no factorization of the form (14).

– The operator

A_{−1} = ∂_{x}^{2} − ∂_{y}^{2} + x∂_{y} + y∂_{x} + \frac{1}{4}(y^{2} − x^{2}) − 1

has a factorization of the form (14) but has no factorization of the form (13).

In fact

A_{1} = [∂_{x} + ∂_{y} + \frac{1}{2}(y − x)] [∂_{x} − ∂_{y} + \frac{1}{2}(y + x)];  (15)

and

A_{−1} = [∂_{x} − ∂_{y} + \frac{1}{2}(y + x)] [∂_{x} + ∂_{y} + \frac{1}{2}(y − x)].

This example shows that the existence of a factorization depends crucially on the term of order zero, and also that an operator may have a left factor of a certain form but no right factor of the same form, or conversely. As we have shown, in two variables the question of existence of a first order left factor of a given form can be settled, and the the factors calculated, in a systematic way using procedures that are simple and (in general) purely algebraic.

• Example 5.2.2 The question of existence of a first order right factor can be reduced to the question for a left factor by taking the formal transpose. As an illustration, let us find a right factor for operator A_{1} in this way. In order to take advantage of what is already known, we make a coordinate change y → −y, which converts A_{1} to

Ã_{1} = ∂_{x}^{2} − ∂_{y}^{2} − y∂_{x} − x∂_{y} + \frac{1}{4}(y^{2} − x^{2}) + 1.

The transpose is

Ã_{1}^{t} = ∂_{x}^{2} − ∂_{y}^{2} + y∂_{x} + x∂_{y} + \frac{1}{4}(y^{2} − x^{2}) + 1 = A_{1} = L_{1}^{+}(L_{1}^{−}).

Therefore

Ã_{1} = (L_{1}^{−})^{t}(L_{1}^{+})^{t} = (−L_{1}^{−})^{t}(−L_{1}^{+})^{t}.  (16)

Now according to (15),

(−L_{1}^{+})^{t} = −[∂_{x} + ∂_{y} + \frac{1}{2}(y − x)]^{t} = ∂_{x} + ∂_{y} − \frac{1}{2}(y − x).

Undoing the coordinate transformation y → −y converts this last operator to

∂_{x} − ∂_{y} + \frac{1}{2}(y + x).

According to (16), this is a right factor for A_{1}, which is consistent with (15).

5.3 Classes of factorizable LPDOs

Putting some restrictions on the a_{i,j} as functions of x and y, one can describe all factorizable operators in a given class of functions. Let us exhibit explicit factorizations for one case of hyperbolic operators:

A_{2} := ∂_{xx} − ∂_{yy} + a_{10}∂_{x} + a_{01}∂_{y} + a_{00},  (17)

i.e. a_{20} = 1, a_{11} = 0, a_{02} = −1.

In this case there are only two possible factorizations:

• case 1: if 4a_{00} = 2(∂_{x} + ∂_{y})(a_{10} + a_{01}) + (a_{10}^{2} − a_{01}^{2}) holds, then we have the following factorization:

A_{2} = (∂_{x} + ∂_{y} + \frac{a_{10} − a_{01}}{2} ) (∂_{x} − ∂_{y} + \frac{a_{10} + a_{01}}{2} )  (18)

• case 2: if 4a_{00} = 2(∂_{x} + ∂_{y})(a_{10} − a_{01}) + (a_{10}^{2} − a_{01}^{2}) holds, then we have the following factorization:

A_{2} = (∂_{x} + ∂_{y} + \frac{a_{10} + a_{01}}{2} ) (∂_{x} − ∂_{y} + \frac{a_{10} − a_{01}}{2} ) .  (19)

Let us consider as a simple example the case when the lower order coefficients a_{10}, a_{01}, a_{00} are linear functions of one variable x, i.e.

\begin{cases}a_{10} & = b_{10,1}x + b_{10,0}, \\ a_{01} & = b_{01,1}x + b_{01,0}, \\ a_{00} & = b_{00,1}x + b_{00,0}. \end{cases}

Then the condition of factorization (24) is

4b_{00,1}x+4b_{00,0} = 2b_{10,1}+2b_{01,1}+b_{10,1}^{2}x^{2}+2b_{10,1}b_{10,0}x+b_{10,0}^{2}−b_{01,1}^{2}x^{2}−2b_{01,1}b_{01,0}x−b_{01,0}^{2} which leads to

zero term:

4b_{00,0} = 2b_{10,1} + 2b_{01,1} + b_{10,0}^{2} − b_{01,0}^{2},

linear term:

4b_{00,1} = 2b_{10,1}b_{10,0} − 2b_{01,1}b_{01,0},

quadratic term:

0 = b_{10,1}^{2} − b_{01,1}^{2}.

The last equality gives immediately b_{10,1} = ±b_{01,1}. Let us consider, for instance, the case b_{10,1} = −b_{01,1} = b. Then

\begin{cases}4b_{00,0} = b_{10,0}^{2} − b_{01,0}^{2} \\ 2b_{00,1} = −b(b_{10,0} + b_{01,0}), \end{cases}

or in parametric form:

b_{01,0} = \sqrt{2}t_{1} sinh t_{2},
b_{10,0} = \sqrt{2}t_{1} cosh t_{2},
b_{00,0} = \frac{1}{2}t_{1}^{2},
b_{00,1} = −\sqrt{2}t_{1}t_{3}(sinh t_{2} + cosh t_{2}),
b_{10,1} = t_{3}
b_{01,1} = −t_{3}

with arbitrary real parameters t_{1}, t_{2}, t_{3} ∈ R. Varying these parameters one can get coefficients b_{01,0}, b_{10,0}, ..., b_{00,0} and, therefore, coefficients of the initial operator a_{10}, a_{01}, a_{00}, that belong to some prescribed field. For instance, with t_{1}, t_{3} ∈ Q and t_{2} = 0,

\begin{cases}b_{01,0} = 0, \\ b_{10,0} = \sqrt{2}t_{1}, \\ b_{00,0} = \frac{1}{2}t_{1}^{2}, \\ b_{00,1} = −\sqrt{2}t_{1}t_{3}, \\ b_{10,1} = t_{3} \\ b_{01,1} = −t_{3} \end{cases}

and we get 2-parameter factorization in Q(\sqrt{2}):

A_{2} = (∂_{x} + ∂_{y} + \frac{2t_{3}x + \sqrt{2}t_{1}}{2} ) (∂_{x} − ∂_{y} + \frac{\sqrt{2}t_{1}}{2} )

The case b_{10,0} = b_{01,0} can be treated in a way analogous to the previous one and case 2 in a way analogous to case 1. Combining all the results one will give a description of all factorizable LPDO of the form (22) with corresponding restrictions on the coefficients.

6 Conclusions

The elimination order-reduction procedure has been presented for explicit factorization of LPDO of arbitrary order. It was shown that in generic case, i.e. in case when characteristic polynomial has a distinct root, this procedure is wholly algebraic and no differential equation has to be solved. Otherwise the factorization problem is equivalent to solving some Riccati equation. It was also demonstrated that different types of LPDO have different factorization properties, namely, algebraic factorization of hyperbolic operators, when possible, could be done over R, while factorization of elliptic operators is only possible over C. In case of parabolic operator no algebraic factorization is possible.

7 Acknowledgements

The authors are very much obliged to the organizing committee of the Workshop “Integrable Systems” which provided an excellent opportunity to meet and work together. We also highly appreciate the hospitality of Centro Internacional de Ciencias, Cuernavaca, Mexico where most of this work was done. We express special gratitude to Prof. A. B. Shabat and Prof. S. P. Tsarev for very stimulating discussions during the writing of this paper. Author^{∗} gratefully acknowledges support of the Austrian Science Foundation (FWF) under projects SFB F013/F1304.

References

[1] J. Apel. Theory of Involutive divisions and an applications to Hilbert Functions, J. Symb. Comp., 25, pp.683-704 (1998)

[2] E. Beke. Die Irreducibilität der homogenen linearen Differentialgleichungen, Math. Annalen 45 , pp.278-300 (1894)

[3] M. Bronstein. An improved algorithm for factoring linear ordinary differential operators. Proc. ISSAC’94, J. von zur Gathen (ed.), pp.336-340 (1994)

[4] D. Grigoriev, F. Schwarz. Factoring and Solving Linear Partial Differential Equations. J. Computing 73 , pp.179-197 (2004)

[5] E. Landau. Ein Satz über die Zerlegung homogener linearer Differentialausdrücke in irreducible Faktoren. J. für die reine and angewandte Math. 124 , pp.115-120 (1901/1902)

[6] A. Loewy. Über vollständig reduzible lineare homogene Differentialgleichungen, Math. Annalen 62 , pp.89-117 (1906)

[7] A. Loewy. Über reduzible lineare homogene Differentialgleichungen, Math. Annalen 56 , pp.549-584 (1903)

[8] P. J. Olver. Applications of Lie groups to differential equations. 2nd ed. Springer-Verlag, New York (1993)

[9] J. F. Pommaret. Partial differential equations and Lie pseudogroups, Gordon & Breach, New York (1978)

[10] F.Schwarz. A factorization algorithm for linear ordinary differential equations. Proc. ISSAC’89, G.H. Gonnet (ed.), pp.17-25 (1989)

[11] S.P.Tsarev. An algorithm for complete enumeration of all factorizations of a linear ordinary differential operator. Proc. ISSAC’96 (1996), Y.N. Lakshman (ed.), pp.226-231 (1996)

[12] S.P.Tsarev. Factorization on linear partial differential operators and Darboux integrability of nonlinear PDEs, Poster at ISSAC’98 (1998)

[13] S.P.Tsarev. Factorization of overdetermined systems of linear partial differential equations with finite dimensional solution space, In Proc. 4th Int. Workshop on Computer Algebra Scient. Comput. (CASC-2001) V. Ganzha, E.Mayr, V. Vorozhtsov (eds), pp-529-539. Springer-Verlag (2001)

[14] Z.Li, F.Schwarz, S.P.Tsarev. Factoring systems of linear PDEs with finite dimensional solution space, J. Symbolical Computation, V.36, pp.443-471 (2003)