The Parametric Solution of Underdetermined linear ODEs

Thomas Wolf

Department of Mathematics, Brock University 500 Glenridge Avenue, St.Catharines, Ontario, Canada L2S 3A1

Abstract

The purpose of this paper is twofold. An immediate practical use of the presented algorithm is its applicability to the parametric solution of underdetermined linear ordinary differential equations (ODEs) with coefficients that are arbitrary analytic functions in the independent variable. A second conceptual aim is to present an algorithm that is in some sense dual to the fundamental Euclid’s algorithm, and thus an alternative to the special case of a Gröbner basis algorithm as it is used for solving linear ODE-systems. In the paper Euclid’s algorithm and the new ‘dual version’ are compared and their complementary strengths are analyzed on the task of solving underdetermined ODEs. An implementation of the described algorithm is interactively

accessible under http://lie.math.brocku.ca/crack/uode.

Key words: control theory, non-commutative Gröbner Bases, differential Gröbner Bases, linear ordinary differential equations, underdetermined differential equations 2000 MSC: 34H05, 49J15, 13P10, 68W30

1. Introduction

Underdetermined ordinary and partial differential equations (ODEs/PDEs) are typical objects of investigation in control theory but also overdetermined systems of equations, where originally the number of equations is larger than the number of functions, may turn into underdetermined problems in the course of their solution. Examples are the conditions for Lie-symmetries and conservation laws of linearizable differential equations. These overdetermined systems of conditions have solutions that involve arbitrary functions of one or more variables, i.e. in the process of solving these systems towards its conclusion, underdetermined problems often occur.

Email address: twolf@brocku.ca (Thomas Wolf).

Preprint submitted to Elsevier Science 18 May 2009

Algorithms for solving underdetermined linear equations and systems (ODEs, PDEs, multidimensional discrete,..) are known ([1],[2],[3]). There also exists efficient software by Daniel Robertz et al. [4], [5]. The purpose of this paper is to describe an alternative algorithm that is elegant, efficient and in some sense complementary to the fundamental Euclid’s algorithm which is the basis for the Gröbner basis algorithms that are used in other implementations. An open question is whether the new algorithm can be generalized to partial differential equations (PDEs) and thus not only be a complimentary algorithm for the non-commutative differential Gröbner basis algorithm for solving ODEs but also for solving PDEs.

An implementation of algorithms described in this paper is accessible online under [7] and is part of the package Crack ([8]) for solving overdetermined systems. As we will describe, a new feature of these programs are flags which allow to prevent denominators or that allow to reduce the number of terms in the solution. This can be achieved by absorbing explicit x-dependent factors into new functions that are introduced during the computation.

We consider a single linear ODE

0 = ∑_{i}^{r} D_{i}f_{i} + c_{00}(x), r > 1,  (1)

D_{i} = ∑_{j=0}^{n_{i}} c_{ij}(x) (\frac{d}{dx})^{j} , n_{i} > 0

for functions f_{1} … f_{r} of the independent variable x. The equation may be homogeneous (c_{00} = 0) or inhomogeneous (c_{00} ≠ 0). With the requirement that the ODE is underdetermined we only assume that at least two functions f_{i} are involved, i.e. r > 1. The differential order of each f_{i} is n_{i}. For simplicity we assume that all coefficients c_{ij} are sufficiently often differentiable, at least ∑_{i} n_{i} times. An upper index in round brackets indicates the number of differentiations, i.e. \frac{d^{n}f_{i}}{dx^{n}} = f_{i}^{(n)} and low order derivatives are denoted by apostrophe, e.g. \frac{d^{2}f_{i}}{dx^{2}} = f′′_{i}.

The task is to find the general solution of the ODE (1) in the form of explicit differential expressions F_{i}

f_{i} = F_{i}(x, h_{1}(x), …, h_{r−1}(x)), i = 1, …, r (2)

in terms of parametric functions h_{1}(x), …, h_{r−1}(x) which are either all free or one of them, say h_{1}(x), having to satisfy an ODE and the other h_{j}(x) being free.

For example, the general solution of the ODE

0 = f′′x^{2} + g′′x − g′x^{2} + f + 3x

for f = f(x), g = g(x) can be written in the form

f = \frac{x}{x^{8} − 2x^{6} + 7x^{4} − 6x^{2} + 9} ((x^{5} − x^{3} + 3x)h′′ − (x^{6} + x^{4} + 3x^{2} − 6)h′  (3)
+(3x^{5} + 3x^{3} + 17x)h − 3x^{8} + 3x^{6} − 16x^{4} + 9x^{2}) ,

g = \frac{x}{2(x^{8} − 2x^{6} + 7x^{4} − 6x^{2} + 9)} ((−2x^{6} + 2x^{4} − 6x^{2})h′′  (4)
+(8x^{5} − 4x^{3})h′ − (14x^{4} + 14x^{2} + 6)h + 4x^{7} + x^{5} + 3x^{3} − 27x)  (5)

where h = h(x) is an arbitrary function. The form of the solution is not unique. Because h(x) is arbitrary, replacing h by a differential expression in one or more arbitrary functions would give a solution too, but this solution would involve expressions of higher differential order and more arbitrary functions without being more general than the original solution. One naturally seeks a general solution which involves only as few as possible arbitrary functions and lowest order derivatives of them. But even this requirement does not give a unique solution. One might want to minimize the highest order of all parametric functions or, for example, the sum of all orders of all parametric functions.

Furthermore, the occurring function(s) could be scaled to modify the form of the solution, for example, to make it denominator free. In the above case replacing h(x) by (x^{8} − 2x^{6} + 7x^{4} − 6x^{2} + 9)^{3}p(x) with an arbitrary function p(x) would make the solution polynomial but increase its size, i.e. the total number of terms in the coefficients of the parametric functions on the right hand sides.

In the following section the algorithm is given, followed by comments on its characteristics. In section 3 and 4 we formulate the new algorithm and Euclid’s algorithm in vector notation in order to compare them in detail in section (5). We show that both are essentially different but also that they compliment each other in the sense that one can give a criterion under which circumstances which of both methods is better suited. Another option is to combine both algorithms in a hybrid version.

Finally, in section (7) an application is described which arose from a classification of hyperbolic evolutionary vector PDEs.

2. The Algorithm

2.1. Outline

The essence of the algorithm is

• to partition the ODE into a total derivative and an algebraic remainder,

• to introduce a new function f_{r+1} such that the total derivative part is df_{r+1}/dx and thus to write the ODE as a system of 2 equations: one equation defining f_{r+1} and one re-formulating the ODE in terms of f_{r+1},

• to use the re-formulated ODE to eliminate and substitute another function and thus to arrive again at a single ODE for the same number of functions which in some sense is closer to be solvable algebraically than the ODE before.

These steps are repeated until either the ODE involves only one function and thus is not underdetermined anymore, or until one function occurs purely algebraically and thus allowing the ODE to be solved for that function. The following is a more detailed pseudo code description.

2.2. Pseudocode

Input • list of functions f_{1}, …, f_{r} of x,

• linear ODE 0 = ω in f_{i} (like (1)) Body L := {} % L will be a list of substitutions s := r while (the ODE 0 = ω involves at least two f_{i}) and

(differential orders n_{i} > 0 ∀f_{i}) do

• Factor out \frac{d}{dx} from ω once as far as possible by introducing a new function f_{s+1}(x) and by computing expressions b_{i}:

0 = ω = fs+1^{′} + ∑_{i=1}^{s} f_{i}b_{i} + a_{00} (6)

where f_{s+1}, b_{i} are given through:

f_{s+1} = ∑_{i=1}^{s} ∑_{j=0}^{n_{i}−1} f_{i}^{(n_{i}−1−j)} ∑_{k=0}^{j}(−1)^{(j+k)}a_{i,n_{i}−k}^{(j−k)}  (7)

b_{i} = ∑_{k=0}^{n_{i}}(−1)^{(n_{i}+k)}a_{i,n_{i}−k}^{(n_{i}−k)}  (8)

• if b_{i} = 0 for all i then the ODE 0 = ω is exact (apart from a_{00}), i.e. consider new ODE 0 = ω̂ := f_{s+1} + ∫ a_{00} dx where f_{s+1} is only an abbreviation defined in (7) else ∘ regard f_{s+1} as a new unknown function, ∘ solve (6) for one function f_{j} of f_{1}, .., f_{s} that

− has a non-vanishing coefficient b_{i} in (6), and − is of lowest possible order in ω and thus in (7),

f_{j} = −\frac{1}{b_{j}} (fs+1^{′} + ∑_{i≠j} f_{i}b_{i} + a_{00}) (9)

∘ use (9) to substitute f_{j} in (7) to get a new ODE 0 = ω̂ := −f_{s+1} + ∑_{i=1}^{s} ∑_{j=0}^{n_{i}−1} … for functions f_{1}, .., f_{j−1}, f_{j+1}, .., f_{s}, f_{s+1}

∘ update ω := ω̂, s := s + 1, L := {f_{j} = ..} ∪ L end if the ODE involves a function f_{j} purely algebraically then solve for f_{j} and add it to L : L := {f_{j} = ..} ∪ L h(x) := f_{s}(x) % to remember the last introduced function P := L % to remember the complete list of substitutions L % finally: substitutions (each as stored in the first % element of P) are performed in the rest of P while s > r do P := rest(P)|_{f_{s}=... % as given in first(P)} s := s − 1 end end % of Body

Output • l % list of new functions • L % the complete list of substitutions

• P % list of all initial functions f_{1}, .., f_{r} that are not % parametric, each given in terms of h • 0 = ω(x, h(x)) % only if the first while loop

% terminates due to ω not

% being underdetermined anymore

2.3. Comments

Before comparing the algorithm with Euclid’s algorithm a few comments are necessary.

• All steps involve only algebra or differentiations with only one exception: if the homogeneous part of the ODE is exact then the integral of the inhomogeneous part is taken. But this integral does not have to be evaluated, i.e. to be expressed in terms of elementary functions. It can stay in a symbolic unevaluated form for the algorithm to continue.

• The coefficients a_{ij} can be arbitrary explicit functions of x, and involve, for example, sin or log with the only condition that it must be decidable whether expressions involving these functions and their derivatives are zero or not.

• The transformation of the ODE in each step of the first while loop is reversible, i.e. the new ODE is equivalent to the previous one. Therefore the obtained solution is the general one.

• The algorithm terminates.

Proof: We consider how the total sum ∑_{i=1}^{s} n_{i} of differential orders n_{i} of all functions f_{i} changes during execution. After substitution of f_{j} with (9) in (7) the new differential orders n̂_{i} of functions f_{i} that occur in the new ODE are:

n̂_{i} \{ \begin{matrix}= & n_{j} & for i = s + 1 \\ ≤ & max(n_{j} − 1, n_{i} − 1) = n_{i} − 1 & for i ≤ s, b_{i} ≠ 0 \\ = & n_{i} − 1 & for i ≤ s, b_{i} = 0. \end{matrix}

The function f_{j} of order n_{j} gets replaced by a new function f_{s+1} which then occurs with same order n_{j}, but the order of all other functions is lowered by at least one because we choose j such that n_{j} = min(n_{k}) (∀k with b_{k} ≠ 0), i.e. we have n̂_{i} ≤ n_{i} − 1 also for the f_{i} with b_{i} ≠ 0. Because the ODE has at least two functions, the total sum of derivatives is decreasing. The algorithm is therefore finite.

• The algorithm results in an ODE for a single function iff the differential operators D_{i} in (1) have a common differential factor. The remaining ODE is of the same order as the common factor, i.e. its order can not be higher than the order of the original ODE (1).

• The algorithm naturally splits into two parts, a part A) establishing a list L of substitutions (9) in the first while loop and part B) performing the substitutions in the second while loop to obtain an explicit solution. The first part is executed very fast (see next section for more details), the second may take longer for higher order ODEs because expressions typically grow with each substitution, often exponentially. An example is given in the appendix.

For many applications the list L of substitutions may even be of higher practical value than the explicit formulas f_{i} = F_{i}(x, h_{1}(x), …, h_{r−1}(x)), i = 1, …, r resulting from B). In general, the list L is a much shorter representation of the parametric solution of the ODE than the explicit solution itself and thus is more useful as a solution, like in the following scenario.

Let us assume that we have a large algebraic expression in terms of the original functions f_{1}, .., f_{r} that has to be evaluated modulo the underdetermined ODE (1). Instead of replacing the f_{i} directly by their large explicit expressions as given in the list P it is usually much better to perform successively the substitutions stored in L in the order they were derived, allowing cancellations to happen after each individual substitution. Also, if numerical computations are to be done, it is much faster to compute the sequence of substitutions than the explicit expressions.

3. A vector representation

To perform the iteration process efficiently, i.e. to replace functions by algebraic combinations of other functions (9) in the ODE (7) without needing to perform any differentiations, the ODE has to be represented in a form where D := \frac{d}{dx} is factored out as far as possible. As we will see further down, this is also the appropriate representation for using the Euclidean algorithm to solve the underdetermined ODE. We write the ODE in the form

0 = ∑_{i}^{r} A_{i}f_{i} + a_{00}(x)  (10)
A_{i} = DÃ_{i} + a_{i0}(x)  (11)
Ã_{i} = D^{n_{i}−1}a_{in_{i}}(x) + … + Da_{i2}(x) + a_{i1}(x)  (12)

Replacing, for example, f_{1} = w(x)f̂_{1} would only require multiplications â_{1k} := wa_{1k} to update this representation and no differentiations.

In this notation a single step in the first while loop of the new method consists of • introducing a new function f_{r+1} through

f_{r+1} = ∑_{i} Ã_{i}f_{i} (13)

giving the ODE the form

0 = Df_{r+1} + ∑_{i} a_{i0}(x)f_{i} + a_{00}. (14)

• regarding the defining relation (13) for the new function as the new ODE and using the old ODE (14) for a substitution of an ‘old’ function. For that we choose one a_{j0} of the non-vanishing a_{i0} in (11) for which the order n_{j} of the corresponding f_{j} is minimal. Thus, performing the substitution

f_{j} = −\frac{1}{a_{j0}} (Df_{r+1} + ∑_{i≠j} a_{i0}f_{i} + a_{00}) (15)

in the definition (13) gives the new ODE

0 = −f_{r+1} + ∑_{i≠j} Ã_{i}f_{i} − ∑_{i≠j} Ã_{j} (\frac{a_{i0}}{a_{j0}} f_{i}) − Ã_{j} (\frac{1}{a_{j0}} Df_{r+1}) − Ã_{j} (\frac{a_{00}}{a_{j0}} ) . (16)

In vector notation the above two steps change the differential operators (assuming for simplicity of notation j = 1) according to:

( \begin{matrix}A_{1} \\ A_{2} \\ : \\ A_{r} \end{matrix} ) → ( \begin{matrix}A_{r+1} \\ A_{2} \\ : \\ A_{r} \end{matrix} ) = ( \begin{matrix}−1 − Ã_{1} \frac{1}{a_{10}} D \\ Ã_{2} − Ã_{1} \frac{a_{20}}{a_{10}} \\ : \\ Ã_{r} − Ã_{1} \frac{a_{r0}}{a_{10}} \end{matrix} ) (17)

a_{00} → −Ã_{1} (\frac{a_{00}}{a_{10}} ) (18)

and add (15) to the list of substitutions. Updating components 2 … r in (17) is done purely by multiplications and additions, only updating the first component takes differentiations when factoring out D.

A more conventional algorithm solving underdetermined linear ODEs tries to lower the order of the ODE by performing Euclid’s algorithm until the ODE is algebraic for one function and thus can be solved algebraically. In the following section we formulate such an Euclidean step to be able to compare it with the new algorithm above.

4. The Euclid version in vector notation

For applying the (right) Euclid algorithm the ODE is given as well in a form with D completely factored out:

0 = ∑_{i}^{r} A_{i}f_{i} + a_{00}(x)  (19)
A_{i} = D^{n_{i}}a_{in_{i}}(x) + … + Da_{i1}(x) + a_{i0}(x).  (20)

One iteration step is performed by

• choosing two functions f_{i}, f_{j} (w.l.o.g. n_{i} ≥ n_{j}) and introducing a new function f_{r+1} through

f_{j} = f_{r+1} − D^{n_{i}−n_{j}} (\frac{a_{in_{i}}}{a_{jn_{j}}} f_{i}) (21)

• and performing this substitution in the ODE (19). The substitution is chosen such that the differential order of f_{i} is lowered by at least one. This iteration process stops like in the new algorithm when the first function appears purely algebraically or when the ODE involves only one function. Therefore, in order to minimize the number of steps one would choose f_{j} as one of the lowest order functions and f_{i} as one of the lowest order functions from the remaining ones. (From all these choices it is recommendable to choose a pair so that multiplication is minimized and A_{i}, A_{j} have as few as possible terms.) For j = 1, i = 2 the vector notation gives

( \begin{matrix}A_{1} \\ A_{2} \\ A_{3} \\ : \\ A_{r} \end{matrix} ) → ( \begin{matrix}A_{r+1} \\ A_{2} \\ A_{3} \\ : \\ A_{r} \end{matrix} ) = ( \begin{matrix}A_{1} \\ A_{2} − A_{1}D^{n_{2}−n_{1}} \frac{a_{2n_{2}}}{a_{1n_{1}}} \\ A_{3} \\ : \\ A_{r} \end{matrix} ) (22)

a_{00} → a_{00}. (23)

The update of the 2^{nd} component in (22) requires differentiations when D is factored out as far as possible.

5. Relations between both algorithms

5.1. Differences

In this section we want to justify our claim that both algorithms differ significantly from each other and that they are not merely variations of one and the same procedure.

Both algorithms differ conceptually in that Euclid’s algorithm pairs two differential operators to get a new differential operator of lower order whereas in the new algorithm the relation that defines a new function becomes the ODE and the old ODE is used for substituting one of the ‘old’ functions.

Another difference is that Euclid’s algorithm lowers in one step the order of only one differential operator whereas the new algorithm lowers the order of all but one differential operators. Both types of steps are performed at about the same cost as for the new method only the update of the first component in (17) is potentially size increasing and for Euclid’s method it is only the update of the second component of (22). Using Euclid’s algorithm one can of course take one of the lowest order differential operators and decrease the order of all others but that would take r − 1 computations, each potentially size increasing. This difference does usually not affect the number of necessary iterations (until one function occurs purely algebraically), because Euclid’s method can operate just on two of the lowest order operators and ignore all others. So, by applying Euclid’s method to lower the order of the two lowest order operators A_{i} with each other, all other higher order operators keep their high order. Moreover, when the algorithm stops because the final equation contains one function algebraically, then back substitutions start (part B of the algorithm) which even increase the highest orders.

Differently with the new method where the differential order of all functions but one gets lowered in one step. Thus the obtained parametric solution is typically of lower order than the solution obtained by Euclid’s algorithm if this is executed by pairing only lowest order derivatives.

For example, the size of the parametric solution of the ODE

0 = x^{3}f_{ax} + (x − 1)g_{bx} + h_{5x}, for f, g, h of x (24)

and the differential order of parametric functions in it depend strongly on the differential orders a, b of f and g in the ODE and the method that is used. As a, b increase from a = b = 1 to a = b = 3 one can see the following trends:

• In Euclid’s solution the differential order of parametric functions increases from 5 to 7 whereas in the solution of the new algorithm the order decreases from 4 to 3 in the expression for g, from 4 to 2 in the expression for f and only increases from 0 to 1 in the expression for h.

• The size of expressions in Euclid’s solution steadily increases and in the solution of the

new method it decreases.

For a = b = 1 the solutions have comparable size:

- When explicit x-dependent factors are not absorbed (for absorbing factors see section 5.3) then h is parametric and the rational expressions for f, g have the form (5 terms)/(2 terms) in Euclid’s solution and in the solution of the new method f = (4 terms)/(2 terms), g = (3 terms)/(2 terms).

For a = b = 3:

- When factors are not absorbed then in Euclid’s solution h is parametric, f = (85 terms)/(15 terms), g = 13 terms and in the solution of the new algorithm we have f = 1 term, g = 7 terms, h = 4 terms.

- When factors are absorbed then Euclid’s solution gives f = 27 terms, g = 16 terms and in the solution of the new algorithms we have f = 1 term, g = 7 terms, h = 4 terms. Not only the solutions of both algorithms may differ significantly but also the number

of steps to reach them. Starting, for example, with

0 = f′ + a(x)h^{(n)} = 0

to be solved for f = f(x), h = h(x) this takes only one iteration step with the new algorithm giving, for example for n = 5, a solution of the form f = (27 terms)/(1 term), h = (1 term)/(1 term), whereas it takes n steps for Euclid’s algorithm giving, for example for n = 5, a more spacious solution of the form f = (40 terms)/(1 term), h = (2 terms)/(1 term).

If the ODE is inhomogeneous then another difference between both algorithms becomes apparent. With non-zero inhomogeneity a_{00} in the ODE (10) the transformation (15) becomes inhomogeneous and changes a_{00} in (18) whereas Euclid’s algorithm does not change a_{00} in the homogeneous transformation (21).

Two ODEs that highlight the duality between both algorithms are

0 = f′ + f + g′ + (ah)^{(20)}  (25)
0 = f′ + f + g′ + a(h)^{(20)}  (26)

where a = a(x) is a given function of x and the ODE has to be solved for f, g, h. As shown in table 1 for the ODE (25) the new algorithm is to be preferred whereas for ODE (26) Euclid’s algorithm gives a shorter solution.

algorithm  equation (25)  equation (26)
new  f=1 term  g=22 terms f=2 terms g=23 terms
Euclid  f=22 terms g=23 terms f=2 terms g=3 terms

Table 1. The size of solutions given by both algorithms (both executed without absorbing of factors).

An explanation of the above behaviour comes from the fact that Euclid’s method is better suited when factors \frac{a_{in_{i}}}{a_{jn_{j}}} are small, i.e. when the coefficients of the leading derivatives in the representation (10) - (12) effectively cancel each other in a quotient, and less suited when these factors are large prime expressions in x. Differently, the new method works best if factors \frac{a_{i0}}{a_{j0}} are small, i.e. when the coefficients of the algebraic terms effectively cancel each other, and less suited when these factors are large prime expressions in x.

Although both ODEs look similar, they are rather different as (ah)^{(20)} has many terms if the product rule of differentiation is fully applied and a(h)^{(20)} has many terms if \frac{d}{dx} is factored out as far as possible. Thus both algorithms differ in their suitability for both ODEs.

The above computations and any other tests can easily be performed online (see [7]), where access to the computer algebra system REDUCE and to the procedures uode and print uode solution is provided. The procedures have been written and tested in the PSL (Portable Standard Lisp) version of REDUCE. At the time of submitting this paper (May 2009) the CSL version of REDUCE is freely downloadable from [9].

5.2. Hybridization

The advantage of having these two, in some sense, complementary algorithms lies in the fact that both operate on the same data structure, i.e. both their input and output consists of an ODE in the form (10) - (12) and both add one substitution to a list of substitutions of functions in terms of newly introduced functions. Consequently, both algorithms are interchangeable, i.e. one step performed with one algorithm could be followed by another step performed with the other algorithm, in order to minimize the size of x-dependent factors and thus to minimize growth. A different strategy could be to perform each step with both algorithms in parallel, to compare the size of the ODE and/or size of derived substitution that both methods give and to choose which one to adopt for this step. This strategy does at most double the amount of computation but more likely will lead to savings due to working with smaller expressions.

Any hybrid algorithm performing a mixture of Euclid steps and new algorithm steps is still finite because in any such step the sum of differential orders of all functions is decreasing.

5.3. Absorbing factors

The following efficiency improving measures work for both methods. A representation (10) - (12) (and identically (19) - (20)) of the ODE where D is maximally factored out has the advantage that any change of functions f_{i} → h(x)f̄_{i} does require only multiplications and is done very easily in both algorithms. This freedom of efficiently multiplying functions with x-dependent factors can be used for different purposes.

If in a newly generated ODE all coefficients of a function f_{i} have a non-trivial GCD ĉ(x) := GCD(a_{in_{i}}, …, a_{i0}) ≠ 1 then this can be absorbed into a new function f_{r+1} and the ODE be simplified by performing the substitution

f_{i} = f_{r+1}/ĉ (27)

and adding it to the accumulating list of substitutions. Such non-trivial GCDs occur relatively frequently as the GCD is taken only from the few coefficients a_{ij} of any single one operator A_{i} that changed in an iteration step.

Another purpose of introducing new functions multiplied with an explicit x-depending factor is to avoid a denominator (den). This would arise in the first component of (17). It can be prevented by introducing another new function f_{r+2} and by performing the substitution

f_{r+1} = a_{10}^{2}f_{r+2}. (28)

Similarly, in the other components of (17) one can scale

f_{i} = den (\frac{a_{i0}}{a_{10}} ) f_{r+i}. (29)

In Euclid’s algorithm the appearance of denominators in the second component of (22) can be prevented by introducing the scaling

f_{2} = den (\frac{a_{2n_{2}}}{a_{1n_{1}}} ) f_{r+1}. (30)

Another situation where a dominator occurs is at the end of both algorithms when one function occurs purely algebraically

0 = a_{i0}f_{i} + ∑_{j≠i} A_{j}f_{j} + a_{00}(x). (31)

To avoid the denominator a_{i0} one can scale all f_{j} depending on their differential order and on their coefficients in A_{j}. Avoiding this denominator is especially helpful, as f_{i} = … is the last substitution and thus substituted successively in all previous substitutions leading easily to a substantial growth of denominators. By avoiding denominators in the above way both methods produce denominator free solutions if the ODE is homogeneous.

5.4. Embedded ODEs

The case that an underdetermined ODE factorizes, i.e. that the differential operators A_{i} in (10) have a non-trivial GCD, or in other words, that the ODE can be written in form of two nested ODEs 0 = Ω(x, ω(x, f_{i})) is discovered by both algorithms. The inner ODE 0 = ω is only determined up to a linear change ω = α(x)ω̂ and both algorithms will usually find Ω’s that differ by some α(x). A slight advantage of the new algorithm is that the D^{0} part of the ODE is used within the algorithm, thus it is automatically recognized if this vanishes, i.e. if the ODE is exact (up to an inhomogeneity).

6. ODE-systems

The introduced methods of solving a single underdetermined ODE (or converting it to an ODE for a single function) can be used to convert an ODE system into an equivalent set of fully decoupled ODEs, each for a single function.

After treating any one equation of the original ODE system the original functions in it are expressed in terms of fewer functions which are either all free or at most one has to satisfy a single ODE. All functions can be replaced in the remaining ODEs. This can go on as long as we have equations of at least two functions. When the procedure stops we have solved the system or got ODEs, each containing only one function. If more than one ODE contain only one and the same single function, then these ODEs form an over-determined subsystem which can be treated by a Gröbner bases computation and result either in the explicit solution for this function or a single ODE for this function of an order not higher than the lowest order of the ODEs for this function. This whole procedure terminates with either the explicit parametric solution of the original system or a decoupled set of ODEs, each ODE for a single function.

7. An application

A class of applications where underdetermined linear ODEs occur frequently is the classification of hyperbolic evolutionary PDE systems. The aim of such an investigation is to find integrable systems of PDEs by determining those systems which have a higher order symmetry (see below). More information about the mathematical background is given in [6] where a classification of hyperbolic vector PDEs is discussed.

Let us look at two hyperbolic scalar PDEs for functions u(x, t, τ), v(x, t, τ) (where x, t are the usual independent variables and τ is a symmetry parameter). The following ansatz for the system and symmetry is generated based on homogeneity considerations. By requiring the same homogeneity weights as the potential nonlinear Schrödinger equation we get for the system the ansatz

u_{tx} = a_{10}u_{3x} + a_{6}u_{x}v_{x}^{2}u^{2} + a_{8}u_{2x}u_{x}v + a_{9}u_{2x}v_{x}u + a_{1}u_{x}^{2}v_{x} + a_{3}u_{x}^{3}v^{2} + a_{7}u_{x}v_{2x}u
v_{tx} = a_{19}v_{3x} + a_{15}u_{x}^{2}v_{x}v^{2} + a_{18}v_{2x}v_{x}u + a_{17}u_{x}v_{2x}v + a_{13}u_{x}v_{x}^{2} + a_{16}v_{x}^{3}u^{2}.  (32)

Assuming the same differential order for the symmetry we get the ansatz

u_{τ} = b_{10}u_{3x} + …, v_{τ} = b_{19}v_{3x} + …, (33)

with right hand sides identical to those of (32), only with coefficients b_{i} instead of a_{i}. All coefficients a_{i}, b_{i} are undetermined functions of the product uv.^{1}

The relations (33) are considered to be a symmetry of the system (32) if the symmetry conditions

∂_{τ} (u_{tx}) − ∂_{t}∂_{x}(u_{τ}) = 0, ∂_{τ} (v_{tx}) − ∂_{t}∂_{x}(v_{τ}) = 0 (34)

are fulfilled identically in u, v and in derivatives of u and v, all modulo substitutions based on (32) and (33). In performing the differentiations in (34), doing repeatedly substitutions (32), (33) and finally setting all coefficients of different products of powers of derivatives of u, v individually to zero gives 27 ODEs with a total of 1334 terms for 26 functions a_{i}, b_{j} of z := uv. The length of equations ranges from 2 to 266 terms. In the course of solving this system the program Crack performs integrations, substitutions, splittings (i.e. separations when z occurs only explicitly) and a number of case distinctions. In one of the sub cases the resulting conditions can be integrated successively up to the underdetermined linear ODE

0 = 3b13^{′}z − 6b15^{′}z^{2} − 2b17^{′′}z^{2} + b17^{′}z − 6b_{15}z + 2b_{17}. (35)

Step 1: To start, partial integration gives

0 = (3b_{13}z − 6b_{15}z^{2} − 2b17^{′}z^{2} + 5b_{17}z)′ − 3b_{13} + 6b_{15}z − 3b_{17}. (36)

which can be written as

0 = c1^{′} − 3b_{13} + 6b_{15}z − 3b_{17} (37)

by introducing c_{1}(z) through

c_{1} = 3b_{13}z − 6b_{15}z^{2} − 2b17^{′}z^{2} + 5b_{17}z. (38)

From the functions that occur only algebraically in (37) (i.e. b_{13}, b_{15}, b_{17}) the ones that have lowest derivatives in (35) are b_{13}, b_{15}. Eliminating one of them from (37), say b_{13} and substituting it in (38) gives

0 = 2b17^{′}z^{2} − 2b_{17}z − c1^{′}z + c_{1} (39)

which is not algebraic in any function yet, but already of first order, so one more step has to be performed.

1 For hyperbolic systems the homogeneity weights include negative values, for example, in this case weight(u)=1, weight(v)=-1, so that uv has weight zero and therefore we have no limitation on the degree of powers of uv and thus all unknown coefficients are arbitrary functions of uv.

Step 2: Partial integration of (39) and introduction of

c_{2} = 2b_{17}z^{2} − c_{1}z (40)

results in

0 = c2^{′} − 6b_{17}z + 2c_{1} (41)

which allows to solve for c_{1} and replace it in (40) giving

0 = \frac{1}{2}c2^{′}z − c_{2} − b_{17}z^{2}. (42)

This condition is purely algebraic for one function, b_{17}, and therefore the algorithm stops.

Cleanup: We need explicit solutions of (35), so what remains to be done are back substitutions: (42) provides

b_{17} = \frac{1}{2z} c2^{′} − \frac{1}{z^{2}} c_{2}. (43)

The second substitution expressing b_{13} in terms of the parametric function c_{2} is obtained after backward substituting c_{1} from (41) into (37) and b_{17} from (43) in (37) to get the explicit solution consisting of (43) and

b_{13} = \frac{1}{3}c2^{′′} − \frac{3}{2z} c2^{′} + 2b_{15}z + \frac{2}{z^{2}} c_{2} (44)

involving the free function c_{2}(z).

Another underdetermined equation resulting in this integrability problem is

0 = 3b1^{′}z − 9b3^{′}z^{2} + 3b6^{′}z^{2} − 4b8^{′′}z^{2} − 12b8^{′}z − 12b_{3}z + 6b_{6}z − 3b_{8} (45)

which has the solution

b_{1} = (−8b8^{′}z^{2} + 3c3^{′}z + 6b_{6}z^{2} − 5b_{8}z − 2c_{3})/(3z)
b_{3} = (−4b8^{′}z^{2} + c3^{′}z + 3b_{6}z^{2} − 3b_{8}z − c_{3})/(3z^{2}).

We finally obtain as a system with higher order symmetries:

u_{tx} = \frac{a_{1}}{2uv^{3}} ((uv)_{x}^{3}v − uv(uv)_{x}^{2}v_{x})
v_{tx} = \frac{a_{13}}{2uv} (uv)_{x}^{2}v_{x}.

8. Appendix

The following example illustrates the comments made in the last paragraphs of section 2.3. It shows two possible representations of a solution to an underdetermined ODE. For the equation

(x − 1)^{3}f_{1}^{(5)} + 3f_{1}^{(3)} + xf′′_{1} + (1 − x^{2})f′_{1} + f_{1} − (x − 2)(x − 3)f′′_{2} − xf′_{2} = 0

the generated list L of substitutions is given through

f_{6} = (7289/9f′_{8}x^{6} − 39400/3f′_{8}x^{5} + 725945/9f′_{8}x^{4} − 691667/3f′_{8}x^{3} + 2665016/9f′_{8}x^{2}
−463541/3f′_{8}x + 40582f′_{8} − 82543/9f_{8}x^{5} + 393212/3f_{8}x^{4} − 1997777/3f_{8}x^{3}
+1404610f_{8}x^{2} − 9064364/9f_{8}x + 863648/3f_{8})/(x^{1}2 − 27x^{1}1 + 339x^{1}0 − 2551x^{9}
+12566x^{8} − 43294x^{7} + 111667x^{6} − 221121x^{5} + 332143x^{4} − 447477x^{3} + 625912x^{2}
−422746x + 117948)
f_{7} = (−197/419f′_{8} − 2227/6704f_{6}x^{6} + 40655/6704f_{6}x^{5} − 282747/6704f_{6}x^{4}
+937545/6704f_{6}x^{3} − 737203/3352f_{6}x^{2} + 236215/1676f_{6}x − 15636/419f_{6})/
(x^{5} − 5988/419x^{4} + 30423/419x^{3} − 64170/419x^{2} + 46012/419x − 13152/419)
f_{1} = (−288/197f′_{7} − 450/197f_{6})/(x^{5} − 2988/197x^{4} + 15683/197x^{3} − 32571/197x^{2}
+18956/197x − 6093/197)
f_{5} = −9/16f′_{6} + 37/32f_{1}x^{4} − 415/32f_{1}x^{3} + 347/8f_{1}x^{2} − 815/32f_{1}x + 261/32f_{1}
f_{4} = −4/9f′_{5} − 17/9f_{1}x^{3} + 124/9f_{1}x^{2} − 176/9f_{1}x + 11f_{1}
f_{3} = −1/4f′_{4} + 5/4f_{1}x^{2} − 17/4f_{1}x − 11/4f_{1}
f_{2} = f′_{3} + 2f_{1}x + f_{1}

which is a much shorter representation of the solution than the 10 page explicit form which results from substituting f_{6}, f_{7}, f_{1}, f_{5}, f_{4}, f_{3} in this order into each other and takes the form f_{1} =(42 terms)/(25 terms), f_{2} =(304 terms)/(61 terms) involving up to 29-digit integers. The list of 7 substitutions is not only shorter than the list of 2 substitutions for f_{1}, f_{2}, it also is much faster to derive and more useful if a differential expression is to be simplified modulo the solution of the above ODE by substituting f_{2}, f_{3}, f_{4}, f_{5}, f_{1}, f_{7}, f_{6} in this order.

The difference in size of both solution representations can be arbitrarily amplified by having an input ODE of higher order and higher degree polynomials as coefficients.

9. Summary

We present an algorithm for the solution of underdetermined linear ODE that is compatible but structurally different from the (right) Euclidean algorithm. Because both algorithms operate on the same data structure and because both complement each other in the sense that each one is most efficient for ODEs of different form, the combination of both algorithms is superior to each individual one.

Acknowledgements

I would like to thank Sergey Tsarev for many comments. Daniel Robertz is thanked for comparative runs with the OreModules package.

References

[1] Chyzak, F., Quadrat, A., Robertz, D. (2005). ”Effective algorithms for parametrizing linear control systems over Ore algebras”, Rapport de Recherche INRIA , Applicable Algebra in Engineering, Communications and Computing 16, no 5, 319-376.

[2] Pommaret, J.-F., Quadrat, A. (2004). ”A differential operator approach to multidimensional optimal control”, International Journal of Control 77, 821-836.

[3] Fröhler, S., Oberst, U. (1998). ”Continuous time-varying linear systems”, Systems & Control Letters 35, 97-110.

[4] Chyzak, F., Quadrat, A., Robertz, D. OreModules project

http://wwwb.math.rwth-aachen.de/OreModules

[5] Chyzak, F., Quadrat, A., Robertz, D. ”OreModules: A symbolic package for the study of multidimensional linear systems” in: J. Chiasson, J.-J. Loiseau, ”Applications of Time-Delay Systems”, Springer, to appear.

[6] Anco, S. and Wolf, T. (2005). ”Some symmetry classifications of hyperbolic vector evolution equations”, JNMP, Volume 12, Supplement 1, p 13-31 (also nlin.SI/0412015).

[7] Wolf, T. (2007). Online demo for solving underdetermined ODEs.

http://lie.math.brocku.ca/crack/uode

[8] Wolf, T.: Applications of Crack in the Classification of Integrable Systems, CRM Proceedings and Lecture Notes, vol 37 (2004) pp. 283-300. (arXiv nlin.SI/0301032) and online under [7].

[9] REDUCE - a portable general-purpose computer algebra system

http://reduce-algebra.sourceforge.net/