If Lie had known about Mathematica:

The author presents the method of first integrals with the application MathLie

by Gerd Baumann

This article discusses a general procedure for solving ordinary differential equations of arbitrary order. The method used is based on symmetries of the given differential equation. The known symmetries allow the derivation of first integrals of the equation. Knowledge of at least r symmetries of an ordinary differential equation of order n with r ≥ n  is the basis for deriving the solution. Our aim is to show that Lie's theory is instrumental in solving an ordinary differential equation of higher-order. We discuss the solution of first order, second order  and higher order equations.

The procedure presented extends the capabilities of Mathematica, especially of DSolve[]. The technique of first integrals is one of several methods for integrating an ODE, provided  the symmetries of the equation are known.

1 Lie's theory and differential equations

The basis for our discussion is a theory conceived  by Sophus. Lie (1842-1899). Lie developed a general theory dealing with symmetries and group properties of differential equations. Lie's theory is a valuable tool for solving ODEs and partial differential equations. In the current text we restrict our considerations to ordinary differential equations. The discussion of partial differential equations can be found in standard books by Olver  [Olver 1986], Bluman and Kumei [Bluman and S. Kumei 1989],  Ibragimov  [Ibragimov 1994-1996] and in a Mathematica package by Baumann [Baumann 1997]. 

The impact of Lie's theory is not restricted to differential equations but is omnipresent in the field of mathematics. Lie deals with fundamental areas of modern mathematics, combining group theory with analysis. A century after Lie, continuous groups are today a central part of modern mathematics. Many newer theories, for example, those on super symmetry or on string theory, are based on aspects of Lie's work.

This paper is organized as follows. Following the introduction, in section two we discuss the concept of symmetry in connection with ordinary differential equations. In section three, we recall some results for first order ODEs. Section four deals with second-order equations. The fifth section gives a generalization of the method to third and higher order equations. The last section summarizes our results.

2 Symmetries

Before we apply Lie´s method to ordinary differential equations, we will briefly discuss symmetries in connection with functions. This section introduces the main concepts of the theory.

One of the remarkable aspects of Lie's theory of symmetry transformations is the invariance of a function under the change of coordinates. When dealing with differential equations, often we find the equation is simplified by an appropriate change of variables. This transformation generally involves both the independent and the dependent variables. In Mathematica, we can represent such transformations by the following list of rules

\(\(rule1 = {x -> Function[{x, u}, X[x, u]], 
		
        u -> Function[{x, u}, U[x, u]]}; \)\)

This kind of transformation—involving the original independent and dependent variables (x,u)—is called a point transformation. A point (x, u) of the manifold 𝔐 is transformed into another point (X, U). Point transformations only take into account changes of the coordinates. The point transformations considered by Lie are transformations depending on at least one parameter ϵ. As we will see, the parameter ϵ is the parameter of the related group. Thus we call ϵ the group parameter of the transformation. A one-parameter transformation is thus given  by

\(rule2 = {x -> Function[{x, u, ϵ}, X[x, u, ϵ]], 

        		u -> Function[{x, u, ϵ}, U[x, u, ϵ]]}\)

Definition:  Symmetry Transformation

Such transformations have the following properties: They are invertible if a corresponding Jacobi determinant exists. Their repeated applications yield transformations of the same kind. An identity of the transformation for ϵ = 0 is given. These three properties are the basis of the definition of a Lie group. The three properties are summarized in the definition of a symmetry transformation.

A set 𝔊 of transformations given, for example by rule2, is a one-parameter group if it contains the identical transformation I = (Subscript[T, 0]) and includes the inverse (Subscript[T, ϵ]^(-1)) and the composition (Subscript[T, ϵ]) ⊗ (Subscript[T, β]) ∈ 𝔊. By choosing the group parameter to be ϵ, the main group property (Subscript[T, ϵ]) ⊗ (Subscript[T, β]) ∈ 𝔊  can be written as

(Subscript[T, ϵ]) ⊗ (Subscript[T, β]) = (Subscript[T, ϵ + β])				(1)

that is

X(X(x,u,ϵ),U(x,u,e),β)  =  X(x,u,ϵ+β)		(2)

U(X(x,u,ϵ),U(x,u,ϵ),β)  ==  U(x,u,ϵ+β).	(3)

In special applications, conditions (2) and (3) only hold for sufficiently small values of ϵ and β. In such cases, we arrive at what is called a local one parameter group 𝔊 or an infinitesimal group.

These properties ensure that the transformations given in rule2 form a one-parameter group of point transformations.

A favorite example of Lie is the one-parameter group of rotations in the (x, u) plane. This type of point transformation is represented in Mathematica  by the rules

rotation  = 
       {X -> 
          Function[{x, u, ϵ}, 
            x Cos[ϵ] - u Sin[ϵ]], 
        U -> 
          Function[{x, u, ϵ}, 
            x Sin[ϵ] + u Cos[ϵ]]};

The one-parameter symmetry transformation defined in rotation can be graphically displayed. By changing the parameter ϵ and plotting the function, we can animate the transformation. In general, the image of an initial point will move along some curve. The change in ϵ moves an initial point ((Subscript[x, 0]), (Subscript[u, 0])) to another point ((Subscript[x, 1]), (Subscript[u, 1])) along the curve. In the present case, the point moves along a circle. The initial point ((Subscript[x, 0]), (Subscript[u, 0])) follows from rotation if we set ϵ = 0. Thus, the motion of the point occurs when ϵ changes. Mathematica offers us the ability to show the motion of the point by an animation of sequential pictures. The animation can be started with Ctrl-Y.

\(\(	\(Table[
      ParametricPlot[
        Evaluate[
	
          \({X[x, u, ϵ], U[x, u, ϵ]} /. rotation\) /. {
              x -> 1, u -> 1}], 
		{ϵ, 0, end}, 
	
        AspectRatio -> Automatic, 
        PlotRange -> {{\(-1.5\), 1.5}, {\(-1.5\), 1.5}}, 
        PlotStyle -> RGBColor[0, 0, 1], 
		
        AxesLabel -> {"x", "u"}], 
      	{end, N[2  π\/20], N[2  π], N[2  π\/20]}]; \)\)\)

The above animation shows that an initial point ((Subscript[x, 0]), (Subscript[u, 0])) moves along a circle if the group parameter ϵ changes continuously. If we repeat this transformation for different initial points ((Subscript[x, 0]), (Subscript[u, 0])), we obtain a picture representing the global action of the transformation in the (x, u) plane. The figure below contains the orbits of the transformation given in rotation for several initial points. The figure is created with initial points located along the x-axis.

\(\(ParametricPlot[
      Evaluate[
		
        Table[\({X[x, u, e], U[x, u, e]} /. rotation\) /. {x -> xo, 
              u -> 1}, 
			{xo, 1, 5, 1\/2}]], {e, 0, 2\ π}, 
	
      AspectRatio -> Automatic, 
      PlotStyle -> Table[Hue[x], 
			{x, 0, 1, 1\/5}], 
      AxesLabel -> {"x", "u"}]; \)\)

Figure 1

The circles in different colors represent the point transformation of rotation. Each color corresponds to an initial point located along the x-axis. All circles represent the global transformation of the manifold under rotation.

The above figure shows that the given transformation will move a point along a circle if ϵ is changed. Each curve represents points that can be transformed into each other by the given transformation. This is a typical example of a global representation of symmetry transformations.

2. Infinitesimal transformations

One of the most important aspects of Lie's theory is the discovery that a symmetry transformation as given above can be represented in a much simpler form. Lie called this simplified representation an infinitesimal transformation. The main idea behind this simplification is the existence of an identity transformation. An identity transformation denotes the lowest approximation of a transformation around the identity. Correspondingly, the finite transformation can be expanded into a Taylor series around the identity transformation. In Mathematica we have direct access to such an expansion when applying the function Series[] to the global representation. 

The following is an example for a general point transformation in the (x, u) plane. We derive the infinitesimal representation by expanding the one-parameter transformation around ϵ = 0 up to the first order in ϵ:

\(infinitesimalTrafo = 
	
      Series[{X[x, u, ϵ], U[x, u, ϵ]}, {ϵ, 0, 1}]
        \)

If we use the group properties of the transformation, especially the identity X(x, u, ϵ = 0) = x and U(x, u, ϵ = 0) = u, we can reduce the expression to

\(infinitesimalTrafo = 
      MatrixForm[
		
        infinitesimalTrafo /. 
          \(({∂\_ϵ X[x, u, ϵ] -> 
                  ξ[x, u], 
					
                ∂\_ϵ U[x, u, ϵ] -> 
                  ϕ[x, u], 
					X[x, u, 0] -> x, 
					
                U[x, u, 0] -> u} /. ϵ -> 0)\)]\)

The result shows that the transformation is represented by the original variables plus linear terms in ϵ. The coefficients of ϵ are called the infinitesimals of the transformation and are denoted by ξ and ϕ. Both expressions are the expansion coefficients of the Taylor expansion \(TraditionalForm
      \`ξ = 
        ∂\_ϵ X  
            |\_\(ϵ = 0\)\), and \(TraditionalForm
      \`ϕ = 
        ∂\_ϵ U  
            |\_\(ϵ = 0\)\), respectively.

The above is summarized in Lie's first theorem. The theorem considers the inverse problem of an infinitesimal representation. Since, in this case the infinitesimals ξ and ϕ are known, we are looking for a global transformation.

Theorem:  Lie´s first theorem

Let us assume that there exists a parameterization of a transformation such that the global transformation like rotation is equivalent to the solution of the initial value problem for the system of first order differential equations

\(TraditionalForm
      \`dX\/dϵ\  = \ ξ(X(ϵ), U(ϵ))\)		(4)

\(TraditionalForm
      \`dU\/dϵ\  = \ ϕ(X(ϵ), U(ϵ))\)		(5)

with the initial conditions

X(0)  =  x and U(0)  =  u.	(6)

Lie's theorem states that the finite or global transformation represented by X and U is derived from an initial value problem. By integrating these equations with respect to ϵ and eliminating the group parameter, we end up with a global transformation.

We assume we know the infinitesimals ξ and ϕ to which we apply Lie's first  theorem in a straightforward way. Lie's first theorem can be implemented in a Mathematica function allowing the determination of the global symmetry transformation. All steps occurring in the theorem are combined by this function. The function GlobalSymmetry-Transformation[]  calculates the finite symmetry transformation if the infinitesimals of the transformation are known. The input variables for this function are the infinitesimals ξ and ϕ. We also need the dependent and independent variables for our calculation. The function is designed to calculate a general representation of the global transformation for an arbitrary number of independent and dependent variables.

\(GlobalSymmetryTransformation[xi_List, phi_List, 
	depend_List, 
        independ_List] := 
      Block[
		{variables, Vars, infini, dVars, equations, 
			
          initial, sol}, 
		variables = Join[depend, independ]; 
		
        infini = Join[phi, xi]; 
		
        Vars = ToExpression/@\(ToUpperCase/@\(ToString/@variables\)\); 
		
        dVars = Table[
            Vars[[i]]@@{ϵ}, {i, 
              1, Length[Vars]}]; 
		
        equations = 
          Table[∂\_ϵ dVars[[i
                    ]] == 
                infini[[i]], 

              				{i, 1, Length[Vars]}] /. 
            Thread[variables -> dVars]; 
		
        initial = Thread[dVars == variables] /. ϵ -> 0; 
		
        equations = Join[equations, initial]; 
		
        sol = DSolve[equations, Vars, ϵ]]\)

The function can be verified with the known symmetry transformation of rotation. The infinitesimal form of this symmetry transformation  is given by  ξ  = -u and ϕ = x. Applying the function GlobalSymmetryTransformation[] on these infinitesimals, we find

\(GlobalSymmetryTransformation[{\(-u\)}, {x}, 
	{u}, {x}]\)

The result is identical with the transformation given by rotation, realizing a rotation of an initial point ((Subscript[x, 0]), (Subscript[u, 0])) around the origin.

So far we have demonstrated that a global transformation is equivalent to its infinitesimal representation. However, the question of how symmetries and transformations are connected remains. Their link is the subject of the following subsection.

2. Group invariants

In this subsection a criterion of invariance is formulated. This criterion is instrumental for solving ordinary differential equations. We will show here that a symmetry transformation does not change the form of the differential equation and thus opens a way toward constructing the solution of the equation.

We know that a one-parameter group of transformations is represented by infinitesimal transformations. As the transformation is applied to a point of the manifold 𝔐, it moves along a path in such a way that the path maps into itself. This basic concept of invariance needs a mathematical formulation. For the moment we restrict our considerations to curves in the plane.

Group invariants are the basic quantities of a symmetry analysis. A group invariant is defined as follows

Definition: Group invariant

A function F(x, u) is an invariant of the group of transformations, if

F(X(x,u,ϵ),U(x,u,ϵ)) = F(x,u)		(7)

is satisfied identically in x,  u for all values of the group parameter ϵ.

This definition contains the basic information for an invariant in symmetry analysis in a nut shell. The definition states that a function is an invariant if it has the same representation in the original and in the new coordinates. The definition also shows us a way how the invariant F may be calculated.

We know that any curve in the (x, u)-plane can be represented by a function F(x, u) = constant. Following our definition this curve is considered to be invariant under a transformation if the curve remains the same in both coordinate systems. As already mentioned, global transformations need only an infinitesimal representation such as

\(\(infini = {
        X -> Function[{x, u, ϵ}, x + ϵ\ ξ[x, u]], 
        
		U -> 
          Function[{x, u, ϵ}, u + ϵ\ ϕ[x, u]]}; \)\)

where ξ and ϕ are the infinitesimals. The analytical expression for the invariance is simply stated by the relation

invar = F[x, u] == F[X[x, u, ϵ], U[x, u, ϵ]]

showing that the form of the curve is the same in the old and new coordinates. Applying the transformation to this expression delivers the invariant condition in an implicit form

ivarI = invar /. infini

The explicit representation of the invariance condition is obtained by a Taylor expansion of the equation of invariance with respect to parameter ϵ.

Thread[Series[ivarI, {ϵ, 0, 1}], Equal]

The invariance condition is satisfied once all terms containing the infinitesimal parameter ϵ have vanished. Introducing Lie's symbol, \(TraditionalForm
      \`v\&⇀\  = \ 
        ξ ∂\_x\(+ϕ\) ∂\_u\), today known as the tangent vector field or group generator, we can write the invariance condition as

\(TraditionalForm
      \`F(x, u)\  = \ 
        F(x, u)\  + \ v\&⇀\ \(F(x, u)\)\ ϵ\  + \ 
          o(ϵ\^2)\).		(8)

This calculation shows that a sufficient condition of invariance is (Overscript[v, ⇀]) F =  0. This condition is accessible by an algorithmic procedure and thus applicable to any function F.

Assuming that the tangent vector field applied to the curve F has to vanish, we can immediately conclude that all higher terms in  (8) containing multiple applications of the vector field (Overscript[(v ), ⇀]) also have to vanish. As a result, we have a sufficient condition of invariance. This condition is central to the determination of symmetries and can be generalized to a manifold of q independent variables  x = ((Subscript[x, 1]), (Subscript[x, 2]), …, (Subscript[x, q])) and p dependent variables  u = ((u^1), (u^2), …, (u^p)). In the q × p-dimensional space, the invariance condition reads

\(TraditionalForm
      \`∑\+\(i = 1\)\%q ξ\_i\ 
              \(∂F \((x, u)\)\)\/∂x\_i\  + \ 
          ∑\+\(α = 1\)\%p ϕ\_α\ 
              \(∂F \((x, u)\)\)\/∂u\^α\  = \ 0
        \)			(9)

where (Subscript[ξ, i]) and (Subscript[ϕ, α]) are the unknown infinitesimals of a point transformation.

The term of invariance boils down to the point  that each  one-parameter group of point transformations in the plane has one independent invariant which can be taken to be the left-hand side of any first integral J(x, u)  =  C of the characteristic equation (9) of the tangent vector field.

Let us consider the example of a circle, obviously invariant under the symmetry group of rotations. As we know from (9) and Lie's first theorem the tangent vector field of corresponding group of rotation is given by the operator

Trot[f_] := -u (∂_x) f + x (∂_u) f

A circle in the coordinates x and u can be represented by

circ = (x^2) + (u^2) - const;

The application of the operator Trot[] to the circle  gives us

Trot[circ]

The result demonstrates that the tangent vector field applied to a curve invariant under the given symmetry vanishes. Thus each circle is mapped into itself. This condition is very useful in the determination of symmetries not only for curves but also for differential equations.

2. Symmetries of ordinary differential equations

So far we have introduced symmetry transformations of curves. The essential difference between curves and ODEs is that the manifold 𝔐 is extended or prolonged by the derivatives. The manifold  of an ODE does not only contain the independent and dependent variable but also the derivative u' = du/dx  =  p. Thus we have in the simplest case a three-dimensional space with coordinates (x, u, p). We are now looking for the symmetries of this kind of manifold. Again we can define the invariance of a differential equation in a similar way as for functions. Let us  state the definition for the most general case of an nth order ODE.

Definition: Symmetry group of a differential equation

A group 𝔊 of point transformations for variables in the manifold 𝔐 is a symmetry group of an nth-order ordinary differential equation

\(TraditionalForm
      \`Δ(x, u(x), du\/dx, …, \(\(d\^n\) u\)\/dx\^n)
          \  = \ 0\)		(10)

admitting 𝔊 if the nth extended manifold (𝔐^n) is invariant with respect to the nth prolongation (Subscript[𝔊, n]) of the group 𝔊.

This definition contains the special case for first-order ordinary differential equations.  A first-order differential equation

Δ(x, u(x), u' (x))  =  0 	(11)

admits a group 𝔊 if the once extended manifold (𝔐^1), the surface in the space x, u, u', is invariant with respect to the first prolongation (Subscript[𝔊, 1]) of 𝔊. This means that equation (11) is invariant under the coordinate transformations X = X(x, u, ϵ) and U = U(x, u, ϵ). The invariance condition can be formulated as

\(TraditionalForm
      \`Δ(X, U, U')\  = \ Δ(x, u, u')\).		(12)

This kind of relation is also true for the general case of an nth-order ODE. Consequently, by knowing the transformations X and U, we also know the symmetries of the equation. Lie demonstrated that the knowledge of the infinitesimal transformations is sufficient to solve the equations. We will  state here only the fact that the knowledge of the infinitesimals is sufficient to solve the equations. We do not show how the infinitesimals can be calculated with Mathematica. The interested reader may consult the function Lie[] published by Baumann [MathSource 0202-622]. Once we know the symmetries, it is quite easy to derive the solutions  of the equations. The following sections will illustrate solutions for first-order, second-order, and higher-order equations.

3 First Order ODEs

The general form of a first order ordinary differential equation is given by

F(x, u(x), (du/dx))  =  0		(13)

where F is an arbitrary function combining  the independent and the dependent variables x, u, and the derivatives  du/dx  in a general way. Our aim is to use point symmetries for the solution. Knowing the infinitesimal transformation of equation (13) allows two basic strategies for a solution. The first strategy is based on canonical variables extensively discussed by Ibragimov  [Ibragimov 1994-1996]; Bluman and Kumei [Bluman and S. Kumei 1989]; as well as Olver[Olver 1986]. The second strategy, originally invented by Lie [Lie 1874] uses the well-known procedure of an integrating factor. Here we restrict our considerations to the second strategy.

The common belief in literature is that the method involving an integrating factor is useful only in connection with a first-order ordinary differential equation. In contrast, Lie already developed a generalization of the procedure in the late 1870s.

In the previous section, we demonstrated that a curve is invariant under a symmetry transformation if the tangent vector applied to the curve vanishes. Let us here assume that the curves discussed are integral curves or solutions of a first-order differential equation. In this case, we can state the following theorem

Theorem: Integrating Factor

If the first-order ordinary differential equation

(du/dx) - ω(x, u) = 0		(14)

allows a one-parameter group given by the vector field (Overscript[v, ⇀]) = ξ (∂_x)+ ϕ (∂_u) , then and only then the following function

μ  =  (1/(ξ ω  - ϕ))			(15)

is an integrating factor with ξω - ϕ ≠ 0. If this is the case, the original equation is solved by a quadrature

\(TraditionalForm
      \`∫\((ωdx - du)\)\ μ\  = \ \(constant . \)\)		(16)

The theorem states that the knowledge of the infinitesimals ξ and ϕ allows the derivation of a solution just by integration. According to Lie, relation (16) can be symbolically simplified if we introduce two determinants of the form

dΔ = det(({{dx, du}, {1, ω(x, u)}}))	(17)

and

Δ = det(({{ξ, ϕ}, {1, ω(x, u)}})) = (1/μ).		(18)

Then relation (16) reduces to the simple relation

∫ (d(Δ/Δ)) = constant . 		(19)

Equation (19) combines all information necessary for our solution in a nut shell. All we need to know are the infinitesimals and the right-hand side of the equation (14).  Integration over the manifold provides us with the solution.

As a first example of integration, we discuss an ODE which Ibragimov [Ibragimov 1994-1996] used to demonstrate the procedure of canonical transformations. The equation under consideration is a Riccati equation of the form

riccati  =  (∂_x) u[x] + u[x]^2 - (2/(x^2)) == 0;

The infinitesimal symmetries of this equation are scalings of the coordinates given by ξ = x and ϕ = -u. These two expressions and the right hand side of the equation are the main components for deriving the Δ matrix of equation (19). With the help of the MathLie package  we are able to create the Lie-matrix by using the independent and the dependent variable, the right hand side of the equation ω, and the infinitesimals ξ and ϕ.

Δmatrix = 
      DeltaMatrix[x, u, -u[x]^2 + (2/(x^2)), 1, {{{x}, {-u[x]}}}]; 
    TableForm[Δmatrix]

The above matrix is the basis for deriving a first integral. The integral is the result of a series of integrations along a path in the (x, u)-plane. FirstIntegral[] of MathLie carries out the calculation

int1  =  FirstIntegral[x, u, Δmatrix, 1] ==  c1

The input parameters for function FirstIntegral[] are the independent and dependent variables, the Lie matrix, and the number of the first integral. The number of first integrals is given by the order of the ODE. In case of first-order ODEs there exists only one integral. For second-order ODEs there exist two integrals and for nth-order ODEs there may exist n first integrals. If we solve the above expression with respect to the unknown function u, we end up with the solution of the problem.

sriccati  =  Solve[int1, u[x]]

Here c1 is an arbitrary constant of integration. We arrive at the same solution when applying  DSolve[] to the problem.

DSolve[riccati, u, x]

This simple example was solved by Mathematica. However, things change if we examine higher-order equations. The following section gives a generalization of the integrating factor method for a second-order ODE.

4 Second Order ODEs

Second-order ordinary differential equations are very important for applications in physics and engineering. All equations based on Newton's second law are second-order equations. Thus, mechanics, for example, is mainly based on second-order ordinary differential equations. The integration theory of second-order ODEs was developed by Lie during the years 1871-1874. In 1891 Lie's theory of integration was published by Lie and Scheffers in Vorlesungen über Differentialgleichungen mit bekannten infinitesimalen Transformationen. Scheffers, a student of Lie,  assembled all of Lie's work on ordinary differential equations in a beginner's book. In this book Lie described the problem  of integrating a differential equation as follows:

"I observed that a large number of ordinary differential equations integrated by older integration methods are invariant under easyly derivable classes of transformations. The older integration methods are all based on the transformation properties of the equation. In other words, I realized that the term differential invariant of a finite continuous group is contained implicitly in every textbook on ordinary differential equations. Discovering the connection between transformation groups and older integration strategies I started to develop a general integration theory based on the finite or infinitesimal transformation of the equation. In my investigations it was clear from the beginning that the related transformations always created a group for each case.

I noticed that the majority of ordinary differential equations which were integrable by the old methods were left invariant under certain transformations, and that these integration methods consisted in using that property. Once I had thus represented many old integration methods from a common viewpoint, I set myself the natural problem: to develop a general theory of integration for all ordinary differential equations admitting finite or infinitesimal transformations." [Lie 1874].

We will exemplify Lie's line of thought below. The most general form of a second-order ODE  is given by

F(x, u, u', u'')  =  0. 		(20)

For our purposes, we assume that equation (20) is solvable with respect to the second-order derivative. Thus we consider equations in the form

u''  =  ω(x, u, u')		(21)

where ω is a given function of x, u, and u'.

According to Lie we can now state: If a second-order equation admits a finite symmetry of dimension r ≥ 2, it can be integrated by a group-theoretic quadrature method. This can be done in various ways, one of which is given by the following algorithm:

1. compute the Lie algebra (Subscript[L, r]). A basis of  (Subscript[L, r]) is the set  (Subscript[(Overscript[v, ⇀]), 1]),  (Subscript[(Overscript[v, ⇀]), 2]), …, (Subscript[(Overscript[v, ⇀]), r]). The tangent vector fields (Subscript[(Overscript[v, ⇀]), i]) follow by appropriately specifying the group constants.

2. If r = 2, go to the next step;
     if r > 2,  distinguish any two-dimensional subalgebra (Subscript[L, 2]) of (Subscript[L, r]).
     If r = 1, the order of the equation may be lowered;
     if r = 0, the group method is not useful.

3. Calculate the Lie determinants d(Subscript[Δ, i]) and Δ and determine the two first integrals by integration. The Lie determinants are defined by

d(Subscript[Δ, 1]) = det(({{dx, du, du'}, {(Subscript[ξ, 2]), (Subscript[ϕ, 2]), \(ϕ\_2\%'\)}, {1, u', ω}})),         d(Subscript[Δ, 2])= det(({{(Subscript[ξ, 1]), (Subscript[ϕ, 1]), (Subscript[ϕ, 1])'}, {dx, du, \(du\^'\)}, {1, u', ω}}))		(22)

and

Δ = det(({{(Subscript[ξ, 1]), (Subscript[ϕ, 1]), \(ϕ\_1\%'\)}, {(Subscript[ξ, 2]), (Subscript[ϕ, 2]), \(ϕ\_2\%'\)}, {1, u', ω}}))		(23)

where (Subscript[ξ, i]) and (Subscript[ϕ, i]) are the infinitesimals corresponding to the vector field (Subscript[(Overscript[v, ⇀]), i]), and \(TraditionalForm\`\(ϕ\_i\%'\)\) denote the first extensions of the infinitesimals (Subscript[ϕ, i]). The first integrals (Subscript[ψ, i]) related to the Lie determinants d(Subscript[Δ, i]) are given by

\(TraditionalForm
      \`ψ\_i\  = \ 
        \(∫dΔ\_i\/Δ\  = \ c\_i\), 
      \ \ \ \ \ \ \ \ \ i = 1, 2. \)		(24)

The constants (Subscript[c, i]) denote the integration constants of the ODE.

4. Solve one of the two integrals with respect to the first order derivative and substitute the result into the remaining relation.

5. If we can solve the resulting relation with respect to the unknown function u, we end up with an explicit solution. Otherwise, we get the solution in an implicit form.

These five steps are completely algorithmic and are implemented among others in the package MathLie. The functions of MathLie carry out the needed calculations completely automatically. To show how the solution procedure works interactively, we demonstrate the algorithm by two  examples.

Example: 1

The first example considers a nonlinear second order equation in which the nonlinearity is given by the square of the first derivative. This nonlinear term is multiplied by a real constant α. The second order ODE published by Kamke [Kamke 1977] as equation 6.1 and 6.2 follow from our equation by specifying α = 1 or α = 6, respectively. The generalized equation we examine reads

\(firstExample\  = \ 
      ∂\_{x, 2}u[x] - α\ \((∂\_x u[x])\)\^2\)

The infinitesimal symmetries are derived by the function Infinitesimals[] which is part of the package MathLie.

\(infi\  = \ 
      Infinitesimals[firstExample, u, x, {α}, 
        SubstitutionRules -> {∂\_{x, 2}u[x]}]\)

The result is a symmetry group with 8 group parameters ki, i=1,2,…8. The infinitesimals ξ = xi[1] and ϕ =  phi[1] are represented in a pure function form.  The group constants ki characterize the symmetries of the equation. We note that this equation allows a transformation which reduces the original equation to the simple form y''  = 0.  This reduction is always possible if a second order equation allows a symmetry group of order eight (cf. [Lie 1874].

To solve the nonlinear second order ODE, we select a solvable subgroup of order two from the eight-dimensional group. The last argument of the function DeltaMatrix[] contains the groups related to the solvable algebra. We choose  subgroups k8=β and k2=1.

inf1  =  
      {{xi[1][x, u]}, {phi[1][x, u]}}  /. infi /. {
            k1 -> 0, k2 -> 0, k3 -> 0, k4 -> 0, 
            k5 -> 0, k6 -> 0, k7 -> 0, k8 -> β}
              /. u -> u[x];

inf2  =  
      {{xi[1][x, u]}, {phi[1][x, u]}}  /. infi /. {
            k1 -> 0, k2 -> 1, k3 -> 0, k4 -> 0, 
            k5 -> 0, k6 -> 0, k7 -> 0, k8 -> 0}  /. 
        u -> u[x];

Lie's matrix replacing the integrating factor is now calculated by using the function DeltaMatrix[]

Δmatrix = 
      DeltaMatrix[x, u, α ((((∂_x) u[x]))^2), 
        2, {inf1, inf2}]; TableForm[Δmatrix]

The result is a 3×3 matrix containing the infinitesimals and the first prolongation of the two subgroups and the right hand side of the equation. The determinant of this matrix

Det[Δmatrix]

is a non-vanishing expression containing a second order polynomial in the first derivatives. The coefficients of this polynomial depends on the independent variable x and the parameters α and β. Inserting the Lie matrix  into equation (24), we are able to calculate the first integrals of the nonlinear ODE

integrals  =  
      Thread[FirstIntegral[x, u, Δmatrix, {1, 2}] ==  {c1, c2}]

The result contains two expressions for the integrals combining the variable u and its first order derivative in an algebraic way. These two first integrals define two manifolds in the space (x, u, u'). The projection of the intersection onto the (x, u)-plane of these two manifolds defines the solution we are looking for. The following figure represents a case with fixed values c1 and c2.

Figure 2

Figure 2 shows the intersection of  the two integrals  for c1=c2=1. The parameters α and β  take the values α = 1/10 and β=1. The intersecting line represents the solution of the equation if we project the intersection to the (x, u)-plane.

Since we know two first integrals, we are able to eliminate analytically the derivatives from the two integrals. We solve the first relation with respect to u'

sol1  =  
      Solve[integrals[[1]], 
          u'[x]] // Simplify

Substituting this result into both expressions for the integrals gives us

int21  =  
      integrals  /. sol1[[1, 1]] // 
          PowerExpand // Simplify

The resulting list contains True and an implicit representation of the solution. The explicit solution follows from the second relation if we solve it with respect to u.

solution  =  
      Solve[int21[[2]], u[x]]

The constants c1 and c2 denote integration constants. α and β are the model constant and the group parameter, respectively. The occurrence of the group parameter β as a multiplier  reminds us that the additive integration constant is closely related to the group of translations.

The same solution of the original ODE is derived by the function DSolve[] in the form

DSolve[firstExample == 0, u, x]

This example shows that the solution steps of Lie's method of first integrals results into the same solution as  that of Mathematica.

Example: 2

The second example for a second-order ODE is related to the problem of a generalized suspended cable equation. The problem of the suspended cable is discussed by Ames [Ames 1968]. We examine a generalization of the cable equation with an arbitrary exponent in the second term. This problem is similar to the first example. We only add α to the square of the derivative and rise it to the νth power. These small changes lead to substantial changes in the solution.

secondExample  =  
      (∂_(x, x))u[x]  - 
          α (((1 +  ((((∂_x) u[x]))^2)))^ν)  == 0

The original cable equation follows from our equation with ν = 1/2. We first calculate the point symmetries of this equation for arbitrary ν.

The infinitesimals follow by applying the MathLie function Infinitesimals[] to the equation.

\(infi\  = \ 
      Infinitesimals[secondExample, u, x, {α, ν}, 
        SubstitutionRules -> {∂\_{x, 2}u[x]}]\)

The result is a two-dimensional symmetry transformation representing translations in the independent and dependent variables. The first difference compared to example 1 is that the group is much smaller. This reduction of the group order has consequences for the solutions.

The two subgroups necessary for integration are derived by setting each of the group parameters k1 or k2 independently to unity. Since the symmetry group is two-dimensional, we know from Lie that the algebra is always solvable. In the worst case, the solution  may be represented in the end by an implicit representation. The infinitesimals for the two subgroups follow from

inf1  =  
      {{xi[1][x, u]}, {phi[1][x, u]}}  /. infi /. {
            k1 -> 1, k2 -> 0}  /. u -> u[x];

inf2  =  
      {{xi[1][x, u]}, {phi[1][x, u]}}  /. infi /. {
            k1 -> 0, k2 -> 1}  /. u -> u[x];

The related Lie matrix for the case with ν = 3 is calculated by

Δmatrix = 
      DeltaMatrix[x, u, 
        α (((1 +  ((((∂_x) u[x]))^2)))^3) , 
        2, {inf1, inf2}]; TableForm[Δmatrix]

The determinant of this matrix is given by a polynomial of sixth order in u'

Det[Δmatrix]

The first integrals for the generalized cable equation follow from

integrals  =  
      Thread[FirstIntegral[x, u, Δmatrix, {1, 2}] ==  {c1, c2}]

Solving the first relation with respect to u' allows us to eliminate this term in the second integral.

sol1  =  
      Solve[integrals[[1]], 
          u'[x]] // Simplify

The solution for u' consists of four expressions containing square roots of u. The differences of the four solutions are the signs in front of the first and second term. Inserting, for example, the second solution into the second first integral, we get one of the final solutions in an implicit representation

int21  =  
      integrals[[2]]  /. 
        sol1[[2, 1]]

The explicit solution of this expression is impossible since it contains transcendental functions.

solution  =  Solve[int21, u[x]]

At this point we have to accept that the solution can  only be represented in an implicit form. Trying to solve the original equation by Mathematica, we learn that DSolve[] is not capable to resolve the relation for the first integral. A glance at the result of DSolve[] explains the reason why.

DSolve[secondExample  /.  ν -> 3, u, x]

The above shows that a transcendental function is given in an essential non-algebraic way. The steps presented above demonstrate that the solution of the generalized cable equation is solvable in an implicit form.

The two examples show that the technique of the integrating factor can be generalized to second order equations. We also observe that the presented procedure is capable of deriving solutions for cases in which Mathematica fails. The following section will show you the extension of the integrating factor technique to a general nth order ODE. Lie called this procedure the method of generalized multipliers.

5 Higher Order ODEs

Differential equations of higher order arise naturally in physics. For example, third-order ODEs show up in fluid dynamics and fourth-order ODEs in elasticity. For general higher order equations there exist hardly any techniques for obtaining explicit symbolic solutions. Higher order ODEs are thus not well studied in the literature. In the following, we will present a symbolic technique for producing explicit solutions independent of the order of the equation. The method described in the preceding sections can, without essential changes, be applied to the solution of differential equations of higher order. The essential change in the method is the extension of the Lie matrix to higher prolongations. Extending the Lie matrix to higher prolongations is the key step for generalizing the procedure of integrating factors. Lie called this extension the determination of the multiplier of the differential equation.

The five steps of integration discussed for second order equations remain the same for higher order equations. However, the dimension of the Lie matrix changes from a 3×3 matrix to a (n + 1)×(n + 1) matrix. Before we describe the algorithm to solve higher order ODEs let us first state the general settings for higher order equations.

The most general form of an nth-order ODE is given by

\(TraditionalForm
      \`F(x, u, u', \(u'\)', …, u\^\((n)\))\  = \ 0. \)		(25)

where (u^((n))) =  (d^n) u/(dx^n) denotes the nth derivative of u with respect to x. In the following, we assume that equation (22) can be solved with respect to the nth derivative. The actual equation under consideration is thus

\(TraditionalForm
      \`u\^\((n)\)\  = \ ω(x, u, u', …, u\^\((n - 1)\))\)		(26)

where ω is a given function of x, u, u', …, (u^((n - 1))).

If an nth-order equation admits a finite symmetry of dimension r ≥ n, then the equation can be integrated by the group theoretic quadrature method. This again can be done in various ways. For a detailed discussion of  other procedures, compare parts II and III of this series of papers. One of the group theoretic algorithms is based on first integrals. The algorithm is summarized as follows:

1. compute the Lie algebra (Subscript[L, r]). A basis for  (Subscript[L, r]) is the set  (Subscript[(Overscript[v, ⇀]), 1]),  (Subscript[(Overscript[v, ⇀]), 2]), …, (Subscript[(Overscript[v, ⇀]), r]). The tangent vector fields (Subscript[(Overscript[v, ⇀]), i]) follow from appropriately specifying the group constants.

2. If r = n, go to the next step;
     if r > n,  distinguish any n-dimensional subalgebra (Subscript[L, n]) of (Subscript[L, r]).
     If r = n - 1, the order of the equation may be lowered;
     if r = 0, the group method is not useful.

3. Calculate the Lie determinants d(Subscript[Δ, i]) and Δ and if  possible determine the n first integrals by integration. The Lie determinants are defined by

d(Subscript[Δ, 1]) = det(({{dx, du, du', …, (du^((n - 1)))}, {(Subscript[ξ, 2]), (Subscript[ϕ, 2]), \(ϕ\_2\%'\), …, (Subscript[ϕ, 2]^((n - 1)))}, {⋮, ⋮, ⋮, ⋱, ⋮}, {1, u', u'', …, ω}})),        
 	(27)

and

d(Subscript[Δ, i])= det(({{(Subscript[ξ, 1]), (Subscript[ϕ, 1]), (Subscript[ϕ, 1])', …, (Subscript[ϕ, 1]^((n - 1)))}, {⋮, ⋮, ⋮, ⋱, ⋮}, {dx, du, \(du\^'\), …, (du^((n - 1)))}, {⋮, ⋮, ⋮, ⋱, ⋮}, {1, u', u'', ⋯, ω}}))	(28)

where i denotes the row of the Lie matrix in which infinitesimals are replaced by differentials. The n×n Lie determinant  reads

Δ = det(({{(Subscript[ξ, 1]), (Subscript[ϕ, 1]), \(ϕ\_1\%'\), …, (Subscript[ϕ, 1]^((n - 1)))}, {(Subscript[ξ, 2]), (Subscript[ϕ, 2]), \(ϕ\_2\%'\), …, (Subscript[ϕ, 2]^((n - 1)))}, {⋮, ⋮, ⋮, ⋱, ⋮}, {1, u', u'', …, ω}})) 		 (29)

where (Subscript[ξ, i]) and (Subscript[ϕ, i]) are the infinitesimals of the vector field (Subscript[(Overscript[v, ⇀]), i]) and \(TraditionalForm\`\(ϕ\_i\%'\)\) denote the first and (Subscript[ϕ, i]^((n - 1))) the (n-1)th extension of the infinitesimals (Subscript[ϕ, i]). The corresponding first integrals (Subscript[ψ, i]) are given by

\(TraditionalForm
      \`ψ\_i\  = \ 
        \(∫dΔ\_i\/Δ\  = \ c\_i\), 
      \ \ \ \ \ \ \ \ \ i = 1, 2, …, \(n . \)\)		(30)

The constants (Subscript[c, i]) denote the integration constants of the ODE.

4. Solve one of the n integrals with respect to the (n-1)th order derivative and substitute the result into the remaining relations. Repeat this procedure until no derivative remains in the relations.

5. If we can solve the resulting relation with respect to the unknown function u, we have found an explicit solution. Otherwise, our solution is implicit.

This procedure becomes very cumbersome with increasing orders of the equation if done by hand. In principle, the procedure can be  applied to any kind of linear or nonlinear ODE. How the algorithm works in particular examples is demonstrated below.

Example: 1

The first example is a third-order equation listed by Kamke [Kamke 1977] as No. 7.13.

\(\(thirdOrderExample\  = \ 
      ∂\_\(x, x\)u[x] ∂\_{x, 3}u[x] - 
        α 
          \@\( 1 + \(β\^2\) \((∂\_\(x, x\)u[x])\)\^2\)\ 
    \)\)

The parameters α and β are real constants. According to Kamke this third-order ODE is solvable and the solution can only be represented in parametric form. We will show that an explicit solution of the equation is possible. First let us check if Mathematica can solve the third-order equation.

DSolve[thirdOrderExample == 0, u, x]

The above line shows that Mathematica is not capable of solving the equation. The question now is can we derive the necessary number of symmetries in order to integrate the equation? The derivation of the symmetries is the first step in the general algorithm. The calculation of symmetries is carried out by the MathLie function Infinitesimals[].

\(infi\  = \ 
      Infinitesimals[thirdOrderExample, u, x, {α, β}, 
        SubstitutionRules -> {∂\_{x, 3}u[x]}]\)

The result is a symmetry group of order three. The specific symmetries are denoted by the group constants k1, k2, and k3. Each of these parameters is related to a vector field (Subscript[(Overscript[v, ⇀]), i]), i = 1, 2, 3. Since the number of vector fields is equal to the order of  equation, we can go to step three of the algorithm. In the third step we determine the Lie matrix by inserting the prolongations of the infinitesimals and the equation itself

inf1  =  
      {{xi[1][x, u]}, {phi[1][x, u]}}  /. infi /. {
            k1 -> 1, k2 -> 0, k3 -> 0}  /. u -> u[x];

inf2  =  
      {{xi[1][x, u]}, {phi[1][x, u]}}  /. infi  /. {
            k1 -> 0, k2 -> 1, k3 -> 0}  /. u -> u[x];

inf3  =  
      {{xi[1][x, u]}, {phi[1][x, u]}}  /. infi  /. {
            k1 -> 0, k2 -> 0, k3 -> 1}  /. u -> u[x];

The Lie matrix is derived by

Δmatrix = 
      DeltaMatrix[x, u, 
        ((α 
            (Sqrt[ 1 + 
                (β^2) 
                  ((((∂_(x, x))u[x]))^2)]))/(∂_(x, x
              )))u[x], 3, {inf1, inf2, inf3}]; 
    TableForm[Δmatrix]

One of the three first integrals is

integ1  =  
      FirstIntegral[x, u, Δmatrix, 1] ==  c1 // Simplify

The result depends on u'' and now allows us to rewrite all terms containing u'' in the Lie matrices. Next we solve the first integral with respect to u''. Since the integral depends quadratically on u'', we get two solutions

sol1  =  Solve[integ1, (∂_(x, x))u[x]]

The first of the two solutions is used to replace u'' in the Lie matrix

dmat  =  
      Δmatrix /. 
        sol1[[1]]

The simplified Lie matrix again is used to calculate the second first integral of the third-order equation.

integ2  =  FirstIntegral[x, u, dmat, 3] ==  c2 // Simplify

As expected we find the integral depending only on first derivatives of u. Since this integral is linear in u', it is uniquely solvable in u'

sol2  =  Solve[integ2, u'[x]] // Simplify

The resultant expression contains radicals of quadratic polynomials in x. Inserting this result into the Lie matrix, we are able to eliminate the dependencies on u'. We  find

dmat1  =  
      dmat  /.  sol2[[1, 1]] // 
        Simplify

The last step of integration is inserting the Lie matrices into the third first integral depending on u and x

integ3  =  FirstIntegral[x, u, dmat1, 2] ==  c3 // Simplify

The integral contains the dependent variable u again in a linear form. In turn, we end up with a unique solution for the Kamke equation 7.13 which is

solution  =  Solve[integ3, u[x]] // Simplify

The solution depends on three constants c1, c2, and c3, all of which are constants of integration. The parameters α  and β are the parameters of  the original equation. To get a feeling for how the solution evolves, we plot the solution for different parameter sets c1, c2, c3 at fixed α and β.

s1  =  
      Table[(u[x]  /. solution  /. {c1 -> 1, c2 -> (1/2), 
                c3 -> i, α -> 1, β -> (1/2)})
            [[1]], {i, 0, 5}];

Plot[Evaluate[s1], {x, 5, 7}, 
      PlotStyle 
        -> {RGBColor[0,  0,  0.996109], RGBColor[0,  0,  0.62501], 
          RGBColor[0.500008,  0,  0.500008], 
          RGBColor[0.500008,  0,  0.996109], 
          RGBColor[0.996109,  0,  0.500008], 
          RGBColor[0.500008,  0,  0.250004]}, 
      AxesLabel -> {"x", "u"}]

Figure 3

Figure 3 represents the real valued solutions of the equation \(\(\ \(∂\_\(x, x\)u[x] ∂\_{x, 3}u[x] - \)\)\) α 
          (Sqrt[ 1 + (β^2) ((((∂_(x, x))u[x]))^2)])  = 0
        . The different curves represent the solutions for values of  c3  ∈ {0, 1, 2, 3, 4, 5} and fixed values for c2=1/2,c1=1. The parameters of the equation are α=1 and β=1/2.

We note  that the solution is explicitly represented by a complicated expression containing radicals and polynomials. This result is new, since Kamke offers only a parametric representation of the solution. The example shows that we are able to arrive at a solution for higher order ODEs with Lie's procedure. Mathematica is not yet able to handle such types of equations.

Example: 2

The second example for higher order equations is also a third order equation. In this instance, the equation has no physical origin. It is only used for checking  the integration procedure. The equation has some similarity with hydrodynamic equations if the first term of the equation is replaced by a fourth-order derivative. However, the fourth-order equation does not possess the necessary number of symmetries to start the integration process. By replacing the fourth-order derivative with a third-order derivative, the equation obtains the right amount of symmetries to integrate the equation. The equation under consideration is

thirdOrder  =  
      (∂_(x, x, x))u[x] +  
        ℜ𝔢 
          ((∂_x) u[x] (∂_(x, x))u[x] -  
              u[x] (∂_(x, x, x))u[x])

and has a minimum number  of symmetries.

\(infi\  = \ 
      Infinitesimals[thirdOrder\ , u, x, {ℜ𝔢}, 
        SubstitutionRules -> {∂\_{x, 3}u[x]}]\)

The calculation below shows that the functions of MathLie are not able to find all first integrals related to the vector fields (Subscript[(Overscript[v, ⇀]), i]). The symmetries related to the vector fields are given by

inf1  =  
      {{xi[1][x, u]}, {phi[1][x, u]}}  /. infi /. {
            k1 -> 1, k2 -> 0, k3 -> 0}  /. u -> u[x];

inf2  =  
      {{xi[1][x, u]}, {phi[1][x, u]}}  /. infi  /. {
            k1 -> 0, k2 -> 1, k3 -> 0}  /. u -> u[x];

inf3  =  
      {{xi[1][x, u]}, {phi[1][x, u]}}  /. infi  /. {
            k1 -> 0, k2 -> 0, k3 -> 1}  /. u -> u[x];

The right hand side of the third-order equation ω is given by

ω  =  
      (Solve[thirdOrder == 0, (∂_(x, x, x))u[x]])
          [[1, 1, 2]]

The 4×4 Lie matrix follows by inserting the infinitesimals of the vector fields and ω into the function DeltaMatrix[]

Δmatrix = 
      DeltaMatrix[x, u, ω, 3, {inf1, inf2, inf3}]; 
    TableForm[Δmatrix]

The determinant of the Lie matrix is

Det[Δmatrix]

In order to obtain first integrals for the equation, the second line of the above Δmatrix is replaced by the differentials.

integ1  =  
      FirstIntegral[x, u, Δmatrix, 2] ==  c1 // Simplify

The remaining two integrals (Subscript[ψ, 1]) and (Subscript[ψ, 3]) are not accessible to  FirstIntegral[]. The reason for this is that the function Integrate[] of Mathematica cannot solve certain types of integrals. However, the result so far derived is useful to find a solution of the equation. If we look at the first integral (Subscript[ψ, 2]), we observe that this relation is a second order ODE. The solution for u'' clearly shows

sol1  =  Solve[integ1, (∂_(x, x))u[x]] // Simplify

If we transform Rule to Equal, we get

eqh  =  
      (sol1  /.  Rule -> Equal)[[1, 1
        ]]

a second order equation which is now solvable by DSolve[]

sol2  =  DSolve[eqh, u, x]  /.  c1 -> C[3]

Inserting this solution into our original equation, we can verify that the original equation is satisfied

thirdOrder  /.  sol2 // Simplify

meaning that the left-hand side of the equation vanishes, and, in turn equality has been established. Solving  the original equation  with Mathematica, we get

DSolve[thirdOrder == 0, u, x]

showing us that Mathematica is not capable of handling third-order equations in the present form shown above.

The example shows that we can even find solutions for cases of which we know less first integrals than the order of the equation. The procedure of integration splits heterogeneously using different tools to solve the reduced equation. The pseudo-symmetry method allows us to construct a solution at least in an implicit form.

6 Conclusion

We have demonstrated that the well-known term of an integrating factor can be  generalized by using Lie matrices. These matrices are the bases for the multipliers of differential equations. Knowing the multiplier, we are able to derive first integrals of motion by only using  integrations (quadratures). The integrals of motion allow the elimination of the higher derivatives corresponding to a reduction of order in the differential equation. If we are able to construct n first integrals of an nth-order ODE, we are able to solve the equation implicitly. If we fail to determine all  n first integrals, we can find a solution, if at least one first integral is known. Once a minimal number of symmetries of the equation is known, we are always able to arrive at a solution as demonstrated by examples for first-order, second-order, and third-order equations. The method presented, however, is not restricted to these kinds of equations only.

References

S. Lie, Verh. Gesell. d. Wissenschaften zu Christiania, 1874

S. Lie, Theorie der Transformationsgruppen, Vol. I-II, Leipzig, (1888, 1890, 1893), reprinted by Chelsea Publ. Com., New York , (1970)

S. Lie und G. Scheffers, Vorlesungen über Differentialgleichungen mit bekannten infinitesimalen Transformationen, B.G. Teubner, Leipzig, 1891.

W.F. Ames, Nonlinear Ordinary Differential Equations in Transport Processes, Academic Press, New York, 1968.

G. Baumann, Lie Symmetries of Differential Equations, A Mathematica  program to determine Lie symmetries, Wolfram Research Inc. Champaign, Math-Source 0202-622

G. Baumann, MathLie, University of Ulm, 1997.

G.W. Bluman and S. Kumei, Symmetries and Differential Equations, Springer, New York, 1989.

N.H. Ibragimov, CRC Handbook of Lie Group Analysis of Differential Equations, Vol. 1-3, CRC-Press, Boca Raton, 1994, 1995, 1996.

N.H. Ibragimov, Sophus Lie and Harmony in Mathematical Physics, on the 150th Anniversary of His Birth, The Math. Intel. 16, 20-28, (1994).

E. Kamke, Differentialgleichungen, G.B. Teubner, Stuttgart, 1977.

P.J. Olver, Applications of Lie groups to differential equations, Springer, Berlin, 1986.