1.1 Introduction

The magic of Lie symmetries comes from the idea that if non trivial symmetry for an ode can be found, the symmetry makes solving the ode much simpler. For first order differential equation, the symmetry (or rather the tangent vectors associated with the symmetry) are used to convert the ode to canonical coordinates where the original ode becomes a quadrature ode solved by just integration, and this is regardless of how complicated the original ode was.

The ode in canonical coordinates is solved by integration and the solution is converted back to the original natural coordinates giving the solution to the original ode.

For higher order ode’s, the symmetry is used to reduce the order of the ode by one. Making it easier to solve.

Let us start by looking at first order ode’s first. Next section will consider second order ode’s once I learn that subject more.

Given an ode in the space \(\left ( x,y\right ) \,\)  which is called the natural coordinates to distinguish these from the canonical coordinates \(\left ( X,Y\right ) \) used later in which the ode becomes a quadrature. Assume the ode to solve is

\begin{equation} \frac {dy}{dx}=\omega \left ( x,y\right ) \tag {1}\end{equation}

We look for one parameter group of symmetries that when applied to the above ode leaves it in same form. i.e. invariant. This is called Lie continuous group of symmetries which depends on one parameter. It is continuous group and not discrete symmetry group (as say with the case of symmetries for a triangle) because symmetries depend on parameter \(\varepsilon \) which is a real number and whose value is continuous. This Lie parameter is typically called \(\varepsilon \) but some books call it \(\lambda \). These transformations are given by

\begin{align} \bar {x}\left ( x,y\right ) & =f\left ( x,y;\varepsilon \right ) \tag {2}\\ \bar {y}\left ( x,y\right ) & =g\left ( x,y;\varepsilon \right ) \nonumber \end{align}

Where in the above \(f,g\) are independent functions of each others and continuous and analytic in \(\varepsilon \). The \(;\) before \(\varepsilon \) is used instead of comma, to indicate this is a parameter. Eq (2) is the symmetry transformation that when applied to ode (1) results in an ode of same form as (1) but using the new coordinates \(\left ( \bar {x},\bar {y}\right ) \).

We are interested in non trivial transformation (2). This means a transformation that maps points on one solution curve of (1) to a different solution curve. If the transformation maps point \(\left ( x,y\right ) \) to \(\left ( \bar {x},\bar {y}\right ) \,\) on the same solution curve, then it is called a trivial transformation.

The identity transformation in (2) is when \(\varepsilon =0\). Therefore

\begin{align} \bar {x} & =f\left ( x,y;0\right ) =x\tag {3}\\ \bar {y} & =g\left ( x,y;0\right ) =y\nonumber \end{align}

Every Lie group must have an identity element and this is given by \(\varepsilon =0\). (The Lie group must also have a unique inverse element for each element in the set of transformations, and the set be closed under application of the transformation. But not all Lie groups are commutative or associative (i.e. Abelian). Consult more mathematical book on properties of Lie group if interested. These notes will concentrate on the algorithmic aspect.

Looking back at (2), since \(f,g\) are analytic in \(\varepsilon \), they can be expanded in Taylor series near \(\varepsilon =0\) which results in

\begin{align} \bar {x} & =f\left ( x,y;0\right ) +\left . \frac {\partial f}{\partial \varepsilon }\right \vert _{\varepsilon =0}\varepsilon +O\left ( \varepsilon ^{2}\right ) \tag {4}\\ \bar {y} & =g\left ( x,y;0\right ) +\left . \frac {\partial g}{\partial \varepsilon }\right \vert _{\varepsilon =0}\varepsilon +O\left ( \varepsilon ^{2}\right ) \nonumber \end{align}

Choosing \(\varepsilon \) to be very small, higher order terms \(O\left ( \varepsilon ^{2}\right ) \) can be ignored. Also since \(f\left ( x,y;0\right ) =x\) and \(g\left ( x,y;0\right ) =y\) because \(\varepsilon =0\) is the identity transformation, (4) simplifies to

\begin{align} \bar {x} & =x+\left . \frac {\partial f}{\partial \varepsilon }\right \vert _{\varepsilon =0}\varepsilon \tag {5}\\ \bar {y} & =y+\left . \frac {\partial g}{\partial \varepsilon }\right \vert _{\varepsilon =0}\varepsilon \nonumber \end{align}

Typically, the term \(\left . \frac {\partial f}{\partial \varepsilon }\right \vert _{\varepsilon =0}\) is called \(\xi \left ( x,y\right ) \) and \(\left . \frac {\partial g}{\partial \varepsilon }\right \vert _{\varepsilon =0}\) is called \(\eta \left ( x,y\right ) \,\). Therefore (5) is written as

\begin{align} \bar {x} & =x+\varepsilon \xi \left ( x,y\right ) \tag {6}\\ \bar {y} & =y+\varepsilon \eta \left ( x,y\right ) \nonumber \end{align}

The functions \(\xi \left ( x,y\right ) ,\eta \left ( x,y\right ) \) are the most important quantities in Lie group symmetry. These are called the Lie infinitesimals. Geometrically, these are the tangent vectors at the identity on the tangent plane which is the Lie algebra. If these tangent vectors can be found for the given ode, then the ode can be now transformed to quadrature in the canonical coordinates, regardless of how complicated the original ode was.

The transformation (6) must leave the ode invariant. This is called the Lie invariance condition. Therefore

\begin{equation} \frac {d\bar {y}}{d\bar {x}}=\omega \left ( \bar {x},\bar {y}\right ) \tag {7}\end{equation}

Whenever

\begin{equation} \frac {dy}{dx}=\omega \left ( x,y\right ) \tag {1}\end{equation}

Hence, starting with (7) it can be written as

\begin{equation} \frac {d\bar {y}}{d\bar {x}}=\frac {\frac {d\bar {y}}{dx}}{\frac {d\bar {x}}{dx}}=\omega \left ( \bar {x},\bar {y}\right ) \tag {8}\end{equation}

Where \(\frac {d\bar {y}}{dx}\) is the total derivative with respect to \(x\) and \(\frac {d\bar {x}}{dx}\)is the total derivative with respect to \(x\). But

\begin{align*} \frac {d\bar {y}}{dx} & =\frac {d}{dx}\left ( y+\varepsilon \eta \left ( x,y\right ) \right ) \\ & =\frac {dy}{dx}+\varepsilon \frac {d\eta \left ( x,y\right ) }{dx}\\ & =\frac {dy}{dx}+\varepsilon \left ( \eta _{x}+\eta _{y}\frac {dy}{dx}\right ) \end{align*}

But \(\frac {dy}{dx}=\omega \left ( x,y\right ) \), therefore the above becomes

\begin{equation} \frac {d\bar {y}}{dx}=\omega +\varepsilon \left ( \eta _{x}+\omega \eta _{y}\right ) \tag {9}\end{equation}

Similarly

\begin{align} \frac {d\bar {x}}{dx} & =\frac {d}{dx}\left ( x+\varepsilon \xi \left ( x,y\right ) \right ) \nonumber \\ & =1+\frac {d}{dx}\left ( \varepsilon \xi \left ( x,y\right ) \right ) \nonumber \\ & =1+\varepsilon \left ( \xi _{x}+\xi _{y}\frac {dy}{dx}\right ) \nonumber \\ & =1+\varepsilon \left ( \xi _{x}+\omega \xi _{y}\right ) \tag {10}\end{align}

Substituting (9,10) into (8) gives

\begin{align*} \frac {d\bar {y}}{d\bar {x}} & =\frac {\frac {d\bar {y}}{dx}}{\frac {d\bar {x}}{dx}}\\ & =\omega \left ( \bar {x},\bar {y}\right ) \end{align*}

Hence the above becomes

\begin{align} \frac {\omega +\varepsilon \left ( \eta _{x}+\omega \eta _{y}\right ) }{1+\varepsilon \left ( \xi _{x}+\omega \xi _{y}\right ) } & =\omega \left ( \bar {x},\bar {y}\right ) \tag {11}\\ & =\omega \left ( x+\varepsilon \xi ,y+\varepsilon \eta \right ) \nonumber \end{align}

Looking at left side of (11), the term \(\frac {1}{1+\varepsilon \left ( \xi _{x}+\omega \xi _{y}\right ) }\) can be moved to numerator using binomial expansion. Since \(\frac {1}{1+x}=\left ( 1+x\right ) ^{-1}\), then using

\[ \left ( 1+x\right ) ^{n}=1+nx+\frac {n\left ( n-1\right ) }{2!}x^{2}+\cdots \]

Therefore, using \(n=-1\) the term \(\left ( 1+\varepsilon \left ( \xi _{x}+\omega \xi _{y}\right ) \right ) ^{-1}\) becomes (using \(x\equiv \varepsilon \left ( \xi _{x}+\omega \xi _{y}\right ) \) here)

\[ \left ( 1+\varepsilon \left ( \xi _{x}+\omega \xi _{y}\right ) \right ) ^{-1}=1-\varepsilon \left ( \xi _{x}+\omega \xi _{y}\right ) +\left ( \varepsilon \left ( \xi _{x}+\omega \xi _{y}\right ) \right ) ^{2}-\cdots \]

But since \(\varepsilon \) is very small, all terms with \(O\left ( \varepsilon ^{2}\right ) \) and higher can be ignored. The above becomes

\[ \left ( 1+\varepsilon \left ( \xi _{x}+\omega \xi _{y}\right ) \right ) ^{-1}=1-\varepsilon \left ( \xi _{x}+\omega \xi _{y}\right ) \]

Using the above into (11) gives

\begin{align*} \left ( \omega +\varepsilon \left ( \eta _{x}+\omega \eta _{y}\right ) \right ) \left ( 1-\varepsilon \left ( \xi _{x}+\omega \xi _{y}\right ) \right ) & =\omega \left ( x+\varepsilon \xi ,y+\varepsilon \eta \right ) \\ \omega -\varepsilon \omega \left ( \xi _{x}+\omega \xi _{y}\right ) +\varepsilon \left ( \eta _{x}+\omega \eta _{y}\right ) -\varepsilon ^{2}\left ( \left ( \eta _{x}+\omega \eta _{y}\right ) \left ( \xi _{x}+\omega \xi _{y}\right ) \right ) & =\omega \left ( x+\varepsilon \xi ,y+\varepsilon \eta \right ) \end{align*}

But since \(\varepsilon \) is very small, the above simplifies to

\begin{equation} \omega -\varepsilon \omega \left ( \xi _{x}+\omega \xi _{y}\right ) +\varepsilon \left ( \eta _{x}+\omega \eta _{y}\right ) =\omega \left ( x+\varepsilon \xi ,y+\varepsilon \eta \right ) \tag {12}\end{equation}

Now we expand the right side of the above using Taylor series. Since \(\omega \left ( x,y\right ) \) is function of two variables, we use Taylor series expansion for two variables. The expansion is made around the point \(\left ( x,y\right ) \) as shown in this diagram

Figure 1: Expanding \(\omega (x,y)\) in Taylor series

Hence approximation of the right side of (12) is

\begin{align*} \omega \left ( x+\varepsilon \xi ,y+\varepsilon \eta \right ) & =\omega \left ( x,y\right ) +\varepsilon \xi \omega _{x}+\varepsilon \eta \omega _{y}+O\left ( \varepsilon ^{2}\right ) +\cdots \\ & =\omega \left ( x,y\right ) +\varepsilon \left ( \omega _{x}\xi +\omega _{y}\eta \right ) \end{align*}

Substituting the above in (12) results in

\begin{align*} \omega -\varepsilon \omega \left ( \xi _{x}+\omega \xi _{y}\right ) +\varepsilon \left ( \eta _{x}+\omega \eta _{y}\right ) & =\omega +\varepsilon \left ( \omega _{x}\xi +\omega _{y}\eta \right ) \\ -\varepsilon \omega \left ( \xi _{x}+\omega \xi _{y}\right ) +\varepsilon \left ( \eta _{x}+\omega \eta _{y}\right ) & =\varepsilon \left ( \omega _{x}\xi +\omega _{y}\eta \right ) \\ -\omega \left ( \xi _{x}+\omega \xi _{y}\right ) +\left ( \eta _{x}+\omega \eta _{y}\right ) & =\omega _{x}\xi +\omega _{y}\eta \\ -\omega \xi _{x}-\omega ^{2}\xi _{y}+\eta _{x}+\omega \eta _{y}-\omega _{x}\xi -\omega _{y}\eta & =0 \end{align*}

Or

\begin{equation} \fbox {$\eta _x+\omega \left ( \eta _y-\xi _x\right ) -\omega ^2\xi _y-\omega _x\xi -\omega _y\eta =0$} \tag {13}\end{equation}

The above equation (13) is what is used to solve for \(\xi ,\eta \). It is the linearized symmetry condition. There is an additional constraint not mentioned above which is that we must have

\[ \bar {x}_{x}\bar {y}_{y}\neq \bar {x}_{y}\bar {y}_{x}\]

PDE (13) is solved using ansatz. Examples below show how this is done.

OK, now let us assume that we have found the tangent vectors \(\xi ,\eta \) by solving (13). What to do next? Here comes the main point of using Lie symmetry. To make any first ode \(\frac {dy}{dx}=\omega \left ( x,y\right ) \) solvable by quadrature, we impose transformation of the form

\begin{align} X & =x\tag {14}\\ Y & =y+\varepsilon \nonumber \end{align}

This transformation maps points on one solution curve to points on another solution curve, but by only changing the \(y\) coordinate. Hence it is vertical displacement or shift.  The \(X,Y\) are called the canonical coordinates and these are functions of \(x,y\). In other words, the above can be written as

\begin{align} X\left ( x,y\right ) & =x\tag {14(a)}\\ Y\left ( x,y\right ) & =y+\varepsilon \tag {14(b)}\end{align}

From (14 a) we obtain

\begin{align} \left . \frac {\partial X}{\partial \varepsilon }\right \vert _{\varepsilon =0} & =0\nonumber \\ & =X_{x}\left . \frac {dx}{d\varepsilon }\right \vert _{\varepsilon =0}+X_{y}\left . \frac {dy}{d\varepsilon }\right \vert _{\varepsilon =0} \tag {15}\end{align}

But what are \(\left . \frac {dx}{d\varepsilon }\right \vert _{\varepsilon =0}\) and \(\left . \frac {dy}{d\varepsilon }\right \vert _{\varepsilon =0}\)? These come from the earlier relation we had, which is

\[ \bar {x}=x+\varepsilon \xi \]

Therefore

\[ \frac {d\bar {x}}{d\varepsilon }=\xi \]

But at \(\varepsilon =0\,\), we have \(\bar {x}=x\) since the identity. Therefore

\[ \left . \frac {dx}{d\varepsilon }\right \vert _{\varepsilon =0}=\xi \]

Similarly

\[ \left . \frac {dy}{d\varepsilon }\right \vert _{\varepsilon =0}=\eta \]

Now (14 a) becomes

\begin{equation} 0=X_{x}\xi +X_{y}\eta \tag {16}\end{equation}

Similarly (14 b) gives

\begin{align*} \left . \frac {\partial Y}{\partial \varepsilon }\right \vert _{\varepsilon =0} & =1\\ & =Y_{x}\left . \frac {dx}{d\varepsilon }\right \vert _{\varepsilon =0}+Y_{y}\left . \frac {dy}{d\varepsilon }\right \vert _{\varepsilon =0}\\ & =Y_{x}\xi +Y_{y}\eta \end{align*}

Hence

\begin{equation} 1=Y_{x}\xi +Y_{y}\eta \tag {17}\end{equation}

Equations (16,17) are now solved for \(Y\left ( x,y\right ) \) and \(X\left ( x,y\right ) \). Once these are found, we will find the ode in canonical coordinates to be of the form

\[ \frac {dY}{dX}=G\left ( X\right ) \]

always. In other words, one that can be solved by just integration. Equations (16,17) are again the following

\begin{align} 0 & =X_{x}\xi +X_{y}\eta \tag {16}\\ 1 & =Y_{x}\xi +Y_{y}\eta \tag {17}\end{align}

Each is solved by the standard method of characteristics in PDE. To apply this method, let us start with (16). We assume the characteristics variable is \(\tau \) and say that \(X\left ( \tau \right ) \equiv X\left ( x\left ( \tau \right ) ,y\left ( \tau \right ) \right ) \). Hence

\begin{equation} \frac {dX}{d\tau }=X_{x}\frac {dx}{d\tau }+X_{y}\frac {dy}{d\tau } \tag {16 a}\end{equation}

Comparing (16, 16a) shows that

\begin{align} \frac {dX}{d\tau } & =0\tag {18}\\ \frac {dx}{d\tau } & =\xi \nonumber \\ \frac {dy}{d\tau } & =\eta \nonumber \end{align}

And similarly, we let \(Y\left ( \tau \right ) \equiv Y\left ( x\left ( \tau \right ) ,y\left ( \tau \right ) \right ) \), hence

\begin{equation} \frac {dY}{d\tau }=Y_{x}\frac {dx}{d\tau }+Y_{y}\frac {dy}{d\tau } \tag {17 a}\end{equation}

Comparing (17, 17a) shows that

\begin{align} \frac {dY}{d\tau } & =1\tag {19}\\ \frac {dx}{d\tau } & =\xi \nonumber \\ \frac {dy}{d\tau } & =\eta \nonumber \end{align}

Now we solve (18,19) for \(Y\left ( x,y\right ) ,X\left ( x,y\right ) \) since \(\xi ,\eta \) are known. Examples below show how this is done. At the end of this process,  we will have \(Y,X\) and now we evaluate \(\frac {dY}{dX}\) using

\begin{align*} \frac {dY}{dX} & =\frac {\frac {dY}{dx}}{\frac {dX}{dx}}\\ & =\frac {Y_{x}+Y_{y}\frac {dy}{dx}}{X_{x}+X_{y}\frac {dy}{dx}}\\ & =\frac {Y_{x}+Y_{y}\omega }{X_{x}+X_{y}\omega }\end{align*}

Everything on the RHS above is known, since we have solved for \(Y,X\) from (18,19). This complete the process. All what is left is to solve this resulting ode, which should always be

\[ \frac {dY}{dX}=G\left ( X\right ) \]

This will give solution \(Y\left ( X\right ) \). Now this is converted back to \(y\left ( x\right ) \) using the known \(X\left ( x,y\right ) ,Y\left ( x,y\right ) \) relations. This complete the algorithm using Lie symmetry for first order ode. Here is a summary of the steps.

1. First order ode \(\frac {dy}{dx}=\omega \left ( x,y\right ) \) is given.
2. Solve the Lie linearized symmetry pde \(\eta _{x}+\omega \left ( \eta _{y}-\xi _{x}\right ) -\omega ^{2}\xi _{y}-\omega _{x}\xi -\omega _{y}\eta =0\) for the tangent vectors \(\eta ,\xi \).
3. Now that tangent vectors \(\eta ,\xi \) are found, these are used to find canonical coordinates \(X\left ( x,y\right ) ,Y\left ( x,y\right ) \). without knowing tangent vectors the canonical coordinates can not be found.
4. Impose transformation to canonical coordinates using \(X=x,Y=y+\varepsilon \). Solve the resulting pair of pde which result due to this. These are \(0=X_{x}\xi +X_{y}\eta \) and \(1=Y_{x}\xi +Y_{y}\eta \).
5. Using method of method of characteristics, the above pair of pde’s are solved for \(X\left ( x,y\right ) ,Y\left ( x,y\right ) .\)
6. Set up the ode \(\frac {dY}{dX}=\frac {Y_{x}+Y_{y}\omega }{X_{x}+X_{y}\omega }\) and simplify the RHS. It should come out as function of \(X\) only.
7. Solve the ode \(\frac {dY}{dX}=G\left ( X\right ) \) by integration to find \(Y\left ( X\right ) \).
8. Replace all \(Y,X\) in the solution above by \(y,x\) to convert the solution back from canonical coordinates \(\left ( X,Y\right ) \) to natural coordinates \(\left ( x,y\right ) \).