Size reduction and partial decoupling of systems of equations^{∗}

Thomas Wolf Department of Mathematics, Brock University, 500 Glenridge Avenue, St. Catharines, Ontario, Canada L2S 3A1 email: twolf@brocku.ca

November 24, 2010

Abstract A method is presented that reduces the number of terms of systems of linear equations (algebraic, ordinary and partial differential equations). As a byproduct these systems have a tendency to become partially decoupled and are more likely to be factorizable or integrable. A variation of this method is applicable to non-linear systems. Modifications to improve efficiency are given and examples are shown. This procedure can be used in connection with the computation of the radical of a differential ideal (differential Gröbner basis).

1 Motivation

Algorithms for applying integrability conditions to a system of differential equations in a systematic way in order to generate simplified differential equations are implemented in a number of programs ([3, 4, 7, 8, 10, 11, 12] and more in [6]). Such calculations result in the radical or a (pseudo) differential Gröbner Basis of the differential ideal generated by the original system. A common problem of these algorithms, and consequently their implementations, is an explosive expression swell. Optimizations like Buchberger’s 2^{nd} criterion (section 5.5 in [1]) and their analogue for differential equations aim to reduce the number of steps to reach a characteristic set (one step = computation of an S-polynomial for algebraic systems or a cross-differentiation of two differential equations for differential systems). These optimizations do not cover other ‘obvious’ simplifications. For a very simple example, consider Df to be a leading derivative of a function f in two equations 0 = Df + A, 0 = Df + 2A with A a sum of a large number of terms. A simplification step in the standard procedure would aim at eliminating Df and get as a consequence the system 0 = Df + A, 0 = A where the big expression A occurs twice. An alternative to be described in this paper would be to try to shorten equations and therefore to get at first the system 0 = Df + A, 0 = Df and in a second length reduction step the system 0 = A, 0 = Df.

For a slightly more realistic example, consider two equations 0 = A + C and 0 = B − C, where A, B, C are differential expressions, A, B having only few terms and C involving many terms. If both equations are long and if they involve high derivatives then they would have a low priority to be used in standard algorithms. Typically, each of them would be paired individually with short low order equations for their reduction or the generation of integrability conditions. A simplification of both equations to the system 0 = A + C and 0 = A + B would usually not be found. Instead the expression C would grow in both equations when any substitutions of functions are made that occur in C.

∗This work was partially supported by the Natural Sciences and Engineering Research Council of Canada under CRD Grant # 661-090/93 and partially supported by the Engineering and Physical Sciences Research Council England under Grant # GR/L99036

Both examples are clear cut situations where big expressions in different equations can cancel each other. Although in practical applications to be described in a later section each length reduction step saves only a small number of terms, these reduction steps can often be repeated many times.

The need to reduce the length of equations comes from the danger of dealing with long equations in elimination algorithms. For elimination algorithms to be finite, i.e. to involve only a finite number of steps, they have to eliminate the leading derivative of equations first which involves differentiations and multiplications. Both tend to increase the length of long equations even further. However, elimination algorithms would remain finite if terms other than the leading derivative would be eliminated a finite number of times. Flexibility of this kind could be used to prevent excessive expression swell. In this way not only memory is saved. Long expressions also require an increased time to be computed which slows down any future computations in which long expressions are involved.

The procedure to be described in section 2 is a first step in the direction of an ‘intelligent’ and more efficient computation. It aims at finding equations in the algebraic ideal of the given system with fewer terms. The basic idea is rather straightforward. Any pair of two equations of a given system of equations is checked whether there is a linear combination (with non-vanishing coefficients) of these two equations that is shorter than the longer of the two. If that is the case then the longer one is replaced by the shorter new equation. To find a length reducing linear combination of a pair of two equations, each term of one equation is divided by each term of the other equation, the quotient is simplified, i.e. common factors of numerator and denominator are dropped and a counter of the number of occurrences of this quotient is incremented. The quotient occurring most often is picked, its numerator and its denominator are the multipliers of both equations which are then subtracted from each other. By choosing the quotient that occurs most frequently a maximum number of terms cancel and the result is as short as possible.

The main content of the first part of the paper is to introduce data structures (L, c_{i} below) which allow an efficient implementation and to describe some optimizations that speed up the method and that allow one to consider non-linear equations for length reduction. Beneficial side-effects are discussed such as the increased chance to find ordinary differential equations (ODEs) or to find exact differential equations in a length reduced system of partial differential equations (PDEs).

In the second half, in section 3 the method is applied to determining equations of Killing vectors and Killing tensors in General Relativity. One example is explained in more detail in the appendix where the beneficial side-effects of length reduction become important.

The length reduction module is incorporated in the package Crack for solving over determined PDE-systems. To show that the usefulness of the length reduction module is not just based on special features of Crack but is of a more universal nature a test is described at the end of section 3. In this test other well known programs are used to solve a system of PDEs before and after it has been length reduced. The suitability of length reduction as a pre-processing step is demonstrated.

2 The term reduction method

Introduction The procedure to be described takes as input two expressions E_{1}, E_{2}. Both represent equations 0 = E_{1} = E_{2} and are regarded as sums with n_{1} and n_{2} terms (n_{1}, n_{2} ≥ 1). The aim is to multiply each with a different single term and to add the expressions such that the result has fewer terms than max(n_{1}, n_{2}) and can be substituted for the longer of E_{1}, E_{2}. For a given system of equations this process is repeated with all possible pairs of equations until no pair can produce an equation that is shorter than the longer of both. The restrictions to multiply only with monomials and to combine only two equations at a time are non-trivial constraints (see the discussion in section 4).

Example: The method will be illustrated with the following two expressions

E_{1} = 2xf + 6yf + 4xg + 5x = ∑_{i=1}^{4} E_{1i}  (1)
E_{2} = 3yf − 3xf + 6yg − 7y = ∑_{i=1}^{4} E_{2i}  (2)

f and g are the unknowns that are to be computed from 0 = E_{1}, 0 = E_{2}. To explain the method it does not matter if f and g would be replaced by derivatives of unknown functions or by products of powers of different derivatives as long as f ≠ g. x and y are independent variables or parameters such that 0 = E_{1}, 0 = E_{2} are to be satisfied for any value of x and y.

The treatment of two equations Given are two expressions E_{1}, E_{2} with n_{1} and n_{2} terms, n_{1} ≥ n_{2}. If each expression is multiplied with a single term (monomial) and both expressions are added, then their sum E_{3} can have between n_{1} − n_{2} and n_{1} + n_{2} terms depending on how many terms cancel each other. A way to find the optimal cancellation, i.e. optimal multipliers is to divide each term of expression E_{1} by each term of E_{2} and to collect the simplified quotients (common factors of numerator and denominator dropped) together with the multiplicity they occur.

Example: E_{1}, E_{2} given in (1),(2) have 4 terms each (n_{1} = n_{2} = 4) and the quotients are

(E_{1i}/E_{2j}) = ( \begin{matrix}\frac{2x}{3y} & −\frac{2}{3} & \frac{xf}{3yg} & −\frac{2xf}{7y} \\ 2 & −\frac{2y}{x} & \frac{f}{g} & −\frac{6f}{7} \\ \frac{4xg}{3yf} & −\frac{4g}{3f} & \frac{2x}{3y} & −\frac{4xg}{7y} \\ \frac{5x}{3yf} & −\frac{5}{3f} & \frac{5x}{6yg} & −\frac{5x}{7y} \end{matrix} ) .

A new equation E_{3} would be generated by picking a quotient, say E_{13}/E_{21} = \frac{4xg}{3yf} and using its numerator and denominator to compute E_{3} = 3yf · E_{1} − 4xg · E_{2}. Because this quotient involves f and g, the new equation 0 = E_{3} is non-linear in f, g. As a consequence after replacing E_{1} by E_{3} the new system 0 = E_{3} = E_{2} may not be equivalent to the old system 0 = E_{1} = E_{2}, i.e. E_{1} could not automatically be replaced by E_{3}. The algorithm will therefore not consider quotients that involve any unknowns, here f and g. An effective method to avoid the computation of such quotients will be explained further below.

The quotient occuring most often is \frac{2x}{3y} which is appearing twice and is free of f and g. Each appearance of a quotient means that two terms cancel. We therefore would expect that in the expression E_{3} := 3y · E_{1} − 2x · E_{2} two times two terms cancel and therefore E_{3} has 4 + 4 − 2 · 2 = 4 terms. But E_{3} = 6fx^{2} + 18fy^{2} + 29xy has only 3 terms. When computing E_{3} := 3y · E_{1} − 2x · E_{2} the two terms 3y · E_{14} − 2x · E_{24} = 3y · 5x − 2x · (−7y) = 15xy + 14xy add up to only one term 29xy. We could have forecast the saving of one additional term by realizing that the quotient E_{14}/E_{24} = −\frac{5x}{7y} differs from E_{11}/E_{21} = \frac{2x}{3y} only by a numerical factor.

The example implies that we should record all quotients with the multiplicity they occur and group them into classes c_{i} where all quotients in a class differ by only a numerical factor. Finally, for each class c_{i} the sum M_{i} of all the multiplicities of all quotients in the class is recorded too. All classes c_{i} together with M_{i} are listed in a list L:

L = ((c_{1}, M_{1}), (c_{2}, M_{2}), …, (c_{r}, M_{r}))
c_{i} = ((q_{i1}, m_{i1}), (q_{i2}, m_{i2}), …, (q_{is_{i}}, m_{is_{i}})).

q_{ij} are the different quotients such that two quotients q_{ij}, q_{ik} in the same class c_{i} differ only by a numerical factor.

s_{i} is the number of different quotients in the class c_{i}.

m_{ij} is the number of how often q_{ij} occurs.

M_{i} are defined as M_{i} = ∑_{j=1}^{s_{i}} m_{ij}.

L is the complete (unsorted) list of all classes of quotients.

Disregarding quotients involving f or g in the above example we have

L = ((((\frac{2x}{3y}, 2), (−\frac{5x}{7y}, 1)), 3),

(((−\frac{2}{3}, 1), (2, 1)), 2),

(((−\frac{2y}{x}, 1)), 1)

),

with, for example, q_{11} = \frac{2x}{3y} turning up twice, once as E_{11}/E_{21} and once as E_{13}/E_{23}, therefore m_{11} = 2, and further c_{1} = ((\frac{2x}{3y}, 2), (−\frac{5x}{7y}, 1)), M_{1} = m_{11} + m_{12} = 2 + 1 = 3.

If a quotient q_{ij} is used to combine E_{1}, E_{2} to

E_{3} = denominator(q_{ij}) × E_{1} − numerator(q_{ij}) × E_{2}

then the number n_{3} of terms of E_{3} is

n_{3} = n_{1} + n_{2} − 2 × m_{ij} (due to m_{ij} complete cancellations of 2 terms )
− ∑_{k=1,k≠j}^{s_{i}} m_{ik} (due to savings of one term each time)  (3)
= n_{1} + n_{2} − m_{ij} − M_{i}.

The ∑ m_{ik} in eqn. (3) comes from simplifications like 15xy + 14xy = 29xy which each save one term. In order to be successful and to replace E_{1} by E_{3} we need to find a quotient q_{ij} such that E_{3} has fewer terms than E_{1}, i.e. n_{3} = n_{1} +n_{2} −m_{ij} −M_{i} < n_{1}, hence m_{ij} +M_{i} > n_{2}.

A pre-processing step As argued above only quotients q_{ij} should be considered which do not involve any unknowns, i.e. functions or constants that are to be computed from E_{i} = 0. An effective method to even avoid the computation of those quotients requires to re-write expressions E_{1}, E_{2} in the following way. This initial re-writing step also clarifies how the method works for non-linear expressions E_{1}, E_{2}.

One always can regard non-linear expressions E_{1}, E_{2} as linear homogeneous expressions in some newly defined variables v_{a} which are linear or non-linear constructs of the dependent variables. For example, for independent variables x, y, dependent variables f = const, g = g(x, y) and

E_{1} = 3x + 3 cos(x)fg_{x} − 12xyg + 6xg − yg + sin(g_{y}),  (4)
E_{2} = 1 + 4fg_{x} − 4yg + 2g  (5)

the related system that is homogeneous and linear in v_{i} would be

E_{1} = 3xv_{0} + 3 cos(x)v_{1} − 12xyv_{2} + 6xv_{2} − yv_{2} + v_{3},  (6)
E_{2} = v_{0} + 4v_{1} − 4yv_{2} + 2v_{2}  (7)

with v_{0} = 1, v_{1} = fg_{x}, v_{2} = g, and v_{3} = sin(g_{y}). After this re-writing the method will investigate the system (6),(7) and avoid quotients q_{ij} that involve any v_{l}.

Methods to increase speed The restriction of not multiplying with factors involving dependent variables enables the following major speed up. Instead of investigating the system (4),(5) and computing 6×4 = 24 quotients we investigate the system (6),(7) and compute only quotients between the terms of the coefficients of the same v_{i} in E_{1} and E_{2}. This reduces the number of quotients to 1 × 1 (for v_{0}) +1 × 1 (for v_{1}) +3 × 2 (for v_{2}) +3 × 0 (for v_{3}) = 8 quotients. For large expressions, or more exactly for a high number of different v_{j} the speed up is naturally much higher. If we have r + 1 new dependent variables v_{0}, …, v_{r} and if we denote the number of terms involving v_{j} in E_{i} as n_{ij}, 0 ≤ j ≤ r then instead of computing (∑_{j=0}^{r} n_{1j}) × (∑_{k=0}^{r} n_{2k}) quotients the more efficient method only computes ∑_{j=0}^{r} (n_{1j} × n_{2j}) quotients.

Another way of increasing efficiency is based on knowing the n_{ij} beforehand. For the terms involving a specific v_{j}, an upper bound on the maximal number of cancellations is min(n_{1j}, n_{2j}), saving twice as many terms. This value summed over j = 0 … r (for each v_{j}) gives an upper bound on how many terms can be saved due to cancellations.

This test can be performed without computing any quotients: If the new dependent variables are v_{0}, …, v_{r} and if we have at the beginning

∑_{j=0}^{r} 2 min(n_{1j}, n_{2j}) ≤ n_{2} ( = ∑_{j=0}^{r} n_{2j})

(n_{2} + 1 is the minimum number of terms to be saved to reach a length reduction) then no length reduction is possible.

This test of a necessary criterion is not only possible at the beginning but also during the computation of quotients. We assume that at first all quotients related to v_{0} are computed, then those related to v_{1}, and so on. When the calculation has reached v_{j} and the first w terms in E_{1} that involve v_{j} have been processed, i.e. all quotients between each of them and all terms in E_{2} with v_{j} have already been computed then the following holds. At most min(n_{1j} −w, n_{2j}) more cancellations related to v_{j} are possible and an upper bound B_{j}(w) of the total number of any cancellations still to be found is given as

B_{j}(w) := min(n_{1j} − w, n_{2j}) + ∑_{i=j+1}^{r} min(n_{1i}, n_{2i}).

At this moment any quotient q_{kl} that has a chance to provide a length reduction must satisfy

M_{k} + m_{kl} + 2B_{j}(w) > n_{2}.

All q_{kl} which do not satisfy this condition can be dropped from the list L. For the same reason no new quotients should be added to L as soon as 2B_{j}(w) ≤ n_{2}. If the list L becomes empty at any time during the computation of quotients, then the search can stop. In that case no length reduction is possible. By dropping quotients from the list L, the updating of m_{ij} and M_{i} speeds up.

A speed up not recommended When deleting quotients from L that have no chance to give a length reduction then it is little extra effort to check whether any q_{kl} of the remaining quotients in L already satisfies M_{k} + m_{kl} > n_{2}, and therefore is guaranteed to provide a length reduction. As soon as such a quotient q_{kl} is found the execution could stop. In practical tests it appeared that the negative effects of an early stop dominate. The first length reducing quotient q_{kl} that is found does not have to be the one giving the highest length reduction possible. After using a sub-optimal q_{kl} to compute E_{3} and to substitute E_{1} no further length reductions may be possible, or even if further length reductions were possible, equations tend to have at least intermediately more terms compared with determining always the optimal quotient that gives the highest length reduction. After a suboptimal length reduction, subsequent pairings with other equations would be slower which would result in an overall slow down. It therefore is recommended to complete the computation of all relevant quotients and not to stop early when the first length reducing q_{kl} is found.

Some time tests The following tests can only provide some idea of the running times of the length reduction method. They are measured in a 8 MB Reduce 3.6 session running under Linux on a 133 MHz Pentium PC (Dec. 1998). Equations which have been paired had been generated with the Reduce command RANDPOLY which allows one to specify the number of terms to be generated. In RANDPOLY the randomly generated coefficients may become zero which in that case results in a polynomial with fewer terms than specified. Therefore the number of terms was chosen somewhat larger to be able to drop the surplus terms and get the required size of the polynomial and perform the following statistics.

Performing a length reduction investigation involves no other risk or cost than the computer time that may be lost if no length reduction was possible. Therefore in figure 1 and table 1 a statistics of investigations of each time two equations is shown where no length reduction could be found, i.e. the worst possible result. This will be referred to as unsuccessful pairings in contrast to successful pairings where a length reduction was possible. Both equations are random polynomials of degree up to 5 and with up to 8 variables. Times are obtained by averaging 20 runs. The individual times in these runs differ typically by up to 30%. The results confirm an overall dependence time ∝ (terms of eqn. 1)×(terms of eqn. 2).

sec 50 40 30 20 100  no of terms in 2nd eqn. (x 100) 10 9 8 7 6 5 4 3 2 1 no of terms in 1st eqn. (x 100) 10 9 8 7 6 5 4 3 2 1

Figure 1. Running times for random polynomials of different sizes.

no of terms of 1st eqn. 100 200 300 400 500 600 700 800 900 1000 time in sec (no success) if 2nd eqn. has 10 terms

0.07 0.23 0.45 0.70 0.72 1.01 1.36 1.8 2.2 2.7

time in sec (no success)

if 2nd eqn. has 1000 terms

4.18 6.23 9.35 12.6 15.4 20.3 24.4 30.2 36.4 43.8

Table 1. Timings of unsuccessful length reduction attempts of one equation with varying length and a second equation that has either 10 or 1000 terms.

Test results shown in table 2 are based on pairing equations which are polynomials of 7^{th} degree with each 300 terms but with a varying number of variables. Times are averaged again over 20 runs.

The effect of efficiency improvements as described in the above sub-section to detect the non-existence of length reductions early can be seen clearly from the second row in table 2. As more independent variables occur the number of different quotients q_{kl} increases and the average frequency for each quotient to appear becomes smaller. This in turn rules out many quotients early in the computation and it becomes clear earlier that no quotient will result in a length reduction if that is the case.

Usually length reductions do not happen with random polynomials of that size. In order to measure computing times for pairings when length reductions were possible, pairs of polynomials had been constructed in the following way: a multiple of one random polynomial P_{1} of 300 terms is added to another random polynomial P_{2} and terms in excess of 300 terms are dropped to obtain a polynomial P_{3} with 300 terms. Length reductions between polynomials P_{1} and P_{3} are investigated which produced the third row in table 2. Two trends seem to be present, one lowering the time with an increase of the number of variables (mainly effective between 4 and 5 variables) due to a decrease of potentially successful quotients and another trend slightly increasing the time with the number of variables.

no of variables  4  5  6  7  8
time in msec (unsuccessful) 4720 2461 1673 1656 1640
time in msec (successful)  5706 4968 5137 5388 5728

Table 2: A comparison between average running times of unsuccessful and successful pairings of two equations.

To summarize this sub-section, the main feature found in this study is that computing times are lower in unsuccessful attempts to shorten equations and only when a length reduction becomes possible, i.e. when computing times are of less importance, then they are increased. The speed up measures become increasingly efficient if more unknowns (like f and g in the first example) are present which is the case for partial differential equations with many different partial derivatives acting as different unknowns v_{j} during length reduction.

The order of pairings of equations If more than two equations are given then the question arises in which order they should be paired to search for length reductions. Given that we combine only two equations at a time and multiplying them only with a monomial, we can not expect results that are invariant against combining equations in a different order. The following criteria serve only as a suggestion but they proved to be useful in applications. According to them pairs of equations are picked with the following priorities:

• There should be as few as possible dependent variables v_{i} in the shorter equation which do not occur in the longer equation.

• The shorter of both equations should be as short as possible.

• The longer of both equations should be as short as possible.

The first two rules maximizes the chance to find a reduction of terms. The third rule reduces computation times. The second rule has a higher priority than the third rule because the shorter the equations are, the more useful they are potentially in reducing the length of other equations.

In the following table the above priority list is compared with the same list, only modified by exchanging in the first rule ‘as few as possible’ with ‘as many as possible’. The equations are a set of first order partial differential equations (PDEs) resulting from investigating in General Relativity the Kimura metric [9] with respect to Killing tensors (see section 3). Because sin and cos occur in these differential equations, both length reductions are performed once with the simplification rule cos(x)^{2} ⇒ 1 − sin(x)^{2} and once with the simplification rule sin(x)^{2} ⇒ 1 − cos(x)^{2}. It becomes apparent that these simplifications are not equivalent in their effect.

simplification  original system length red. system length red. system
used using the rule using the rule
‘... as many as ...’ ‘... as few as ...’
no of eqn terms no of eqn terms no of eqn terms
sin(x)^{2} ⇒ 1 − cos(x)^{2}  48  607  25  117  21  81
cos(x)^{2} ⇒ 1 − sin(x)^{2}  48  464  25  108  21  74

Table 3: A comparison of different simplification rules and different rules for the pairing of equations.

As it was to be expected, the ‘... as few as ...’ rule performed better (i.e. resulted in shorter length reduced systems) than the ‘... as many as ...’ rule. What also becomes apparent is that choosing accidentally the less effective simplification rule sin(x)^{2} ⇒ 1 − cos(x)^{2} is less critical when length reduction is performed than without length reduction:

\frac{no of terms of original system using cos-rule}{no of terms of original system using sin-rule} = \frac{607}{464} >

\frac{no of terms of reduced system using cos-rule}{no of terms of reduced system using sin-rule} = \frac{81}{74}.

Beneficial side effects The reduction of length is not only useful for saving memory, and as shorter expressions are quicker to process later on, also for saving time. In this section we want to explain further benefits.

1. Given a set of PDEs, the length reduced system is more likely to contain ordinary differential equations (ODEs) or to contain integrable exact PDEs.

2. Length reduced polynomially non-linear equations are much more likely to be algebraically factorizable.

3. Length reduction of a system of equations has in general the side effect of partially decoupling the system.

If the number of terms of an equation is lowered to, say, n and if the equation is linear then n is necessarily an upper bound for the number of different functions and different derivatives that occur. The length reduction method as described above is indiscriminate to different functions or different derivatives. On average therefore a consequence of a reduction of the number of terms will be a reduction of the number of different functions and the number of different derivatives that occur (see the comparison between tables 6 and 7 below) which in turn provides the above side effects.

About 1: If the number of terms in a length reduced equation got very small (say less than 4) then the chance increases that they involve only one differentiation variable, or that all derivatives can be looked at as derivatives with respect to only one variable of a common partial derivative, like ∂_{y}f, ∂_{xy}f, ∂_{xxy}f are all x-derivatives of ∂_{y}f. In these cases the equation has been reduced to an ODE. (A proper ODE can of course only appear if the differential ideal of the original system does contain ODEs, but that is guaranteed for all the typical sources of over determined PDE-systems, like the computation of infinitesimal symmetries and of conservation laws.)

Similarly the chance increases to obtain an exact differential equation. In order for a differential expression P (f^{i}) that involves functions f^{i} and that satisfies 0 = P to be a total x derivative of some expression I(f^{i}), the identity P = dI/dx has to be satisfied identically in all functions f^{i} and in all their derivatives. The fewer different functions f^{i} and the fewer different derivatives with respect to variables other than x occur, the less restrictive is this assumption of exactness and the more likely it will be satisfied.

The conclusion is not that length reduction is the best way to proceed in order to integrate. We only say that the chance increases to find an integrable PDE in a length reduced system. In the appendix this statement is illustrated by a sequence of 10 integrations which a system of Killing tensor equations admits after being reduced in length.

To explain the usefulness of integrability let us consider an extreme hypothetical example. Gröbner Basis techniques aim at equations with a low differential order. This is a good strategy but not the only way to go. For example, knowing that 0 = ∂^{10}f/∂x^{10} is included in the differential ideal would be very useful as well. After integration, substitution of f and direct separation with respect to x (sometimes called splitting or fragmentation) a highly over determined system for the 10 functions of integration results. Although an expensive Gröbner basis computation might have provided an ODE of lower order than 10, this information is usually gained faster by integrating 0 = ∂^{10}f/∂x^{10} and solving the over determined system that resulted from direct separation for the 10 functions of integration. The key to this speed up would be to sacrifice the minimal differential order for a reduction in the number of terms.

About 2: In a small experiment we want to show that the chance for an algebraic factorization increases with a reduction of the number of terms. Bi-linear polynomials with up to 5 variables and coefficients in the interval −9, …, −1, 1 … 9 have been randomly generated for each length from 2 terms up to 6 terms. The percentage of factorizable polynomials is given in table 4. The chance to have a non-trivial factorization (non-numeric factors) decreases surprisingly quick with an increasing number of terms. Although this test is not proving anything it may serve as an illustration.

no of terms  2  3  4  5  6
factorizable equations in % 55.4 13.6 2.3 0.22 1.3 × 10^{−3}

Table 4: The chances for random quadratic homogeneous polynomials in up to 5 variables to be factorizable in dependence on the number of terms.

About 3: The fewer different functions occur in each equation (on average) the more the system is (at least partially) decoupled and the fewer steps are needed in a subsequent computation to get a differential Gröbner Basis or characteristic system. This means that an elimination algorithm performed afterwards has less work to do, i.e. needs fewer steps. A sparsely occupied system also opens the possibility to choose an appropriate total ordering on which an elimination algorithm is based. For example, if a function turns up in only few equations and with very few different derivatives then this function should get a high priority, i.e. the lexicographical ordering of functions which plays a role in any total ordering should give this function a high priority to be eliminated first because this will probably take the fewest steps to eliminate this function.

3 Applications

One of the applications to which the above length reduction procedure has been applied successfully is the computation of Killing vectors K_{i}(x^{p}) and rank 2 Killing tensors K_{ij}(x^{p}) of given space times with a metric g_{ij} in General Relativity. Killing vectors and Killing tensors provide conservation laws of the form K_{i}u^{i} = const. and K_{ij}u^{i}u^{j} = const. where u^{i} is the 4velocity of geodesic motion in a curved space-time. If enough Killing vectors (i.e. symmetries of the space time) and Killing tensors are found then the equations of free (geodesic) motion of test particles in a curved space can be integrated which is a major step towards the physical interpretation of a space time metric. Killing vectors and their Lie algebras are also widely used to classify space time metrics.

The determining equations for Killing vectors are

K_{i;j} + K_{j;i} = 0, i, j = 1 … 4  (8)

and for Killing tensors they are

K_{ij;l} + K_{jl;i} + K_{li;j} = 0, i, j, l = 1 … 4  (9)

where ‘;’ is the covariant derivative. Conditions (8),(9) are generated for a given space-time with metric tensor g_{ij} by the program Classym described in [13]. Examples are shown in table 5 (not selected from a larger set, but as they appeared in applications). The systems of equations determine either Killing vectors (KV) or Killing tensors (KT) for different space time metrics. One system is extended by integrability conditions (IC) and one system by contractions (CT) with the metric tensor. The number of equations and terms before and after length reduction are shown. Running times in a Reduce 3.6 session running under Linux on a 133 MHz Pentium PC are given in the right column.

equations  original  shortened  no of  time
eqn terms eqn terms red. steps in sec
KV for Taub-NUT (p170 in [5]) + IC 49 1929 44 932  192  120
KT for the Kerr metric (p161 in [5]) 20 1154 20 750  36  87
KT for the Kimura metric [9] + CT  48 464 21 74  93  19.5
KV for the Barnes metric [2]  10 66  10 39  14  2.3
KV for pp waves (p178 in [5])  10 58  10 37  13  1.7

Table 5: The number of length reduction steps and the computation time for the length reduction of different Killing vector (KV) and Killing tensor (KT) determining systems of equations.

The benefits of a length reduction for a program like Crack [14] are manifold. Tables 6 and 7 show the number of terms and the occurring dependent variables in each of the Killing tensor equations for the Kimura metric before and after the reduction of length.

no of  dependent  no of  dependent
terms  variables  terms  variables
2  k01, k00  9  k22, k11, k00, k13, k03, k33, k23
2  k11  9  k22, k11, k01, k00, k02, k03, k33
2  k22, k12  9  k22, k11, k01, k00, k02, k03, k33
3  k13, k33, k23  9  k22, k33, k23
3  k12, k00, k02  9  k11, k01, k00
3  k00, k13, k03  9  k11, k01, k00
3  k11, k01  10  k22, k12, k11, k00, k02, k33, k23
3  k22, k01, k02  11  k22, k12, k13, k33, k23
3  k12, k11  11  k22, k12, k13, k33, k23
3  k11, k13  12  k22, k12, k11, k01, k00, k13, k33
4  k01, k02, k03, k33  14  k22, k12, k11, k01, k13, k02, k03, k33
4  k11, k01, k00  14  k22, k12, k11, k01, k13, k02, k03, k33
4  k12, k01, k02  18  k12, k11, k01, k00, k13, k02, k03
4  k01, k13, k03  18  k12, k11, k01, k00, k13, k02, k03
4  k02, k03, k23  18  k12, k11, k01, k00, k13, k02, k03
4  k22, k12, k11  18  k12, k11, k01, k00, k13, k02, k03
4  k22, k13, k23  19 k22, k12, k11, k00, k13, k02, k03, k33, k23
5  k22, k12, k33, k23  19 k22, k12, k11, k00, k13, k02, k03, k33, k23
5  k12, k11, k13, k33  21 k22, k12, k11, k01, k00, k13, k02, k03, k33
5  k12, k13, k23  21 k22, k12, k11, k01, k00, k13, k02, k03, k33
8  k11, k01, k00  21  k22, k12, k13, k02, k03, k33, k23
8 k12, k11, k00, k13, k02, k03 21  k22, k12, k13, k02, k03, k33, k23
8 k12, k11, k00, k13, k02, k03 21  k22, k12, k13, k02, k03, k33, k23
8  k11, k01, k00  21  k22, k12, k13, k02, k03, k33, k23

Table 6: The original set of conditions for Killing tensors in the Kimura metric.

no of dependent no of  dependent
terms variables  terms  variables
2  k01, k00  4 k22, k13, k33, k23
2  k11  4  k11, k01, k00
2  k22, k12  4 k01, k02, k03, k33
3 k13, k33, k23  4  k12, k01, k02
3 k12, k00, k02  4  k01, k13, k03
3 k00, k13, k03  4  k02, k03, k23
3  k11, k01  4  k22, k12, k11
3 k22, k01, k02  5  k22, k33, k23
3  k12, k11  5 k12, k11, k13, k33
3  k11, k13  5  k12, k13, k23
4 k22, k33, k23

Table 7: The length reduced conditions for Killing tensors in the Kimura metric.

The equations of table 7 (after dropping a single linear dependent equation with 4 terms) are explicitly given in the appendix. It is indicated there how the sparse occurrence of different derivatives leads to integrable ODEs and exact differential equations.

Length reduction as an adjunct to elimination algorithms The main purpose of the term reduction method is to shorten and simplify systems of equations without imposing any risk of an intermediate length increase. This risk is present when using standard Gaussian elimination or standard Gröbner basis methods. For that reason and the benefits discussed above it would not matter if the term reduction method would be comparatively time consuming. The following run time tests show that term reduction can even be time saving if it is used in connections with computing differential Gröbner bases or the radical of a differential ideal.

The program rif [11],[12] has been applied by Allan Wittkopf to the Killing tensor conditions for the Kimura metric (see table 3, the shortened form is given in the appendix). Because we only provide the relative speed up in tables 8,9 the following hard and software specifications are not of much relevance. They are added only for completeness. Computing times of the program rif are measured on a PII 333 MHz Linux system with 384 MB RAM (although only about 3 MB was needed), running in XMaple release 5.0.

Differential Gröbner Bases have been computed, once in the generic case where the constant b is treated like a dependent variable and a second time where b is regarded as an arbitrary constant, ignoring b pivots.

In the tables 8 and 9 kimu1 is the form of the system with 607 terms (see table 3), kimu2 is a pre-optimal form of the system with 128 terms which resulted when a non-optimal pairing of equations was used and finally kimu3 is the system with 74 terms (see appendix).

equations  rif, with b pivots rif, without b pivots
tkimu2 /tkimu1  96.8 %  93.2 %
tkimu3 /tkimu1  86.6 %  72.3 %

Table 8: Relative speed up of computing times (cpu time) for the program rif for different versions of the Killing tensor conditions of the Kimura metric.

The Maple program diffalg has been applied by Evelyne Hubert (see [3, 4, 7, 8] and

the URL http://daisy.uwaterloo.ca/~ehubert/Diffalg). diffalg was run on a PC with

dual Intel 400Mhz CPU’s, 512 MB RAM, 128 MB swap, Asus P2L7-DS motherboard, 4.3 GB Seagate hard drive and Intel EtherExpress Pro 10/100 network card under LINUX.

equations  diffalg
tkimu2 /tkimu1 83.7 %
tkimu3 /tkimu1 40.1 %

Table 9: Relative speed up of computing times (cpu time) for the program diffalg for different versions of the Killing tensor conditions of the Kimura metric.

A potential gain of a reduction of the number of terms is the possibility to choose a more appropriate total ordering in a Gröbner Basis computations when a system is already partially decoupled. For example, dependent variables which turn up in fewer equations could be given a higher lexicographical priority than other dependent variables. So far no advantage has been taken of this opportunity in the package Crack.

4 Possible extensions

The manipulations studied in this paper to achieve a reduction of length seem very special and one might seek a more comprehensive theory, for example, covering all possible linear combinations of any number of equations multiplied not only with monomials but with any number of terms. This problem is harder than it looks. It would include being able to decide whether any given linear algebraic system could be solved without intermediate memory increase. Already the question of combining 3 equations to kill at least 5 terms (4 terms could be killed by combining 2 equations twice) has a much larger search space than combining 2 equations. Although not proven, the author expects that in generic practical applications the frequency of reducing the length of any one of 3 equations which can not be reached by combining repeatedly 2 equations is low (which, of course, may be different in specific applications). To give an example, the system 0 = a+b−c−d+g, 0 = c+d−e−f +h, 0 = e + f − a − b + k gives a shortened equation only if all three equations are combined (added).

It seems to be a general phenomenon that the computational complexity increases drastically when any method of investigating systems of equations is generalized whereas the success rate of the generalized method increases only marginally.

To give an example, in determining Lie-symmetries and conservation laws of generic partial differential equations (not the rare fully integrable PDEs) the success rate in finding symmetries when going from point symmetries to higher order symmetries or the success rate in finding conservation laws when going from zeroth order integrating factors to higher order conservation laws increases only marginally. On the other hand the complexity in solving the related over determined system of conditions grows drastically when the order of the ansatz is increased. The growth in complexity comes from the task to determine an increased number of constants or functions or in determining functions of more independent variables. The decreasing success rate results from a relatively more over determined problem: having to satisfy not only the more restrictive conditions of a more general ansatz but also the conditions that the solution may not be decomposable into simpler cases.

5 Summary

Although conceptually simple, the method explained in this paper proved to be very useful. It is fast and it has zero risk of increasing the size of the system of equations during the computation or as a result of it. This is an important feature for very memory intensive computations. In such a case, computation times are of less concern. Still, it was possible to show that even for other well known programs the time saving using term reduction could be higher than the additional cost.

As discussed in the subsection ‘Other benefits’, size reduced equations on average are more likely to be an ODE or to be a total derivative and therefore easily integrable. They are also more likely to be algebraically factorizable. The fact that the system becomes partially decoupled opens the possibility to select a total ordering according to which the system is already close to a characteristic form. This has the potential to reduce subsequent computations to obtain a characteristic system to a high extend.

6 Acknowledgement

The author thanks Allan Wittkopf for many discussions. He, together with Evelyne Hubert are thanked for running the tests mentioned in the text with their systems. The Symbolic Computation Group at the University of Waterloo is thanked for its hospitality during the authors sabbatical in the fall of 1998 when the manuscript was prepared.

7 Appendix: Length reduced conditions for Killing tensors of the Kimura metric

The size reduced determining conditions for Killing tensors in the Kimura metric read:

0 = k_{22,h} + 2b^{2}r^{3}k_{12}  (10)
0 = rk_{11,r} + 2k_{11}  (11)
0 = k_{00,t} − 2br^{3}k_{01}  (12)
0 = 2 cos(h) sin(h)k_{23} + k_{33,p} + 2 sin(h)^{2}b^{2}r^{3}k_{13}  (13)
0 = rk_{11,p} + 2rk_{13,r} − 2k_{13}  (14)
0 = rk_{11,h} + 2rk_{12,r} − 2k_{12}  (15)
0 = 2k_{02,h} + k_{22,t} + 2b^{2}r^{3}k_{01}  (16)
0 = 2rk_{01,r} + rk_{11,t} − 2k_{01}  (17)
0 = k_{00,p} + 2k_{03,t} − 2br^{3}k_{13}  (18)
0 = k_{00,h} + 2k_{02,t} − 2br^{3}k_{12}  (19)
0 = 6 cos(h) sin(h)k_{23} − sin(h)^{2}k_{22,p} − 2 sin(h)^{2}k_{23,h} + k_{33,p}  (20)
0 = 2 cos(h) sin(h)k_{02} + 2k_{03,p} + k_{33,t} + 2 sin(h)^{2}b^{2}r^{3}k_{01}  (21)
0 = 2rk_{12,h} + rk_{22,r} + 2b^{2}r^{4}k_{11} − 4k_{22}  (22)
0 = 2 cos(h)k_{03} − sin(h)k_{02,p} − sin(h)k_{03,h} − sin(h)k_{23,t}  (23)
0 = rk_{01,p} + rk_{03,r} + rk_{13,t} − 4k_{03}  (24)
0 = rk_{01,h} + rk_{02,r} + rk_{12,t} − 4k_{02}  (25)
0 = rk_{00,r} + 2rk_{01,t} + br^{5}k_{11,r} − 4k_{00}  (26)
0 = 2 cos(h) sin(h)^{2}k_{22} − 4 cos(h)k_{33} − sin(h)^{3}k_{22,h} + 2 sin(h)k_{23,p}
+ sin(h)k_{33,h}  (27)
0 = 2 cos(h) sin(h)rk_{12} + 2rk_{13,p} + rk_{33,r} + 2 sin(h)^{2}b^{2}r^{4}k_{11} − 4k_{33}  (28)
0 = 2 cos(h)rk_{13} − sin(h)rk_{12,p} − sin(h)rk_{13,h} − sin(h)rk_{23,r} + 4 sin(h)k_{23}  (29)

One effect of reducing the number of terms is that fewer different derivatives of different dependent variables k_{ij} occur in each equation. As a consequence the chance for the equations to be exact or to be simple ODEs increases. The package Crack is able to integrate 10 of the above equations and as a consequence to express 10 of the unknown functions k_{ij} in terms of new unknown functions of integration of only 3 variables: integration of (11) to solve for k_{11}, of (14) to solve for k_{13}, of (15) to solve for k_{12}, of (10) to solve for k_{22}, of (17) to solve for k_{01}, of (12) to solve for k_{00}, of (18) to solve for k_{03}, of (16) to solve for k_{02}, of (28) to solve for k_{33} and of (29) to solve for k_{23}. These integrations and substitutions are a first step in a longer calculation which finally gives the following general solution. The 13 free constants c_{1} … c_{13} stand for 13 Killing tensors, i.e. 13 conserved first integrals of geodesic motion in the curved space described by the Kimura metric.

k_{00} = (3b^{2}r^{4}t^{2}c_{1} + 48b^{2}r^{4}tc_{2} + 3br^{2}c_{1} − br^{4}c_{3} − 3br^{4}t^{2}c_{4} − r^{2}c_{4}) /(6b)
k_{01} = (brtc_{1} + 8brc_{2} − rtc_{4}) /(2b)
k_{02} = cos(p)r^{4}c_{12} + sin(p)r^{4}c_{11}
k_{03} = (− cos(h)^{3} cos(p)r^{4}c_{11} + cos(h)^{3} sin(p)r^{4}c_{12} − cos(h)^{2} sin(h)r^{4}c_{13}
+ cos(h) cos(p)r^{4}c_{11} − cos(h) sin(p)r^{4}c_{12} + sin(h)r^{4}c_{13}) / sin(h)
k_{11} = −c_{4}/(3b^{2}r^{2})
k_{12} = 0
k_{13} = 0
k_{22} = (6 cos(p)^{2}r^{4}c_{6} + 6 cos(p) sin(p)r^{4}c_{5} − 3b^{2}r^{4}c_{1}t^{2} − 48b^{2}r^{4}tc_{2} + 3br^{4}t^{2}c_{4}
−6r^{4}c_{7} − 2r^{2}c_{4}) /6
k_{23} = (2 cos(h)^{2} cos(p)r^{4}c_{8} − 2 cos(h)^{2} sin(p)r^{4}c_{9} + 2 cos(h) cos(p)^{2} sin(h)r^{4}c_{5}
−2 cos(h) cos(p) sin(h) sin(p)r^{4}c_{6} − cos(h) sin(h)r^{4}c_{5} − 2 cos(p)r^{4}c_{8}
+2 sin(p)r^{4}c_{9}) /2
k_{33} = (12 cos(h)^{5} cos(p)r^{4}c_{9} + 12 cos(h)^{5} sin(p)r^{4}c_{8} + 6 cos(h)^{4} cos(p)^{2} sin(h)r^{4}c_{6}

+6 cos(h)^{4} cos(p) sin(h) sin(p)r^{4}c_{5} − 6 cos(h)^{4} sin(h)r^{4}c_{7} + 6 cos(h)^{4} sin(h)r^{4}c_{10}

−24 cos(h)^{3} cos(p)r^{4}c_{9} − 24 cos(h)^{3} sin(p)r^{4}c_{8} − 6 cos(h)^{2} cos(p)^{2} sin(h)r^{4}c_{6}
−6 cos(h)^{2} cos(p) sin(h) sin(p)r^{4}c_{5} + 3 cos(h)^{2} sin(h)b^{2}r^{4}t^{2}c_{1} − 12 sin(h)r^{4}c_{7}
+48 cos(h)^{2} sin(h)b^{2}r^{4}tc_{2} − 3 cos(h)^{2} sin(h)br^{4}t^{2}c_{4} − 6 cos(h)^{2} sin(h)r^{4}c_{6}
+18 cos(h)^{2} sin(h)r^{4}c_{7} − 12 cos(h)^{2} sin(h)r^{4}c_{10} + 2 cos(h)^{2} sin(h)r^{2}c_{4}
+12 cos(h) cos(p)r^{4}c_{9} + 12 cos(h) sin(p)r^{4}c_{8} − 3 sin(h)b^{2}r^{4}t^{2}c_{1} + 6 sin(h)r^{4}c_{10}
−48 sin(h)b^{2}r^{4}tc_{2} + 3 sin(h)br^{4}t^{2}c_{4} + 6 sin(h)r^{4}c_{6} − 2 sin(h)r^{2}c_{4}) /(6 sin(h))

What makes this result and the Kimura metric interesting is that two of these Killing tensors are non-trivial, i.e. they are not just symmetrized products of the 4 Killing vectors of this metric (see, for example [9]).

References

[1] Becker, T., Weispfennig, V. Gröbner Bases, Springer-Verlag, Graduate Texts in Mathematics, 141, (1993) p 222 ff.

[2] Barnes, A., On space-times admitting a three-parameter isometry group with twodimensional null orbits, J. Phys. A: Math. Gen. 12, no 9 (1979), 1493-97.

[3] Boulier, F., Lazard, D., Ollivier, F., Petitot, M. Computing representations for radicals of finitely generated differential ideals. Proceedings of ISSAC’95, ACM Press (1995), 158-166.

[4] Boulier, F., Lazard, D., Ollivier, F., Petitot, M. Computing representations for radicals of finitely generated differential ideals. Technical Report IT-306, LIFL. (1997).

[5] Hawking, S. W., Ellis, G. F. R., The large scale structure of space-time, Cambridge Univ. Press (1973).

[6] Hereman, W., Symbolic Software for Lie Symmetry Analysis. Chapter 13 in CRC Handbook of Lie Group Analysis of Differential Equations, vol. 3, ed. NH Ibragimov. Boca Raton, Florida: CRC Press, 1996, pp. 367–413

[7] Hubert, E., Essential Components of algebraic differential equations, J. of Symb. Comp., 28 (1999) no 4-5, pp 657-680.

[8] Hubert, E., Factorisation free decomposition algorithms in differential algebra, J. of Symb. Comp., 29 (2000) no 4-5.

[9] Kimura, M., On quadratic first integrals in static spherically symmetric space-times, having spacial parts of non-constant curvature, I., Tensor N.S., 30 (1976), 27-43.

[10] Mansfield, E.L., The differential algebra package diffgrob2, Mapletech 3, (1996) 33-37.

[11] Reid, G.J., Wittkopf, A.D. and Boulton, A., Reduction of systems of nonlinear partial differential equations to simplified involutive forms, European Journal of Applied Mathematics, Vol 7. (1996) 604-635.

[12] Reid, G.J., Lin, P., and Wittkopf, A.D., Differential-Elimination Completion Algorithms for DAE and PDAE, Studies in Applied Mathematics 106(1) (2001) 1-45.

[13] Wolf, T., Grebot, G., Automatic symmetry investigation of space-time metrics, Int. J. of Mod. Phys. D, 3, (1994) 323-326.

[14] Wolf, T., The program CRACK for solving PDEs in General Relativity, in F.W. Hehl, R.A. Puntigam, and H. Ruder (Eds.): Relativity and Scientific Computing: Computer Algebra, Numerics, Visualization Springer Verlag, (1996), 241-258.