ReLie: a Reduce program for Lie group analysis of differential equations

Lie symmetry analysis provides a general theoretical framework for investigating ordinary and partial differential equations. The theory is completely algorithmic even if it usually involves lengthy computations. For this reason, many computer algebra packages have been developed along the years to automate the computation. In this paper, we describe the program ReLie, written in the Computer Algebra System Reduce, which since 2008 is an open source program (http://www.reduce-algebra.com) and is available for all platforms. \relie is able to perform almost automatically the needed computations for Lie symmetry analysis of differential equations. Its source code is freely available at the url http://mat521.unime.it/oliveri. The use of the program is illustrated by means of some simple examples; nevertheless, it is to be underlined that it provides effective also for more complex computations where one has to deal with very large expressions.


Introduction
A general unified approach to deal with ordinary as well as partial differential equations is provided by the study of their continuous symmetries, i.e., transformations of their solution manifold into itself. Symmetry analysis of differential equations originated in the nineteenth century with Sophus Lie [66,65] who recognized that many special integration theories of differential equations are a consequence of the invariance under one-parameter continuous groups of transformations. Such transformations establish a diffeomorphism on the space of independent and dependent variables, mapping solutions of the equations to other solutions. Any transformation of the independent and dependent variables in turn induces a transformation of the derivatives. As Lie himself showed, the problem of finding the Lie group of point transformations leaving a differential equation invariant consists in solving a linear system of determining equations for the components of its infinitesimal generators. Lie's theory (many textbooks and monographies are available, for instance [92,51,87,12,102,55,56,57,88,3,16,9,72,10,13]) allows for the development of systematic procedures leading to the integration by quadrature (or at least to lowering the order) of ordinary differential equations [6], to the determination of invariant solutions of initial and boundary value problems [87,96,55,56,81,82,83,78,84], to the derivation of conserved quantities [74,60,8], to the algorithmic construction of invertible point transformations mapping the differential equations into equivalent forms easier to handle [61,30,31,32,27,79,80,39,40], or to design efficient numerical schemes [89,94,4], . . .
Lie's classical theory is a source for various generalizations. Among these generalizations there are the nonclassical method first proposed by Bluman and Cole [11], and now part of the more general method of differential constraints [35,37], the potential symmetries [7], the nonlocal symmetries [45,63], the generalized symmetries [87], which in turn generalize contact symmetries introduced by Lie himself, the equivalence transformations [92,67,53,54,103,104,93], to quote a few. A recent introduction is that of approximate symmetries [2,52,36,29,38] for differential equations containing small terms.
The application of Lie's theory to differential equations is completely algorithmic; nevertheless, it usually involves a lot of cumbersome and tedious calculations. Today, many powerful Computer Algebra Systems (CAS) (either commercial or open source) are available, and the needed algebraic manipulations can be rapidly done almost automatically. In fact, many specific packages for performing symmetry analysis of differential equations are currently available in the literature [98,99,47,49,100,57,50,3,17,105,15,23,101,24,95,59].
In this paper, we describe the package ReLie, written in the computer algebra system Reduce [48]. The routines contained in ReLie allow the user to easily compute (exact and approximate) Lie point symmetries and conditional symmetries, as well as contact transformations, variational sym-metries, and equivalence transformations for classes of differential equations containing arbitrary elements.
Reduce is a general purpose computer algebra system whose development started by Anthony Hearn. It is written in a Lisp dialect (Portable Standard Lisp); since December 2008 it is an open source program (available for all platforms at the url http://www.reduce-algebra.com).
In the following, we assume that every potential user of ReLie is familiar with Lie symmetry analysis of differential equations. Nevertheless, in order to fix the notation and keep the paper self-contained, a brief sketch of the underlying theory is given. The use of ReLie requires the user to have a minimal experience in working with the algebraic mode of Reduce. Therefore, in this paper only the necessary details of the use of ReLie will be presented. The source code of ReLie, as well as the user's manual, can be found at the url http://mat521.unime.it/oliveri.
The origin of this package dates back to 1994 when the author developed some routines useful to manage the lengthy expressions needed to determine the Lie point symmetries of differential equations; this set of routines constantly grew through the years providing new capabilities, and now constitutes a package able to perform almost automatically much of the work. Moreover, the program can be used also in interactive mode mimicking the steps one has to apply with pencil and paper but with the benefits of using a computer algebra system: this can represent a useful support in a higher course on symmetry analysis of differential equations or in all those situations (e.g., when one looks for conditional symmetries) where the determining equations can not be automatically solved and some special assumptions are needed.
The plan of the paper is the following. In Section 2, we review the basic elements of Lie group theory of differential equations with reference to point symmetries (subsection 2.1), conditional symmetries (subsection 2.2), contact symmetries (subsection 2.3), variational symmetries (subsection 2.4), approximate point, conditional, contact and variational symmetries (subsection 2.5), and equivalence transformations (subsection 2.6). Then, in Section 3, we consider some examples in order to illustrate the use of Re-Lie. Finally, Section 4 contains a description of the global variables and the main routines contained in the package.
2 Basic elements of the theory 2

.1 Lie point symmetries
Within the framework of Lie group analysis, given a system of differential equations, say ∆ x, u, u (r) = 0, where x ∈ X ⊆ R n is the set of the independent variables, u ∈ U ⊆ R m the set of the dependent variables, and u (r) the set of all (partial) derivatives of the u's with respect to the x's up to the order r ≥ 1, one is interested to find the admitted group of Lie point symmetries. A Lie point symmetry is characterized by its infinitesimal generator [92,87,12] where ξ i and η α are the infinitesimals of the group. In dealing with differential equations, we need to extend (or prolong) the action of the Lie group to the r-th order jet space, whose coordinates are the independent and dependent variables as well as the derivatives of the latter with respect to the former up to the order r. Transformations for derivatives are obtained by requiring that the contact conditions are preserved by the transformation (roughly speaking, the transformed derivatives have to be the derivatives of the transformed dependent variables with respect to the transformed independent variables). First order prolongation of the infinitesimal generator reads and, in general, higher order extended infinitesimal generator has the form where The point transformations leaving system (1) invariant are found by means of the straightforward Lie's algorithm, which requires that the r-th order prolongation Ξ (r) of the vector field (2) acting on (1) is zero along the solutions of (1), i.e., Notice that the invariance condition (4) involves polynomials in the variables u (r) if all these derivatives occur polynomially in the differential equations: in fact, the prolonged infinitesimals are always polynomials in the derivatives u (r) . Since we evaluate the invariance condition on the differential equations, the derivatives which appear in the invariance condition are independent. Therefore, it vanishes if and only if all the coefficients of such polynomials (involving the infinitesimals, their derivatives, the independent and the dependent variables) are zero, whereupon we obtain an overdetermined set of linear and homogeneous partial differential equations (determining equations) for the infinitesimals ξ and η, whose integration provides the generators of Lie point symmetries admitted by (1). The infinitesimal generators of a Lie group, whose components are solutions of a linear homogeneous system of partial differential equations, span a vector space which can be finite or infinite-dimensional; by introducing an operation of commutation between two infinitesimal generators, this vector space gains the structure of a Lie algebra [28,34].

Q-conditional symmetries
It is well known that some differential equations of interest in the applications possess very few Lie symmetries. In 1969, Bluman and Cole [11] introduced a generalization of classical Lie symmetries, and applied their method (called nonclassical) to the linear heat equation. The basic idea was that of imposing the invariance to a system made by the differential equation at hand, the invariant surface condition and the differential consequences of the latter [90,91,25,64,76,1,75,26,97]. This method leads to nonlinear determining equations whose general integration is usually difficult. Nonclassical symmetries are now part of conditional symmetries, i.e., symmetries of differential equations where some additional differential constraints are imposed to restrict the set of solutions; some recent applications of conditional symmetries to reaction-diffusion equations can be found in [18,22,19,20,21].
Given a differential equation and considering the invariant surface conditions Q = 0, where the Q-conditional symmetries are found by requiring where M is the manifold of the jet space defined by the system of equations where 1 ≤ |k| = k 1 + . . . + k n ≤ r − 1.
Of course, a (classical) Lie symmetry is a (trivial) Q-conditional symmetry. However, differently from Lie symmetries, all possible Q-conditional symmetries of a differential equation form a set which is not a Lie algebra in the general case. Moreover, if the infinitesimals of a Q-conditional symmetry are multiplied by an arbitrary smooth function we have still a Q-conditional symmetry. This means that we can look for Q-conditional symmetries in n different situations: When m > 1 (more than one dependent variable), one can look for conditional symmetries of q-th type (1 ≤ q ≤ m) [19,20,21], by requiring the condition Ξ (r) (∆) where M q is the manifold of the jet space defined by the system of equations It is obvious that we may have in principle m q different sets of Q-conditional symmetries of q-th type, whereas for q = 0 we simply have classical Lie symmetries. Moreover, when q < m it may result easier to find solutions to the nonlinear determining equations. Q-conditional symmetries allow for symmetry reductions of differential equations, and may provide explicit solutions not obtainable with classical symmetries.

Contact transformations
Besides point transformations whose infinitesimals depend at most on independent and dependent variables, contact transformations [33,87,62] play an important role in many applications. Nevertheless, true contact transformations exist for differential equations involving only one dependent variable (for more than one dependent variable, they are prolongations of point transformations), and their infinitesimal generator has the form The operator (11) characterizes a group of contact transformations if and only if there exists a function Ω(x, u, u (1) ) (characteristic function) such that Prolongations of contact transformations for higher order derivatives are computed with the usual formulas.
Proper contact transformations (i.e., transformations that are not prolongations of point transformations) are the ones determined by a characteristic function Ω which is not linear in the first order derivatives.

Variational symmetries
In dealing with differential equations, conservation laws have a deep relevance, since they express the conservation of physical quantities such as mass, momentum, angular momentum, energy, electrical charge. They are also important due to their use in investigating integrability, existence, uniqueness and stability of solutions, or in implementing efficient numerical methods of integration [46,58].
In 1918, Emmy Noether [74] presented her celebrated procedure (Noether's theorem) to find local conservation laws for differential equations arising from a variational principle. Noether proved that a point symmetry of the action functional (action integral) provides a local conservation law (i.e., a divergence expression) through an explicit formula that involves the infinitesimals of the point symmetry and the Lagrangian itself.
Consider a functional given by an integral over a domain Ω ⊂ R n , say where L(x, u, u (r) ) is a Lagrangian function. By imposing the functional J to be stationary (for suitable variations of u and u (r) vanishing on the boundary of the domain Ω), we derive the Euler-Lagrange equations where (15) is the Euler operator with respect to u α ; here and in the following, if needed, we use the Einstein convention of sums over repeated indices.
Noether's theorem [74] states that if we know a Lie group of point transformations with generator Ξ leaving the action integral invariant, and this holds true if the condition for some functions φ i (x, u, u (r−1) ) (i = 1, . . . , n) is satisfied, then the conservation law where holds for any solution u(x) of Euler-Lagrange equations. For n = 1, the conservation law provides a first integral of ordinary differential equations. Boyer [14] extended Noether's theorem in order to construct conservation laws arising from invariance under generalized symmetries [87], i.e., symmetries with infinitesimals depending on higher order derivatives (see [10]).

Approximate symmetries
Let ∆ x, u, u (r) ; ε = 0 (18) be a differential equation of order r involving a small parameter ε. It is well known that any small perturbation of an equation usually destroys some important symmetries, and this restricts the applicability of Lie group methods to differential equations arising in concrete applications. On the other hand, differential equations containing small terms are commonly and successfully investigated by means of perturbative techniques [73]. Therefore, it is desirable to combine Lie group methods with perturbation analysis, i.e., to establish an approximate symmetry theory. In the literature, there are two widely used approaches to approximate symmetries: one has been proposed in 1988 by Baikov, Gazizov and Ibragimov [2] (see also [52]) and another one in 1989 by Fushchich and Shtelen [36]. Both approaches have pros and cons [29]. To overcome these problems, recently, a consistent approach has been proposed [29], that is consistent with perturbation theory, and allows to extend all the relevant features of Lie group analysis to an approximate context. In what follows this latter approach is described. If one looks for classical Lie point symmetries, in general it is not guaranteed that the infinitesimal generators depend on the parameter ε. Nevertheless, the occurrence of terms involving ε has dramatic effects since one loses some symmetries admitted by the unperturbed equation ∆ x, u, u (r) ; 0 = 0.
Looking for solutions in the form the differential equation writes as Now, let us consider a Lie generator and assume that the infinitesimals depend on the small parameter ε. By using the expansion (20) of the dependent variables only, we have for the infinitesimals where R being a linear recursion operator satisfying product rule of derivatives and such that where k ≥ 0, j = 1, . . . , m, |τ | = τ 1 + · · · + τ m . Thence, we have an approximate Lie generator where Since we have to deal with differential equations, we need to prolong the Lie generator to account for the transformation of derivatives. This is done as in classical Lie group analysis of differential equations, i.e., the derivatives are transformed in such a way the contact conditions are preserved. Of course, in the expression of prolongations, we need to take into account the expansions of ξ i , η α and u α , and drop the O(ε p+1 ) terms.
The approximate (at the order p) invariance condition of a differential equation reads In the resulting condition we have to insert the expansion of u in order to obtain the determining equations at the various orders in ε. The Lie generator Ξ (0) is always a symmetry of the unperturbed equa- give the deformation of the symmetry due to the terms involving ε. Not all symmetries of the unperturbed equations are admitted as the zeroth terms of approximate symmetries; the symmetries of the unperturbed equations that are the zeroth terms of approximate symmetries are called stable symmetries [2].
Remark 1 If Ξ is the generator of an approximate Lie point symmetry of a differential equation, εΞ is a generator of an approximate Lie point symmetry too.
By the same arguments as in classical Lie theory of differential equations it can be proved that the approximate Lie point symmetries of a differential equation are the elements of an approximate Lie algebra [2,52,29].
Of course, if one considers differential equations containing small terms, it is possible to construct approximate conditional symmetries [41,43], the only difference with respect to exact symmetries being the structure of Lie generator which follows the approach described above. In analogy with the exact case, for p-th order approximate conditional symmetries we have n different cases: We may also consider approximate contact transformations for scalar differential equations containing small terms by assuming the characteristic function Ω to depend on the small parameter ε. Limiting to first order approximate contact symmetries, the expansion for Ω reads Finally, also approximate variational symmetries for Lagrangians containing small terms can be defined.

Equivalence transformations
There are situations where the differential equations at hand contain unspecified functions (arbitrary elements); therefore, it is convenient to consider a class E(p) of differential equations involving some arbitrary continuously differentiable functions p k (x, u) (k = 1, . . . , ℓ).
Let us consider a class of differential equations whose elements are given once we fix the functions p k , where p (s) , denotes the set of all derivatives up to the order s of the p's with respect to their arguments; these arguments may be the independent variables, the dependent ones and the derivatives of the latter with respect to the former up to a fixed order q: let us denote the arguments of the arbitrary elements p as z ≡ (z 1 , . . . , z N ).
To face this problem, it is convenient to consider equivalence transformations, i.e., transformations that preserve the differential structure of the equations in the system but may change the form of the constitutive functions and/or parameters [92,67,53,54,103,104,93].
A one-parameter Lie group of equivalence transformations [92] of a family E(p) of PDEs is a one-parameter Lie group of transformations given by a being the parameter, which is locally a C ∞ diffeomorphism and maps a class E(p) of differential equations into itself; thus, it may change the differential equations (the form of the arbitrary elements therein involved) but preserves the differential structure. A common assumption consists in considering transformations of the independent and dependent variables independent of the arbitrary elements p.
In an augmented space A ≡ X ×U ×P ⊆ R n ×R m ×R ℓ [92,67], where the independent variables, the dependent variables and the arbitrary functions live, respectively, the generator of the equivalence transformation, involves also the infinitesimals µ j (x, u, p) accounting for the arbitrary functions p j . The search for continuous equivalence transformations can be exploited by using the Lie's infinitesimal criterion [92]. The first prolongation of Ξ writes as with where p j,γ = ∂p j ∂z γ and ζ γ is the infinitesimal generator of z γ , along with the Lie derivatives D Dx i already defined in (3) and Higher order prolongations are immediately obtained by recursion, In the augmented space A, the arbitrary functions determining the class of differential equations are assumed as dependent variables, and we require the invariance of the class in this augmented space via Lie's infinitesimal criterion [92]: If we project the symmetries on the space Z ≡ X × U of the independent and dependent variables (this is always possible if the infinitesimals of independent and dependent variables are assumed to be independent of p), we obtain a transformation changing an element of the class of differential equations to another element in the same class (same differential structure but in general different arbitrary elements) [85,86,44]. Such projected transformations map solutions of a system in the class to solutions of a transformed system in the same class. If some arbitrary elements do not depend upon some variables, the differential equations at hand have to be supplemented with auxiliary conditions. For instance, the auxiliary conditions for the equivalence transformations of the class of differential equations [54] where f and g are unspecified functions of their arguments, are: More in general, taking the transformations of the independent and dependent variables as functions of the arbitrary elements p too [71], in the expression of prolongation we have to replace the Lie derivative

Lie remarkable equations
In [70], within the framework of inverse Lie problem, strongly and weakly Lie remarkable differential equations have been defined; relevant examples of such equations have been studied in [69,68,42]. Their analysis involves the study of the rank of the distribution of prolongations of a Lie algebra of generators.
In this Section, we illustrate the use of Reduce [48] package ReLie by considering some simple examples; more complex applications can be considered as well. The package ReLie is a collection of (algebraic) routines that allow the user to investigate ordinary and partial differential equations within the framework of Lie symmetry analysis. By using the program, it is possible to compute almost automatically: Lie point symmetries, conditional symmetries, contact symmetries, variational symmetries (all these symmetries may be either exact or approximate) of differential equations, and equivalence transformations of classes of differential equations containing arbitrary elements. Moreover, the program implements functions for computing Lie brackets, the commutator table of a list of Lie generators, and the distribution of an algebra of Lie symmetries (useful in the context of inverse Lie problem [70,69,68,42]). Remarkably, the program can be used interactively in all the cases where the determining equations are not automatically solved (for instance, when one looks for conditional symmetries or in some group classification problems).
After entering Reduce, assuming that we are in the directory containing the source file relie.red, the set of routines implemented in the package ReLie becomes available after the statement in "relie.red" $ The loading of ReLie is successful if the following output is displayed ReLie, version 3.0 A REduce program for Lie group analysis of differential equations (c) Francesco Oliveri (foliveri@unime.it) -2020 Last update August 9, 2020 http://mat521.unime.it/oliveri/ Alternatively, one may include the package in the Reduce image by issuing in a Reduce session the statements faslout 'relie $ in "relie.red" $ faslend $ If all works, ReLie is loaded through the command load_package relie $ According to the task we are interested to, some input data have to be provided to the program. The minimal required set of data is given by a positive integer value for the global variable jetorder (the maximum order of prolongation of the Lie generator), the list xvar of the independent variables, and the list uvar of the dependent variables.
The number of independent as well as dependent variables is arbitrary; moreover, the identifiers of independent and dependent variables are arbitrary provided that no conflict arises with reserved words of Reduce or identifiers used by ReLie (see Section 4 for a list of internal global variables used by ReLie). The value of jetorder is not bounded by ReLie.
The first function one has to call is relieinit(), that defines and initializes all the needed objects required to perform Lie group analysis. It checks the input and possibly displays some warnings. If only jetorder, xvar and uvar have been assigned, then calling relieinit() produces the output: The list 'diffeqs' of differential equation(s) is missing! The list 'leadders' of leading derivative(s) is missing! Check 'diffeqs' and/or 'leadders'! You may only call relieprol() or generatealgebra(k) (k=1,2,3).
In fact, we did not define the list diffeqs of differential equations and the list leadders of leading derivatives. Nevertheless, we can compute the jetorder-th prolongation of the vector field generating a Lie group of point transformation. The infinitesimal generators of the independent and dependent variables are automatically defined by the program (therefore, the user is not requested to set them), and stored in the list allinfinitesimals. The infinitesimals of the independent variables are denoted by xi followed by the identifiers in xvar; analogously, the infinitesimals of the dependent variables are denoted by eta followed by the identifiers in uvar. According to the Example 2, the list allinfinitesimals is a list of two elements, each element being the list {xi t, xi x, eta u} (this redundancy is used in order to have a unified coding for dealing also with other kinds of symmetries, see below). The dependencies of the elements in allinfinitesimals upon the independent and dependent variables (the elements in xvar and uvar, respectively) are automatically set by ReLie.
By calling the ReLie function relieprol(), as a result we get the list prolongation. This list has two elements: the first one is a list containing the coordinates of the jet space, the second one the list of the corresponding infinitesimal generators. Thus, the second element of prolongation contains the components of the prolonged vector field up to jetorder-th order.
ReLie internally represents the derivatives of the dependent variables upon the independent variables by appending to the identifiers of the dependent variable the underscore followed by the identifiers of the involved independent variables. For mixed derivatives the order of the independent variables in the internal representation of derivatives reflects the order in the list xvar. For instance, if xvar is {t,x} and uvar is {u}, the left-hand side of the equation (Benjamin-Bono-Mahoney equation) has to be represented as u_t + u_x + u*u_x -u_txx .

Computing Lie point symmetries of differential equations
To compute the Lie point symmetries admitted by differential equations, in addition to jetorder, xvar and uvar, we have to provide: 1. the list diffeqs of the left-hand sides of the differential equations to be studied which are assumed with zero right-hand sides; 2. the list leadders of some derivatives appearing in the differential equations: when computing the invariance conditions, the elements in the list leadders are removed by solving the differential equations with respect to them.
Both lists, diffeqs and leadders, need to have the same length, otherwise a warning is displayed; moreover, the differential equations must be solvable with respect to the leading derivatives. The user should be careful in the choice of leadders, especially in the case of systems of differential equations, in order to guarantee that the differential equations stored in the list diffeqs can be solved with respect to the elements in leadders.
The list symmetries contains an element (when more than one solution set is possible for particular choices of the involved parameters this list has one element for each solution set, see later) which in turn is a list of 4 elements: 1. the first one is a list of conditions that are still unsolved (in this case the list is empty since no condition remains unsolved); 2. the second one is a list giving the solution to the determining equations, i.e., the expressions of the infinitesimals; 3. the third one is a list containing the parameters involved in the solution (in this case k 1 and k 2); 4. the fourth one is a list of expressions which can not vanish (in this case the list is empty).
The parameters involved in the symmetries of Blasius equation are constant; symmetries of partial differential equations may involve arbitrary functions that ReLie denotes as f 1, f 2, . . . ; we can check the dependencies of such parameters by passing the name of the function to the procedure fargs() (described in Crack [105]).
Since the infinitesimals are solutions of linear homogeneous differential equations, they are linear combinations of fundamental solutions. We can compute these fundamental solutions through the use of the function reliegen(). This function requires two arguments: an integer and a list. If reliesolve() returns only one solution this integer is 1; when more than one solution is available (this may occur in group classification problems or in the search of conditional symmetries), the integer selects the solution. We have a list where each element is the list of the components of the vector field generating a symmetry (the components of the infinitesimals of the independent and dependent variables); the vector fields here characterize a scaling and a translation of the independent variable, respectively: i.e., spanning a 6-dimensional Lie algebra. The symmetry corresponding to the function f 1, renders the algebra of symmetries of linear heat equation infinite-dimensional.
The function reliedet() assumes that the invariance conditions are polynomials in the derivatives: this is true if all derivatives appear polynomially in the differential equations at hand. In the case in which the differential equations we are investigating contain some derivatives not in polynomial form, ReLie constructs correctly the determining equations if we inform the program about the derivatives not occurring in polynomial form. These derivatives have to be the elements of the list nonpolyders.

Remark 3
The use of ReLie is not limited to simple cases. For instance, consider the Navier-Stokes-Fourier equations of a gas in a rotating (with constant angular velocity ω along the vertical axis) frame reference and subject to gravity, where ρ is the mass density, v ≡ (v 1 , v 2 , v 3 ) the velocity, T the temperature, F i the i-th component of the sum of external and inertial force, g the gravitational acceleration, K the Boltzmann constant, m the mass of a single particle, and the coefficient α is an appropriate constant which describes the interaction between the Maxwellian molecules. The computation of the Lie symmetries (spanning a 12th dimensional Lie algebra) is done automatically in few second by means of ReLie.

An example of group classification
ReLie is able to face the problem of group classification of differential equations. In fact, in certain cases it may be interesting to check if special instances of arbitrary constants and/or functions involved in the differential equations lead to different sets of Lie symmetries. This problem is faced by including the parameters we want to check in the list freepars, as illustrated in the following example. Moreover, we may define the list nonzeropars containing parameters (constants or functions involved in the differential equations) or expressions that can not be vanishing.

Commutator table
ReLie is able to compute the commutator table of a set of infinitesimal generators spanning a Lie algebra. Let us illustrate how to do with an example.

Example 7 Consider the Burgers' equation
To compute its Lie symmetries we provide to ReLie the following input: The commutator  The list qcond tells to the program which invariant surface conditions have to be used to restrict the solution manifold:

Computation of conditional symmetries
• qcond:={1,2} if the invariant surface conditions of both dependent variables have to be used; • qcond:={1} if the invariant surface condition of the first dependent variable has to be used; • qcond:={2} if the invariant surface condition of the second dependent variable has to be used.
Looking for conditional symmetries leads to nonlinear determining equations, whereupon it is not unusual that reliesolve() may fail to recover the solution. In such cases, it is necessary to make some special assumptions on the infinitesimals and/or try to solve interactively the determining equations.
of the linear heat equation As a result, we get Therefore, we have: where f 1 (t, x) and f 2 (t, x) satisfy the differential equations Remark 5 When computing conditional symmetries we have to be careful in the choice of leadders, since also the invariant surface conditions and their differential consequences are solved with respect to some derivatives.

Computation of contact symmetries
For computing contact symmetries of a scalar differential equation we need to assign the value 1 to the variable contact. The functions that we have to call do not change, as the following example shows.
Example 9 Let us compute the contact symmetries of the ordinary differential equation d 3 u dx 3 = 0. The code is as follows.

Computation of variational symmetries and associated conservation laws
For computing variational symmetries, besides assigning jetorder, xvar and uvar we need to assign the list lagrangian with only one element, and set the variable variational to 1. The functions that we have to call do not change, as the following example shows.

Computation of approximate symmetries
Approximate Lie point symmetries of differential equations containing small terms are computed according to the approach described in [29]. The small parameter must be denoted by epsilon. In addition to the input data described before, the user has to define the integer approxorder for the order of approximation, and a rule, say let epsilon**(approxorder+1) = 0. The same functions providing the necessary computations for exact Lie point symmetries work in this case too.
To compute approximate conditional, contact or variational symmetries the only difference with respect to the corresponding exact case consists in setting a positive value to approxorder and defining the rule let epsilon^(approxorder+1) = 0 $

Computation of equivalence transformations
If we consider a class of differential equations, and want to determine equivalence transformations some more input data are necessary. In particular, in addition to the data already discussed for the computation of Lie point symmetries, we need to fix: 1. the list arbelem of the arbitrary elements involved in the differential equations; 2. the integer arborder denoting the highest order of the derivatives of the arbitrary elements with respect to their arguments; 3. the integer zorder characterizing the variables the arbitrary elements depend on; for instance, if zorder is 0, then the arbitrary elements depend at most on the independent and dependent variables; if zorder is 1, then the arbitrary elements depend at most on the independent, dependent variables and first order derivatives, . . . .
If we want to remove the dependence of some arbitrary elements on some variables, we need to add auxiliary conditions to the differential equations at hand. The infinitesimal generators of the arbitrary elements are automatically computed by the program (so the user is not requested to set them), and stored, together with the infinitesimals of independent and dependent variables in the list allinfinitesimals. The infinitesimal generator of an arbitrary element is denoted by mu followed by the name of the arbitrary element. Let us illustrate with two different examples how ReLie works in such a situation.
Remark 7 By default, the infinitesimals for the independent and dependent variables do not depend on arbitrary elements; if we are interested to general equivalence transformations where also the infinitesimal generators of the independent and dependent variables depend on the arbitrary elements [71], then we have to add the statement generalequiv:=1 $

Inverse Lie problem
Here we show with a simple example how ReLie can be used for investigating the Lie remarkability of differential equations [70,69,68,42]. The following code shows how one can easily verify that Monge-Ampère equation describing a surface in R 3 with zero Gaussian curvature, is strongly Lie remarkable [70]: jetorder:=2 $ xvar:={x,y} $ uvar:={u} $ relieinit() $ generatealgebra(2) $ After setting the necessary objects (jetorder, xvar and uvar), and calling relieinit(), we generate the vector fields spanning the affine Lie algebra in R 3 . Then, calling testrank(generators), we ascertain that the rank of the second order distribution is maximal (equal to the dimension of second order jet space); the call inverselie(8) computes the list allminors of all minors of order 8 of such a distribution; all minors vanish by evaluating them on the differential equation, that can be verified by calling allminors:=sub(u_yy=u_xy^2/u_xx,allminors) $ As a last check, rank(sub(u_yy=u_xy^2/u_xx,distribution); returns 7, which is the dimension of the manifold in the second order jet space characterized by Monge-Ampère equation.

Inside ReLie: global variables and routines
ReLie uses some global variables to perform the various tasks: they can be distinguished among input variables (that the user needs to set before starting computation), output variables (computed by ReLie and of interest to the user), and intermediate variables (computed by ReLie and in general not of interest to the common user).

Input variables
• approxorder: maximum order of approximate symmetries of equations containing a small parameter; by default it is 0, i.e., exact symmetries; the small parameter involved in the approximate symmetries must be denoted by epsilon; when looking for approximate symmetries the user has to define the rule let epsilon**(approxorder+1)=0; • arbelem: list of the arbitrary elements (only for equivalence transformations); by default it is an empty list; • arborder: maximum order of derivatives of arbitrary elements (only for equivalence transformations); by default it is −1, i.e., point symmetries; • contact: set to 1 for contact symmetries; by default it is 0; • diffeqs: list of the left-hand sides of differential equations (with vanishing right-hand sides); • freepars: list of arbitrary constants or functions involved in the differential equations (for group classification problems); by default it is an empty list; • generalequiv: set to 1 for general equivalence transformations where all infinitesimals depend on independent and dependent variables and arbitrary elements; the default value is 0, meaning that the infinitesimals of independent and dependent variables do not depend upon the arbitrary elements; • jetorder: maximum order of derivatives in differential equations; • lagrangian: a list with only one element corresponding to the Lagrangian (it is necessary to set variational to 1); • leadders: list of the leading derivatives; diffeqs are solved with respect to them; • nonclassical: set to a value between 1 and the number of independent variables (only for conditional symmetries); by default it is 0, i.e., classical symmetries; • nonpolyders: list of derivatives not occurring in polynomial form in the differential equations; by default it is an empty list; • nonzeropars: list of arbitrary constants or functions involved in the differential equations that can not vanish; by default it is an empty list; • qcond: list of indexes of dependent variables whose invariant surface conditions have to be used for computing conditional symmetries; • uvar: list of the dependent variables; • variational: set to 1 if variational symmetries of a Lagrangian are needed; by default it is 0; • xvar: list of the independent variables; • zorder: maximum order of derivatives of uvar with respect to xvar the elements in arbelem depend on (only for equivalence transformations); if zorder is set to 0, the arbitrary elements depend on xvar and uvar; zorder cannot exceed jetorder.

Output variables
• allinfinitesimals: list of two lists; the first sublist is the list of the infinitesimals, in order, of the independent variables, dependent variables and arbitrary elements (the latter in the case of equivalence transformations); the second sublist is the list of various terms (constants and functions) involved in the expression of infinitesimals; • allminors: list of minors of a given order extracted from a matrix; returned by the function minors(m,k), where m is a matrix and k a positive integer, or by the function inverselie(k) that takes the jetorder-th distribution of generators as the matrix from which the minors of order k are extracted; • arbconst: list of arbitrary constants involved in the expression of the symmetries; • arbfun: list of arbitrary functions involved in the symmetries; • cogenerators: list of the functions entering the definition of variational symmetries and corresponding to the infinitesimal generators (produced by reliegen()); • commtable: table of commutators of a list of vector fields; • deteqs: list of the determining equations (produced by reliedet()); • distribution: matrix of the jetorder-th distribution of a list of generators (produced by reliedistrib()), i.e., a matrix where each row is the prolonged vector field evaluated on one of the provided infinitesimal generators; • fluxes: list of the fluxes of the conservation law corresponding to a Lie generator (computed by relieclaw()); • generators: list of the infinitesimal generators of the finite Lie algebra admitted by the differential equations at hand (produced by reliegen()); the list generators may also been obtained by calling generatealgebra(k), where k can be 1 (algebra of isometries), 2 (algebra of affine transformations) or 3 (algebra of projective transformations); of course, it is necessary to set jetorder, xvar and uvar before calling generatealgebra(k); • invcond: list of the invariance conditions of the differential equations at hand (produced by relieinv()); • nzcomm: list of non-zero commutators of a list of vector fields; • prolongation: list of two lists: the first one is the list of the coordinates of the jet space, the second one the list of the corresponding infinitesimals (produced by relieprol()); • splitsymmetries: list of lists: the k-th element is a list containing the infinitesimals corresponding to the k-th element of generators (produced by reliegen()); the list splitsymmetries may also been obtained by calling generatealgebra(k), where k can be 1 (algebra of isometries), 2 (algebra of affine transformations) or 3 (algebra of projective transformations); of course, it is necessary to set jetorder, xvar and uvar before calling generatealgebra(k); the list splitsymmetries is used internally by the functions reliedistrib(), inverselie() and testrank(); • symmetries: list of the infinitesimals admitted by the differential equations at hand (produced by reliesolve()).

Intermediate variables
• jet: list of three lists: indices for computing the infinitesimals and their prolongations, coordinates of jet space and their internal representation; • jetapprox: list of two lists: list of independent variables and expansions of dependent variables and their derivatives, and list of their internal representation (only for approximate symmetries); • jetequiv: list of three lists: indices for computing the infinitesimals and their prolongations for arbitrary elements, arbitrary elements, and their internal representation (only for equivalence transformations); • jetsplit: list of two lists: indices for computing the infinitesimals and their prolongations, list of independent variables, zeroth order dependent variables and their derivatives (only for approximate symmetries); • solutiondedv: solution of the differential equations specified in diffeqs with respect to leadders; for conditional symmetries, the invariant surface conditions and their needed differential consequences are solved too; • steprelie: integer that stores the status of the computation; 0: no computation done; 1: relieinit() has been called; 2: relieinv() has been called; 3: reliedet() has been called; 4: reliesolve() has been called; • zvar: list of two lists: the first one is the list of the variables arbelem depend on, the second one the corresponding infinitesimals (only for equivalence transformations).

Functions
A short description of the main functions the user may call is as follows: • abelian(gens): checks if the generators gens span an Abelian Lie algebra; • commutatortable(gens): returns the commutator table of the generators gens; • essentialpars(gens,vars): takes a list of generators gens of a multiparameter Lie group of transformations for the variables vars, and returns the generators which are not linearly independent; • generatealgebra(k): once jetorder, xvar and uvar have been properly assigned, this function returns a list of generators spanning the algebra of isometries (for k = 1), affine algebra (for k = 2), projective algebra (for k = 3); • inverselie(k): computes all the minors of order k of the jetorder-th distribution generated by the list of vector fields contained in generators; • liebracket(gen1,gen2): returns the Lie bracket of the generators gen1 and gen2; • newordering(lis,ind): returns a list of the elements in lis reordered according to the permutation ind of the integers {1, . . . , n}, where n is the length of list lis; • nonzerocommutators(gens): returns nzcomm, a list of non-zero commutators of generators gens; • offprintcrack(): prevents reliesolve() to display the steps needed for solving determining equations (this is the default configuration); • onprintcrack(): sets a variable used in Crack package (used in the function reliesolve()) in order to display the steps needed for solving determining equations; • relieclaw(k): returns fluxes, a list of the components of the fluxes of the conservation law corresponding to the k-th Lie generator (obtained after calling to reliegen()); • reliedet(): splits the invariant conditions providing the list deteqs of determining equations; • reliedistrib(): returns the matrix distribution, i.e., a matrix whose rows are the prolonged vector fields evaluated in the list splitsymmetries (computed by the function reliegen(), or by the function generatealgebra(), or suitably assigned by the user; • reliegen(k,lis): returns the list generators; k is an integer (less or equal to the length of symmetries) and lis a list that can be empty; if lis is made by as many values as the number of arbitrary constants occurring in symmetries, generators consists of a list of vector fields, where each vector field is obtained replacing the i-th parameter by the i-th element in the list lis (or 1 if lis is empty or its length is different from the number of arbitrary constants entering symmetries) and the other parameters are replaced by 0; if the list lis is {-1}, then generators is a list with only one element containing the components of the infinitesimals in their general form, i.e., the linear combinations of all admitted generators; the function produces also the list splitsymmetries whose k-th element is a list containing the infinitesimals corresponding to the k-th element of generators; • relieinit(): if input data have been correctly defined, the function initializes the objects for doing the computation; • relieinv(): computes the invariance conditions; returns invcond; • relieprol(): returns the prolongation of a general vector field; • reliesolve(): solves the determining equations for the infinitesimals; returns symmetries; in group classification problems (but also in the case of conditional symmetries), the list symmetries may contain different solutions for the infinitesimals according to the values of freepars; as a default reliesolve() does not display the steps made to obtain the solution of determining equations; the user can see these steps by calling onprintcrack(); this is suggested when reliesolve() seems to use too much time to complete its execution; the list symmetries contains a list of the sets of solutions of determining equations; each element of this list in turn is a list of four elements: the first one is a list of conditions (possibly empty) that remained unsolved; the second one is a list giving the solution to the determining equations, i.e., the expressions of the infinitesimals; the third one is a list containing the parameters involved in the solution; the fourth one is a list of expressions which can not vanish (this list can be empty); • solvable(gens): checks if the generators gens span a solvable Lie algebra; • minors(m,k): m is a matrix and k is a positive integer: returns the list of minors of order k of the matrix m; • nodependence(lis1,lis2): removes the dependence of the objects in the list lis1 upon the variables in the list lis2; • removeelement(lis,elem): removes the element elem from the list lis; • removemultiple(lis1,lis2): removes from the list lis1 the elements of the list lis2; • scalarmult(obj,lis): returns a list whose k-th element is the product of the scalar obj and the k-th element of list lis; • scalarproduct(lis1,lis2): returns the sum of the products element by element of two lists with the same number of elements; • sumlist(lis1,lis2): returns a list summing element by element the two lists with the same length; • zerolist(n): returns a list of n zeros, the empty list if n ≤ 0.

Conclusions
The program ReLie allows the user to perform almost automatically much of the computations needed for investigating ordinary and partial differential equations by means of Lie group methods. The package can also be used in interactive computations where some special assumptions need to be made. Remarkably, the program works in a CAS like Reduce which is open source and freely available for all operating systems. The various routines of ReLie have been tested along the years in many different situations and hopefully are sufficiently reliable. At the url http://mat521.unime.it/oliveri the source code of the package, together with the user's manual, can be found. We illustrated the use of the package by choosing some rather simple examples in different areas of Lie symmetry analysis of differential equations; nevertheless, this does not mean that only simple problems can be faced. In fact, ReLie provides effective also when one has to investigate systems of partial differential equations that require to deal with very long expressions. As a last remark, we observe that this is the first program able to compute approximate Lie symmetries of differential equations containing small terms according to the theory recently proposed [29], as well as approximate Qconditional symmetries [41,43].