The Symbolic Integration of Exact PDEs

Thomas Wolf School of Mathematical Sciences Queen Mary and Westfield College Mile End Road London E1 4NS England September 8, 2000

Abstract An algorithm is described which decides if a given polynomial differential expression Δ of multivariate functions is exact, i.e. whether there exists a first integral P such that D_{x}P = Δ for any one x of a set of n variables and to provide the integral P. A generalization is given to allow integration in the case that the exactness is prevented by terms which contain only functions of less than n independent variables.

1 Motivation

The common way to deal with problems that involve the solution of non-linear differential equations is to try different ansätze which are either geometrically motivated or just chosen to simplify computations. Typical examples are the investigation of infinitesimal symmetries, the search for classes of integrating factors and related first integrals/conservation laws or the search for a variational principle equivalent to a given system of equations. In all these cases over-determined systems of partial differential equations (PDEs) have to be solved. If their solution could be automated then not only would erroneous hand calculations be speeded up but also less over-determined PDE-systems could be solved, i.e. less restrictive ansätze be investigated. One commonly used technique is to apply integrability conditions systematically by computing a characteristic system of the differential ideal [4] or a differential Gröbner basis (or a Pseudo differential Gröbner basis) [5, 6, 7]. Other related programs are listed in [3].

The algorithms to be described in this paper aim to extend such algorithms by enabling the integration of classes of equations either in parallel to the computation of integrability conditions or after a basis for the differential ideal/differential Gröbner basis has been computed. It appears to the author that the problem of adding integrations splits naturally into two parts answering the following questions.

1. What is a good trade off between the generality of the class of equations to be integrated and the practicality of the corresponding integration algorithm? In other words, for which class of equations is an effective and finite algorithm available? At which computational price can such an algorithm be generalized? Are there fast necessary tests available for the integration algorithm to be chosen?

2. How can any integration method of some class of equations be incorporated into the Gröbner basis calculation? What are the conditions which keep the extended Gröbner basis computation finite? With which priority should integrations be performed within such a combined machinery? Which one of a set of integrable equations should be integrated first and if an equation can be integrated, should it be integrated as often as possible? Could the Gröbner basis algorithm be modified to cover the same functionality without having to perform integrations explicitly? Which other techniques, like separation (splitting/fragmenting an equation) become necessary, once integrations are performed? Should any equation which could be integrated also be used for substitutions of unknown functions? Can one give general rules for the order in which the different techniques (computing integrability conditions, integrations, separations) should be applied? Do integrations which integrate an equation if possible at least once and integrations - which integrate an equation only if it can be integrated often enough to enable a substitution - have different priorities? Which data structures should be used to represent differential equations in order to organize the application of different modules efficiently?

In this paper we are going to answer the first set of questions. To deal with the second set of questions properly one has to discuss no less than a complete package of modules aiming at the explicit solution of over-determined PDE-systems. This will be done elsewhere [10].

In the following section we will outline alternative integration problems before in the subsequent section the algorithm is explained. In the concluding sections extensions are described.

2 Choosing an integration task

There is a wide variety of integration related tasks from which we choose one that is to be addressed by the algorithm described in this paper. Our aim will be to choose a type of problem that is as general as possible but that can still be answered algorithmically and efficiently.

On the complex end of possible problems we have the computation of conservation laws, i.e. relations

Div P = D_{x^{1}}P^{1} + … + D_{x^{p}}P^{p}  (1)
= ∑_{ν,J} Q_{ν}^{J}D_{J}Δ_{ν} (identically in all x^{i}, u_{J}^{α})  (2)

where x^{1}, …, x^{p} are the independent variables, u^{1}, …, u^{q} are the dependent variables, 0 = Δ_{ν}, ν = 1, …, r are the given differential equations, D_{x^{i}}, D_{J} are total differential operators with _{J} being a multi-index, standing, for example, for _{1223}, i.e. standing for a derivative with respect to x^{1}, x^{2}, x^{2}, x^{3}. The unknowns in equation (2) that are to be determined are the components of a conserved current P^{i} and coefficients Q_{ν}^{J}, all being expressions in x^{i}, u_{J}^{α}. We will keep this notation for the rest of the paper.

Partial integration of the right hand side of relation (2) and transferring divergences to the left hand side can always reach a form

Div P̃ = ∑_{ν} Q̃_{ν}Δ_{ν}.  (3)

A conservation law could be used to express P by the curl of new potentials that are introduced. In two recent papers [12, 13] it is shown that there are very good chances of solving the resulting over-determined systems by computer algebra but the programs have to perform integrations as part of their solution which is what we want to do in this paper. What one would gain from such an involved computation is an invariant statement. The fact whether or not a divergence Div P is contained in the differential ideal of a given system of equations is not altered by a Gröbner basis computation which does not change the differential ideal. (What may change during a Gröbner basis computation is the order of u-derivatives that the Q_{ν}^{J} involve.)

In order to perform integrations without having to solve a system of over-determined equations we will make more restrictive assumptions in this paper. We assume all but one P^{i} to vanish. Restricting ourselves to the integration of a single equation 0 = Δ we assume the single integrating factor Q to be 1. Further we assume the equations 0 = Δ_{ν} to be at most polynomially non-linear in the dependent variables u and their derivatives. A few minor generalizations will be added later.

Even with this specification it appears that we are more general than most other programs that are listed in [3] which aim at solving over-determined systems. Two strong programs to be mentioned are Dimsym of James Sherring, Geoff Prince and Michael Jerie and Lie of Alan Head. The integration techniques used in these programs are partially described in [2], [8] and [9] but to the best of the author’s knowledge not in the published literature.

Other programs either do not integrate or their integration capabilities are restricted to PDEs of the form

0 = ∂_{x^{1}}u(x^{i}) − g(x^{j})

where u is an unknown function and g is a differential expression which involves x^{1} only explicitly and polynomially and not as an argument of another unknown function, or their integration algorithm is a translation (see [1]) from the one in the package Crack with contributions from the package Lie. The implementation in Crack is based on the algorithm to be described in this paper. The implementation in [1] does not contain the generalization given in section 4.

The basic form of the algorithms discussed here will be to determine for a given equation 0 = Δ any independent variable x^{i} and differential expression P (x^{i}, u^{α}) which is functionally dependent only on u^{α} and their partial derivatives such that

D_{x^{i}}P = Δ.

In doing that we lose invariance with respect to steps of a Gröbner Basis computation, as is seen in the following short example. The equation ∂_{x}^{2}u = 0 satisfies our integration criterium but adding the two equations 0 = ∂_{x}∂_{y}u− \frac{1}{x}∂_{y}u, 0 = ∂_{y}^{2}u−u and computing the Gröbner basis of all three as 0 = x∂_{x}u − u, 0 = ∂_{y}^{2}u − u we lose this property.

Although the algorithm to be described is not the most general possible, we will see that it is fast and still frequently applicable. The faster an integration algorithm is, the lower is the risk of wasting computational resources in unsuccessful tries.

3 The algorithm for exact DEs

A skeleton program In the following we consider integration with respect to only one variable x. In the presence of different variables the computation would have to be repeated for each of them. An algorithm for adding constants/functions of integration will be given further below.

The expression Δ which is to be integrated is assumed to be a polynomial in u^{α} and their derivatives. In the following algorithm % denotes the beginning of a comment ending at the end of a line, highpow(A, B) denotes the highest power of B in the expression A, coeff(A, B) denotes the coefficient of B in A and LD(Δ, u^{α}) denotes the leading derivative of the function u^{α} in the expression Δ with respect to some lexicographic total ordering >_{lex} of variables in which the integration variable x has highest priority. We further denote by u_{J−x}^{α} a derivative of u^{α} that has one less x− derivative than u_{J}^{α}. The lower multi-index _{K} is reserved to annotate derivatives u_{K}^{α} which do not contain derivatives with respect to the integration variable x.

1 Algorithm EXACT 2 Input Δ, x % Δ is a differential expression, x is an independent variable 3 Output P % with D_{x}P = Δ, if Δ is exact 4 nil % if Δ is not exact 5 Body 6 P := 0 7 while Δ contains an unknown function u^{α}(x) or its derivatives do 8 u_{J}^{α} := LD(Δ, u^{α}) 9 while u_{J}^{α} involves x-derivatives do 10 if highpow(Δ, u_{J}^{α}) ≠ 1 then return nil % Δ is not exact 11 a := coeff(Δ, u_{J}^{α}) 12 b := ∫ a d(u_{J−x}^{α}) 13 P := P + b 14 Δ := Δ − D_{x}b

15 uα ; J̃ := LD(Δ, u^{α})

16 if uα ; J̃ >_{lex} u_{J}^{α} then return nil % Δ is not exact

17 u_{J}^{α} := uα ; J̃

18  if Δ involves u^{α} then return nil  % Δ is not exact
19  P := P + ∫ Δ dx  % by now Δ contains x only explicitly

20 return P

Scope of applicability The algorithm is applicable for expressions Δ that are polynomially non-linear in the u^{α} but also for non-polynomially non-linear Δ provided the integral in line 12 can be computed and the simplification of Δ in line 14 is effective. This means, the expression for Δ computed in line 14 must be free of u_{J}.

If Δ involves terms that do not involve unknowns u^{α}(x) but are completely explicit in x then in line 19 an integral of them is added to P. Whether that integral can be expressed in terms of elementary functions or not does depend on these terms but this “integrability” is irrelevant for the exactness of Δ.

The overall structure The overall structure of the algorithm is to integrate successively highest x-derivatives of any function of x. The highest x− derivatives must occur linearly which is tested in line 10. After integrating all terms that include an x-derivative of some function u^{α}(x), this function u^{α} and any derivatives of u^{α} other than x-derivatives must not occur any more in the so far unintegrated terms Δ (line 18). Otherwise the remaining terms Δ and consequently the original Δ are not exact with respect to x. The reason for this is, an x-integral of the original Δ would have to involve an x-integral of u^{α} or an x-integral of (non x-) derivatives of u^{α} which is not what we want as an integral for Δ.

Example: For g = g(x), ′ = d/dx, Δ = gg^{′3} + xg^{′4} + 3xgg^{′2}g′′ the computation is completed by going through the inner loop (steps 8-17) only once:

u_{J}^{α} := g′′
highpow(Δ, g′′) = 1
a := 3xgg^{′2}
u_{J−x}^{α} = g′
b := ∫ 3xgg^{′2} d(g′) = xgg^{′3}
P := xgg^{′3}
Δ := Δ − (xgg^{′3})′ = 0.

The integral is P := xgg^{′3}. The generation of constants of integration is described further below.

The non-integrability tests The test in line 16 is motivated with the following simple example. We assume, the single term Δ = u_{xxy}u_{xyy} with u = u(x, y) is to be x-integrated. The lines 11-14 would perform the partial integration ∫ u_{xxy}u_{xyy} dx = u_{xy}u_{xyy} − ∫ u_{xy}u_{xxyy} dx. Going through these lines the second time to compute ∫ u_{xy}u_{xxyy} dx would reverse this partial integration and so on. The purpose of the test in line 16 is to prevent an infinite loop.

In this example the steps are:

u_{J}^{α} := u_{xxy}
highpow(Δ, u_{xxy}) = 1
a := u_{xyy}
u_{J−x}^{α} = u_{xy}
b := u_{xyy} d(u_{xy}) = u_{xyy}u_{xy}
P := u_{xyy}u_{xy}
Δ := Δ − (u_{xyy}u_{xy})_{x} = −u_{xxyy}u_{xy}
uα ; J̃ := u_{xxyy}
u_{xxyy} >_{lex} u_{xxy} ⇒ Δ is not exact.

Proof for the test in line 16 to be a necessary integrability criterion:

If an expression Δ has an integral P with D_{x}P = Δ then for a leading derivative LD based on a lexicographical ordering with x having highest priority we have

D_{x}(LD(P, u^{α})) = LD(D_{x}P, u^{α}) = LD(Δ, u^{α}) = u_{J}^{α}.  (4)

This argument justifies the test in line 16 to be a necessary criterium for integrability. The case uα ; J̃ = u_{J}^{α} can theoretically not occur. (It can occur in practice if expressions Δ are input that are not polynomial in the u^{α} so that simplification routines in the computer algebra system used are not effective for expressions occurring in Δ and consequently in line 14 the terms with u_{J}^{α} do not cancel. In that case it is appropriate to stop in line 14 too but without concluding non-integrability.)

Finally, if a test in lines 10,16,18 gives true, i.e. if the current value of Δ is not a total x-derivative then the original input Δ cannot be a total x-derivative because the two differ by total derivatives (the repeatedly subtracted D_{x}b in line 14).

The termination of the algorithm The algorithm as shown above is finite. Proof:

1. Assume the leading derivative u_{J}^{α} as determined in line 8 (or determined in line 15,17) is u_{J}^{α} = ∂_{x}^{n}u_{K_{i}}^{α} where n is the highest x-derivative of u^{α} in Δ and u_{K_{i}}^{α} does not involve x-derivatives. Then the test in line 16 guarantees that the new leading derivative uα ; J̃ = ∂ñ ; xuα ; K̃_{i} determined in line 15 has either

(a) ñ < n, or

(b) ñ = n, uα ; K̃_{i} <_{lex} u_{K_{i}}^{α}.

2. There are only finitely many different derivatives u_{K_{m}}^{α} (not involving x derivatives) which either turn up themselves in Δ or x-derivatives of them turn up in Δ. Their number does not increase as only x-derivatives are taken in the algorithm. Therefore case (1b) can happen only a finite number of times before in case (1a) the order in x must be lowered. This in turn can also happen only a finite number of times before u^{α} appears without x-derivatives (end of inner loop) or does not occur anymore in Δ (end of outer loop).

Methods of speeding up the algorithm In this paragraph methods are described which proved useful in the author’s Reduce implementation of the algorithm.

1. The proof given above for the test in line 16 to be a necessary criterion also extends to the test

if LD(a, u^{α}) >_{lex} u_{J−x}^{α} then return nil % Δ is not exact

put between lines 11 and 12. It allows us to drop the test in line 16. If this new test shows LD(a, u^{α}) >_{lex} u_{J−x}^{α} then x-differentiation would yield a contradiction to condition (4) and therefore non-integrability. In the previous example Δ = u_{xxy}u_{xyy} the test would give:

u_{J}^{α} := u_{xxy}
highpow(Δ, u_{xxy}) = 1
a := u_{xyy}
LD(a, u) = u_{xyy}
u_{J−x}^{α} = u_{xy}
u_{xyy} >_{lex} u_{xy} ⇒ Δ is not exact.

If this extra test shows non-integrability then it saves the computationally more expensive steps in lines 13, 14 and 15. If it does not find a contradiction then the effort in determining LD(a, u^{α}) was wasted because LD(Δ, u^{α}) has to be determined in line 15 anyway as it is needed in line 9 when going through the inner loop the next time. This extra but quick test pays off especially if expressions Δ are big and are non-integrable.

2. According to line 7 in the algorithm EXACT at first all terms containing some function u^{α} are integrated, then all terms containing some other functions u^{β} and so on. This requirement can be relaxed as long as it is guaranteed that later partial integrations do not reverse earlier integration steps. One way to ensure this would be to determine in line 8 not the leading derivative LD(Δ, u^{α}) of some given u^{α} but the leading derivative LD(Δ) with respect to any u^{α}. Then the outer while loop in line 7 can be dropped. The effect would be to integrate more x-derivatives of different functions u^{α} before the test in line 18 is carried out with a chance to find non-integrability. Therefore the overall effect would be a slowdown. To reach a speed up the opposite measure could be applied. Instead of partially integrating x-derivatives of any function, only x-derivatives of one special non-x-derivative u_{K_{i}}^{α} of one function are integrated in the inner loop. This can be accomplished, for example, by a preprocessor step in which all different non-x-derivatives u_{K_{i}}^{α} of each function u^{α} would be renamed as different functions U^{β} := u_{K_{i}}^{α} so that Δ involves only x-derivatives of an enlarged set of functions U^{β}. For example, the expression Δ = xu_{xxy} + u_{xy} + u_{xy}^{2} + x^{2}u_{xyy} would be rewritten as Δ = xU_{xx}^{1} + U_{x}^{1} + (U_{x}^{1})^{2} + x^{2}U_{x}^{2} by substituting U^{1} := u_{y}, U^{2} := u_{yy}. This renaming does not affect integrability because only x-differentiations and x-integrations are performed and non-x-derivatives therefore do not change. As a result the test in line 18 is already made after all x-derivatives of one U^{β} have been integrated. Non-integrability is therefore recognized earlier.

3. The above measure can be refined further for the case that Δ is not linear but polynomially non-linear in the U^{β} and their x-derivatives. Differentiations with respect to x do not change the power of any U^{α} in a monomial. This allows us to partition Δ even further into smaller partial sums which all have to be exact if Δ is to be exact. This more fine-grained partitioning is accomplished by associating with each monomial M in Δ a monomial M̃ obtained by dropping all x-derivatives from all functions U^{β} in M but keeping the U^{β} themselves and replacing any U^{β}- independent coefficients by 1. Any two monomials in Δ with the same M̃ belong to the same partial sum. For example, if P contains the monomial xU_{x}^{1}(U^{2})^{2} then all monomials in its x-derivative U_{x}^{1}(U^{2})^{2} +xU_{xx}^{1}(U^{2})^{2} +2xU_{x}^{1}U^{2}U_{x}^{2} are associated to the same M̃ = U^{1}(U^{2})^{2} as is the monomial xU_{x}^{1}(U^{2})^{2} in P. For the example from the previous speed up method this would mean that Δ splits not into two but three parts:

Δ = xU_{xx}^{1} + U_{x}^{1}  (=: first partial sum s_{1} with M̃ = U^{1})
+(U_{x}^{1})^{2}  (=: second partial sum s_{2} with M̃ = (U^{1})^{2})
+x^{2}U_{x}^{2}  (=: third partial sum s_{3} with M̃ = U^{2})

The speed up has the following sources.

(a) By splitting up Δ into smaller partial sums, smaller expressions have to be handled within the inner loop.

(b) Instead of having to integrate all x-derivatives of U^{1} in s_{1} +s_{2} before detecting non-integrability, this would be found quicker if s_{2} is tried before s_{1}. If s_{1} is tried before s_{2} then no significant difference from trying both at once would occur. So on average a speed up results.

(c) The partial sums s_{i} could be sorted by their number of terms before trying their integrability beginning with the shortest s_{i}.

One more remark concerning this efficiency improvement. The overall gain in speed does depend on the overhead spent on partitioning. The easy and quick access to low level data in the symbolic mode of the computer algebra system Reduce makes the partitioning effort low compared with the integration itself. This is different in Maple, where the strength lies in handling relatively efficiently the very large expressions for which the partitioning method is supposed to be especially useful. For these two reasons the partitioning method turned out to be more useful in a Reduce implementation than in a Maple implementation.

Constants and functions of integration To complete an integration a constant (or x-independent function if x is not the only independent variable) has to be added to P. If Δ is to be integrated more than once then it is better to add all constants/functions after all integrations have been completed. Assuming Δ is integrated k_{i}–times with respect to each x_{i}, then we have to add to the integral P an expression P_{c} as computed in the following procedure.

Procedure INTCONST

Input r, m, k_{1}, …, k_{m} % r is the next index to be used to create new names

% c_{r}, c_{r+1} … for constants/functions of integration

Output P_{c}, r

I_{c} := 0

for i := 1 … m do s := 0

for j := 1 … k_{i} do
r := r + 1
s := s · x_{i} + c_{r}(x_{0}, .., x_{i−1}, x_{i+1}, .., x_{m})
P_{c} := P_{c} + s

return P_{c}, r If P is a linear algebraic equation for one of the u^{α} and one intends to eliminate and substitute this u^{α} in other equations, then the choice of the c_{r}’s can be refined. By replacing the original c_{r} by a product of a new c̃_{r} with all those nonzero factors of coeff(P, u^{α}) which only depend on variables of c_{r}, these factors will vanish after substitution. For example, for u = u(x, y) the x-integration of

0 = yu_{x} + 2x

to

0 = yu + x^{2} + yc̃(y)

provides simpler expressions after substitution of u than

0 = yu + x^{2} + c(y).

4 A generalization

The generalization of the above integration algorithm to be discussed in this section is not mathematically as clearly defined as the integration of exact differential expressions. Nevertheless this generalization proves to be very useful when the complete integration module is used as part of a program to solve (usually over-determined) systems of differential equations.

If integrations and substitutions are performed in a system of equations then gradually more and more functions occur that depend on few variables. The following generalization concerns equations which are exact up to terms that contain only unknown functions of fewer variables than occur in that equation. The idea is to reach integrability of an equation in n variables at the price of extra conditions (differential equations) in less than n variables. The algorithm EXACT is to be modified as follows line 7: take only functions u^{α} depending on all independent variables line 19: ∫ Δ dx now includes as well the integration of terms involving only functions of fewer variables than occurred in the original Δ. If this ”generalized integration” of all terms involving any unknown function of x as explained below is not possible then Δ cannot be integrated, i.e. nil is to be returned.

The method The remaining expression Δ to be integrated in line 19 is regarded as a sum of products (if necessary as a single term of one or more factors). Each term T_{i} is written as a product T_{i} = V_{i}W_{i} of a product V_{i} of x–independent factors and a product W_{i} of x– dependent factors. The integral is ∫ T_{i} dx = ∫ V_{i}W_{i} dx = V_{i}c_{r} with the extra condition D_{x}c_{r} = W_{i} and a new function c_{r} that depends on exactly all independent variables occurring in W_{i}. The generalized integration is not successful if any W_{i} does depend on all variables and does contain an unknown function of x. Example: If the original Δ has x, y as independent variables and involves a term sin(xy) g with g = g(x) depending only on x then this term cannot be factorized into an x–independent factor V and an x–dependent factor W which does not depend on all variables. The generalized integration is therefore not successful. If one would introduce a new function c_{r} with D_{x}c_{r} = sin(xy) g then c_{r} would have to depend on all variables x and y which would make this integration useless as no information is gained.

Improvements In order to reduce the number of new functions and extra conditions that are added and that have to be solved afterwards, the following refinements have been implemented in the program Crack.

1. Instead of determining V_{i} to be the product of x–independent factors, it now can contain powers of x, i.e. it can be a polynomial of x.

2. Terms T_{i}, T_{j}, … with equal W –factor are added to T := (V_{i} + V_{j} + …)W.

3. For the integration of T only one new function c which depends on all variables occurring in W has to be introduced in the following way. We assume we want to integrate (v_{0} +v_{1}x+…+v_{n}x^{n})W. Only for simplicity in notation in the following equation we use ∂_{x}^{k}c = c^{(k)}:

∫ W dx = c^{(n−1)},
∫ xW dx = −c^{(n−2)} + xc^{(n−1)},
∫ x^{2}W dx = c^{(n−3)} − 2xc^{(n−2)} + x^{2}c^{(n−1)},
…
∫ x^{k}W dx = (−1)^{k}k!c^{(n−k−1)} + … − kx^{k−1}c^{(n−2)} + x^{k}c^{(n−1)}
= ∑_{m=0}^{k}(−1)^{m} \frac{k!}{(k − m)!}x^{k−m}c^{(n−m−1)}
…
∫ x^{n}W dx = ∑_{m=0}^{n}(−1)^{m} \frac{n!}{(n − m)!}x^{n−m}c^{(n−m−1)}

5 Implementations

The first implementation known to the author was written in the computer algebra system Formac in 1987 [11]. A later version was implemented by A. Brand and the author in the computer algebra system Reduce. This version is currently (1999) used in the package Crack for the solution of over-determined PDE systems. In autumn 1998 a Maple version was implemented by M. White in collaboration with the author. Khai Vu independently translated the integration module of Crack into Maple (see [1]) but without the part of section 4.

Acknowledgement

The author thanks Andreas Brand for years of collaboration on algorithmic and implementation issues regarding the package Crack. Malcolm MacCallum is thanked for comments on the manuscript.

References

[1] J. Carminati and Vu K. Symbolic computation and differential equations: Lie symmetries. to appear in J. of Symb. Comp.

[2] A. Head. Lie and BigLie, programs to find lie symmetries of differential equations.

    http://www.cmst.csiro.au/LIE/LIE.htm.

[3] W. Hereman. Symbolic software for lie symmetry analysis. In N.H. Ibragimov, editor, CRC Handbook of Lie Group Analysis of Differential Equations, Vol. 3: New trends in Theoretical Developments and Computational Methods, chapter 13, pages 367–413. CRC Press, Boca Raton, Florida, 1996.

[4] E. Hubert. Essential components of algebraic differential equations. To appear in the Journal of Symbolic Computation, 2000.

[5] E.L. Mansfield. The differential algebra package diffgrob2. Mapletech, 3:33–37, 1996.

[6] G.J. Reid, A.D. Wittkopf, and A. Boulton. Reduction of systems of nonlinear partial differential equations to simplified involutive forms. European Journal of Applied Mathematics, 7:604–635, 1996.

[7] G.J. Reid, A.D. Wittkopf, and P. Lin. Differential-elimination completion algorithms for differential algebraic equations and partial differential algebraic equations. to appear in Studies in Applied Mathematics, 1995.

[8] J. Sherring. Symmetry and computer algebra techniques for differential equations. PhD thesis, La Trobe University, 1993.

[9] J. Sherring, G. Prince, and M. Jerie. Symmetry determination and linear differential equation package.

http://www.latrobe.edu.au/www/mathstats/Maths/Dimsym.

[10] T. Wolf. Combining Gröbner basis computations and integrations. in preparation.

[11] T. Wolf. A package for the analytic investigation and exact solution of differential equations. In Proceedings of EUROCAL 87, Lect. Notes of Comp. Sci., volume 378, pages 479–91, 1989.

[12] T. Wolf. A comparison of four approaches to the calculation of conservation laws. accepted for publication in the Europ. J. of Appl. Math., 1998.

[13] T. Wolf, A. Brand, and M. Mohammadzadeh. Computer algebra algorithms and routines for the computation of conservation laws and fixing of gauge in differential expressions. J. Symb. Comp., 27:221–238, 1999.