Lie point symmetries for first order ODEs

A symbolic algorithm for the maxima CAS Nijso Beishuizen

Introduction

The symmetry method is a powerful method to solve differential equations by finding invariant transformations (point symmetries) of the ODE. For first order ODEs, an integrating factor can be constructed from the symmetry and the ODE can then be solved. For first order ODEs, there is no closed form algorithm to find a symmetry in general, but there are algorithms for finding symmetries of a certain form. In Cheb-Terrab and Kolokolnikov [1], a very general algorithm is presented that searches for linear symmetries of the form

$[\xi = F(x),\eta=P(x)y+Q(x)]$.

In Cheb-Terrab and Roche [2], an algorithm is presented for finding several simpler symmetries of the form

(1) $[\xi=F(x)*G(y), \eta=0]$, $[\xi=0, \eta=F(x)*G(y)]$

(2) $[\xi=F(x)+G(y), \eta=0]$, $[\xi=0, \eta=F(x)+G(y)]$

(3) $[\xi=F(x), \eta=G(x)]$, $[\xi=F(y), \eta=G(y)]$

(4) $[\xi=F(x), \eta=G(y)]$, $[\xi=F(y), \eta=G(x)]$

(5) $[\xi=ax+by+c, \eta=fx+gy+h]$

We see that the linear symmetries presented in [1] are a generalization of the search for symmetries of the form (3) presented in [2].

These methods have been implemented in the function ode1solve and can be found in the file ode1_lie.mac. Additionally, several other methods have been implemented that search for symmetries of specific ODEs. These algorithms are very fast but also limited in applicability. The implemented methods are symmetry searches for ODEs of the type: linear, inverse-linear, separable, exact, Riccati, Bernoulli and Abel. These symmetries can also be found using the search for more general symmetries mentioned above, but these dedicated methods are usually much faster.

Performance of the implementation

As a performance test, the database of Kamke's first order ODEs was used [3]. In the database are 367 ODEs of first order and first degree. A number of these ODEs (Cheb-Terrab and Kolokolnikov mention 20 of the first order - first degree) are unsolvable because they are too general (these are ODEs 47, 48, 50, 55, 56, 74, 79, 82, 202, 219, 234, 235, 237, 250, 253, 269, 331). 178 of the ODEs are of type separable, linear, inverse linear and Bernoulli and therefore considered straightforward to solve. They were not considered in the evaluation of the symmetry method in [2]. These equations (as well as all homogeneous and Ricatti and Abel odes with constant invariants) all have linear symmetries. Of the 246 ODEs of first order and first degree that are left, 133 were solved by the Cheb-Terrab and Roche method. With the method of Cheb-Terrab and Kolokolnikov, another 20 ODEs of first order and first degree with linear symmetries can be found.

Our implementation currently finds (nontrivial) symmetries for 298 of the 347 solvable ODEs of first order and first degree, which amounts to a success rate of 85%. Every nontrivial symmetry of a first order ODE directly leads to an integrating factor, and therefore to a solution of the ODE. Most ODEs for which no symmetries could be found are currently outside the scope of the method. They are Riccati or (nonconstant invariant) Abel ODEs with special function solutions like Bessel functions.

Overview of functions, options and keywords

Function: ode1solve(ode,y,x,[options])

Returns: a list with the analytical solutions of the ODE, or false when no solutions are found.

ode1solve is an Ordinary Differential Equation (ODE) solver for first order ODEs using symmetries. If point symmetries can be found for the input ODE(s), ode1solve will use these symmetries to solve construct and integrating factor and solve the ODE. The input is a first order ode of the form F('diff(y,x),y,x)=0. ode1solve will first try to solve the ode explicitly for 'diff(y,x) and then solve each of the branches. The default is to return the explicit analytical solution(s) as a list. The list may contain more than one solution.



In [ ]:

    
kill(all);batch("~/mathematics/maxima_files/kamke1_1.mac");batch("~/mathematics/maxima_files/ode1_lie.mac");DEBUGFLAG:1;



In [3]:

    
ode1solve('diff(y,x)+x^2+5,y,x);









    Out[3]:





\[\tag{${\it \%o}_{157}$}\left[ y=-\frac{x^3+15\,x-{\it \%c}}{3} \right] \]

Optional arguments are:

(bool) returnSymmetries={true,false}

if returnSymmetries=true, ode1solve will return the symmetries of the ode as a list [$\xi$,$\eta$

(bool) returnIntegratingFactor={true,false}

if returnIntegratingFactor=true, ode1Solve will first try to find the symmetries and then use them to construct and return an integrating factor $\mu$ for the ode as a list [$\mu$]. A symmetry of the ODE has to be found first.

(bool) returnSolution={true,false}

if returnSolution=true, ode1solve will use the integrating factor to find the general analytic solution(s) and will return it as a list [sol].

(bool) returnExplicit={true,false}

If returnExplicit is true, ode1solve will try to find the explicit solution of the ode. A symmetry and an integrating factor have to be found first. When returnExplicit is false, the implicit solution will be returned as a first integral.

(bool) tryInverse={true,false}

When tryInverse=true, after all regular symmetry methods have been tried and when no symmetry could be found, all methods are tried again but with the inverse ode $y'=1/(a(y)x+b(y)$

(bool) useCanonical={true,false}

when useCanonical is true, the ode will first be cast into some canonical form before ode1solve tries to solve it. When the original ODE is exact, the canonical form might not be exact.

(string) useMethod="method"

Use only "method" to find symmetries, integrating factors or solutions. Method can be a string or a list of strings. Options are: "all", "quadrature", "fx","fy","exact","separable", "riccati", "Riccati", "abel", "Abel", "linear", "inverse-linear","bernoulli","Bernoulli","symmetry1","symmetry2","symmetry3","symmetry4","symmetry5","muc","symmetry5b","symmetry6". method can be a list, and the methods in the list will then be tried in the order found in the list. the method "all" is a predefined list: ["separable","exact","linear","inverse-linear","bernoulli","riccati", "symmetry1","symmetry2","symmetry3","symmetry4","symmetry5","symmetry5b","symmetry6"] methods available are :

fx
fy: solves an ode of the form $y'=f(y)$
quadrature: solves an ode that can be solved using quadrature. The ode is of the form $y'=f(x)$ or $y'=f(y)$.
separable: solves separable odes. The ode is of the form $y'=f(x)*g(y)$
linear: solves linear odes. The ode is of the form $y'=f(x)y + g(x)$
inverse-linear: solves inverse-linear odes. The ode is of the form $y'=1/(f(y)x + g(y))$
exact: solves exact odes. The ode is of the form $y'=\frac{Q(x,y)}{P(x,y)}$ and $\frac{dP}{dx}-\frac{dQ}{dy}=0$
abel, Abel: solves Abel odes. The ode is of the form $y'=a(x)y^3+b(x)y^2+c(x)y+d(x)$ or can be written in this form using a transformation. Only Abel odes with constant invariants are solved. These are also solved by symmetry method 6.
bernoulli,Bernoulli: solves Bernoulli odes. The ode is of the form $y'=f(x)*y^n + g(x)*y$.
riccati,Riccati: solves Riccati odes. The ode is of the form $y'=f(x)*y^2 + g(x)*y + h(x)$. Only Riccati odes with constant invariant are solved with this routine. Many other Riccati ode's are solved with the symmetry methods.
symmetry1..4: The symmetry methods of Cheb-Terrab and Roche.
symmetry5: Tries to find symmetries of polynomial form. The default is a third degree multivariate polynomial in x and y (In Cheb-Terrab and Roche, a first degree polynomial was used).
symmetry6: The symmetry method of Cheb-Terrab and Kolokolnikov. Note that this method makes a call to method symmetry4.
all: use all available methods in the following order: quadrature, separable, exact, linear, inverse-linear, bernoulli, riccati, symmetry1, symmetry2, symmetry3, symmetry4, symmetry5, symmetry6.



In [4]:

    
DEBUGFLAG:1;
ode:'diff(y,x)+p(x)*y=q(x)*y^4;
ode1solve(ode,y,x,'returnExplicit=false,'useMethod="bernoulli");









    





Bernoulli ODE: y' +a(x)*y = b(x)*y^n, with [b(x),n,a(x)] =  \(q\left(x\right)\) ,  \(4\) ,  \(-p\left(x\right)\)






    Out[4]:





\[\tag{${\it \%o}_{158}$}1\]






    Out[4]:





\[\tag{${\it \%o}_{159}$}\frac{d}{d\,x}\,y+p\left(x\right)\,y=q\left(x\right)\,y^4\]






    Out[4]:





\[\tag{${\it \%o}_{160}$}\left[ \frac{e^ {- 3\,\int {p\left(x\right)}{\;dx} }\,\left(3\,e^{3\,\int {p\left(x\right)}{\;dx}}\,\left(\int {e^ {- 3\,\int {p\left(x\right)}{\;dx} }\,q\left(x\right)}{\;dx}\right)\,y^3+1\right)}{y^3}={\it \%c} \right] \]

The default of ode1solve is to use all these methods sequentially, in this order, to find a solution of the ode $y'=\Phi(x,y)$. The method abel is not called because all Abel ODEs with constant invariants have linear symmetries and can be solved by the method "symmetry6". Abel ODEs with non-constant invariants cannot be solved in general, althought there are a couple of specific cases where ode1solve is able to find a symmetry. All the available methods can be called individually. When for instance 'useMethod="symmetry1" is chosen, ode1solve searches for symmetries of type $[\xi=F(x)*G(y), \eta=0]$. If the ODE instead has the inverse symmetries $[\xi=0, \eta=F(x)*G(y)]$, then these inverse symmetries can be found by looking for symmetries of the inverse ode $y' =\frac{1}{\Phi(y,x)}$. The inverse of an ODE can be found by using the change of variables $x\rightarrow y^{*}(x^{*})$ and $y(x)\rightarrow x^{*}$ and then renaming $y^{*}$ as $y$ and $x^{*}$ as $x$, see e.g. [2].

keywords for the method of undetermined coefficients

The method of undetermined coefficients is called "symmetry5", an alias is "muc". This method tries to find a symmetry of an ODE using the method of undetermined coefficients. When no keywords are given by the user, we search for polynomial symmetries of the form

$ax+by+cxy+dx^2+...$,

where the default polynomial order is 2. This is determined by the global keyword SYM5DEGREE and can be changed by the user. It turns out that a polynomial degree of 5 or more is computationally too expensive for even the most patient user.

Symmetry5b is another implementation of the method of undetermined coefficients where we search for rational symmetries of the form $\frac{P}{Q}$, with $P$ and $Q$ polynomials in the undetermined variables (usually the dependent variable $y$ and independent variable $x$ in $\frac{dy}{dx}=\Psi(x,y)$). Symmetry5b looks for first for symmetries of the form $[0,F/G]$ where $F$ is a polynomial $F(x,y)$ and $G=P$ is the polynomial expression in the denominator of the ODE $\frac{dy}{dx}=\frac{Q}{P}$. The degree of the polynomial $F$ is fixed at 5.
When these symmetries can not be found, symmetry5b looks for symmetries of the form $[0,F/G]$ and $[F/G,0]$, where $F$ and $G$ are both undetermined polynomials. The degree of the polynomials are set to SYM5DEGREE-1. if SYM5RAT is true, symmetry5b additionally searches for polynomials of the form

$[(ax+by+c)/(dx+ey+f),(gx+hy+i)/(jx+ky+l)]$.

The default for SYM5RAT is false because this method is very expensive and does not uniquely solve any ODE in the kamke database. The user can also define her own predefined shapes for the symmetries. This method requires two keywords: xi and eta. xi and eta define the shape of the symmetry, e.g. $xi=a*x+b*y+c$ or $xi=a*sin(x)$
A separate tutorial on using the method of undetermined coefficients is symmetry5.ipynb
</dl>

Global variables

System variable: (bool) FIX_INTEGRATION_CONSTANT default: true When FIX_INTEGRATION_CONSTANT=true, the integration constant will be fixed at maxima's system variable integration_constant (which is usually %c). When FIX_INTEGRATION_CONSTANT=false, then every call to ode1SolveWithIntegrationFactor (this is called from ode1solve as well) will result in a new integration constant %c1, %c2, etc. where the initial value for the counter is taken from maxima's system variable integration_constant_counter.

System variable: (integer) DEBUGFLAG={0,1,2,3,4,5} default: 1 When DEBUGFLAG=0, no information except the solution or errors are printed. More intermediate information is given for increasing values of DEBUGFLAG

0: no intermediate information
1: warning messages
2: main results from intermediate calculations (symmetries, integrating factor)
3: additional intermediate results
4: More intermediate results, for debugging
5: A lot of output, for serious debugging

System variable: (string) method default: false

method returns the last successful method used to solve an ODE. if the last call to the solver failed, method is false.

System variable: (string) reason

default: "none"

If known, reason returns information on why the solver was not able to find a solution, for instance if the ODE was found to be an Abel ODE with non-constant invariants.

Option variable: (bool) returnSymmetries

default: false

When true, returnSymmetries returns the symmetries found as a list.

Option variable: (bool) returnIntegratingFactor

default: false

When true, returnIntegratingFactor returns the integrating factors found as a list.

Option variable: (bool) returnExplicit

default: true

When true, returnExplicit tries to return an explicit solution in the dependent variable.

Option variable: (bool) tryInverse

default: false

When true, for every symmetry method that ode1solve tries (trying to find a symmetry $X=[\xi,\eta]$), after it has failed it will then try to find the inverse symmetry $X=[\eta,\xi]$.

Option variable: (bool) SYM5RAT

default: false

When SYM5RAT=true, ode1solve will try to find rational-polynomial symmetries of the ODE. This is usually time consuming.

Option variable: (bool) CHECKSYM

default: true

When true, checks if the symmetry and integrating factor is valid.

function: ode1sym(ode,y,x) an alias for ode1solve that will only return the symmetries



In [5]:

    
ode1sym('diff(y,x)=f(x)*g(y),y,x);









    Out[5]:





\[\tag{${\it \%o}_{161}$}\left[ \frac{1}{f\left(x\right)} , 0 \right] \]

function: ode1intfac(ode,y,x) an alias for ode1solve that will only return the integrating factors



In [6]:

    
ode1intfac('diff(y,x)=f(x)*g(y),y,x);









    Out[6]:





\[\tag{${\it \%o}_{162}$}\left[ \frac{1}{g\left(y\right)} \right] \]

function: (list) isBernoulli(ode,y,x) returns a list $[b(x),n,a(x)]$ of the ODE $y'+a(x)y=b(x)y^n$ if the input ode is a Bernoulli ODE, or false otherwise.



In [7]:

    
isBernoulli('diff(y,x)-x*y^2-3*x*y,y,x);









    





Bernoulli ODE: y' +a(x)*y = b(x)*y^n, with [b(x),n,a(x)] =  \(x\) ,  \(2\) ,  \(3\,x\)






    Out[7]:





\[\tag{${\it \%o}_{163}$}\left[ x , 2 , 3\,x \right] \]

function: (list) isRiccati(ode,y,x) returns the coefficients of the Riccati ODE $y' = A(x)+B(x)y+C(x)y^2$.



In [8]:

    
isRiccati('diff(y,x)=x+y^2,y,x);









    





Riccati ODE: y' = a(x) + b(x)*y + c(x)*y^2, with [a(x),b(x),c(x)] =  \(x\) ,  \(0\) ,  \(1\)






    Out[8]:





\[\tag{${\it \%o}_{164}$}\left[ x , 0 , 1 \right] \]

function: (list) isChini(ode,y,x) returns the coefficients [A,n,B,C] if the input ODE is a Chini ODE $y' = C(x) +B(x)y +A(x)y^n$, and false otherwise



In [9]:

    
isChini('diff(y,x)=x*y^4+x^2*y+5,y,x);









    Out[9]:





\[\tag{${\it \%o}_{165}$}\left[ x , 4 , x^2 , 5 \right] \]

function: (list) isAbel(ode,y,x) returns a list of the coefficients of the abel ODE, either in first kind form or second kind form if the input ODE is an Abel ODE, or false otherwise. Abel's ODE of the first kind is $y' = a_0(x)+a_1(x)y+a_2(x)y^2 +a_3(x)y^3$. Abel's ODE of the second kind is $y' = \frac{a_0(x)+a_1(x)y+a_2(x)y^2 +a_3(x)y^3}{a(x)y+b(x)}$.



In [10]:

    
isAbel('diff(y,x)=5+x*y+x^2*y^2+x^3*y^3,y,x);









    Out[10]:





\[\tag{${\it \%o}_{166}$}\left[ {\it \_a}_{0}=5 , {\it \_a}_{1}=x , {\it \_a}_{2}=x^2 , {\it \_a}_{3}=x^3 , {\it \_y}=y \right] \]



In [11]:

    
isAbel('diff(y,x)=(5+x*y+x^2*y^2+x^3*y^3)/(x+y),y,x);









    Out[11]:





\[\tag{${\it \%o}_{167}$}\left[ \left[ {\it \_a}_{0}=5 , {\it \_a}_{1}=x , {\it \_a}_{2}=x^2 , {\it \_a}_{3}=x^3 , {\it \_y}=y \right]  , \left[ {\it \_a}=1 , {\it \_b}=x , {\it \_x}=y \right]  \right] \]

function: (bool) isAbelFirst(ode,y,x) returns true if the input ode is an Abel ODE of the first kind $y' = a_0(x)+a_1(x)y+a_2(x)y^2 +a_3(x)y^3$.

function: (bool) isAbelSecond(ode,y,x) returns true if the input ode is an Abel ODE of the second kind $y' = \frac{a_0(x)+a_1(x)y+a_2(x)y^2 +a_3(x)y^3}{a(x)y+b(x)}$.

function: ode1AbelSecond2First(ode,y,x) converts an Abel ODE of the second kind to an Abel ODE of the first kind

function: ode1AbelFirst2Second(ode,y,x) converts an Abel ODE of the first kind to an Abel ODE of the second kind

function: ode1AbelRNF(ode,y,x) converts an Abel ODE to its rational normal form, either $w'+aw3+bw=0$ or $w'+aw^3+bw+1=0$

function: ode1AbelSymmetries(ode,y,x) Computes the symmetries of an Abel ODE.

function: (bool) isIntegratingFactor(mu,ode,y,x) returns true if the input expression mu is an integrating factor of the ode.

function: ode1_CanonicalForm(ode,y,x) Write the input ODE into a 'canonical' form, meaning that we try to write the input ODE as an explicit ode $y'=\omega(x,y)$ and we try to simplify the right hand side.

function: ode1PfaffianForm(ode,y,x) writes the input ODE $\frac{dy}{dx}=\frac{P(y,x)}{Q(y,x)}$ in Pfaffian form $P(y,x)dx-Q(y,x)dy=0$ and returns $[P,Q]$.



In [12]:

    
ode1PfaffianForm('diff(y,x)=f(x)/g(y),y,x);









    Out[12]:





\[\tag{${\it \%o}_{168}$}\left[ g\left(y\right) , f\left(x\right) \right] \]

function: checkSymmetries(X,ode,y,x) substitutes the input symmetry into the linearized symmetry condition of the ode and returns the result. If the result is 0, then X is a symmetry of the ODE.



In [13]:

    
X:ode1sym('diff(y,x)=f(x)*g(y),y,x);
checkSymmetries(X,'diff(y,x)=f(x)*g(y),y,x);









    Out[13]:





\[\tag{${\it \%o}_{169}$}\left[ \frac{1}{f\left(x\right)} , 0 \right] \]






    Out[13]:





\[\tag{${\it \%o}_{170}$}0\]

function: reverseODE1(ode,y,x) returns the reverse ODE of the input ode, that is we interchange variables: $x^{*}=y$ and $y^{*}=x$ and return an ode $\frac{dy^{*}}{dx^{*}}=\omega(x^{*},y^{*})$.



In [14]:

    
reverseODE1('diff(y,x)=x*y,y,x);









    Out[14]:





\[\tag{${\it \%o}_{171}$}\left[ \frac{d}{d\,x}\,y=\frac{1}{x\,y} \right] \]

function: rootsexpand(expr) splits an expression containing square roots and factorizes the square roots such that $sqrt((p_1*p_2*...p_n)/(q_1*q_2*..q_m))= sqrt(p_1)*sqrt(p_2)*...sqrt(p_n)/(sqrt(q_1)*sqrt(q_2)*..sqrt(q_m))$. This is the inverse of rootscontract.



In [15]:

    
rootsexpand(sqrt(x^2*y*a^3));









    Out[15]:





\[\tag{${\it \%o}_{172}$}a^{\frac{3}{2}}\,\left| x\right| \,\sqrt{y}\]

function: determiningEquations(ode,y,x) returns the determining equations of the first order ODE $y' = \phi(x,y)$, which is given by: $\eta_x + \eta_y - \xi_x*\phi - \xi_y*\phi^2 - \xi*\phi_x - \eta*\phi_y$



In [16]:

    
determiningEquations('diff(y,x)=f(x)*g(y),y,x);









    Out[16]:





\[\tag{${\it \%o}_{173}$}-{\it \_eta}\,f\left(x\right)\,g_{\left(1\right)}(y)-\frac{d}{d\,y}\,{\it \_xi}\,f\left(x\right)^2\,g\left(y\right)^2-{\it \_xi}\,f_{\left(1\right)}(x)\,g\left(y\right)-\frac{d}{d\,x}\,{\it \_xi}\,f\left(x\right)\,g\left(y\right)+\frac{d}{d\,y}\,{\it \_eta}\,f\left(x\right)\,g\left(y\right)+\frac{d}{d\,x}\,{\it \_eta}\]

function: canoni(X,y,x) compute canonical coordinates from the input symmetry $X$.



In [17]:

    
canoni([1,0],y,x);









    Out[17]:





\[\tag{${\it \%o}_{174}$}\left[ r=y , s=x \right] \]

function: ode1FromSym(X,y,x) Generate the most general ODE under the group $X$. The result will contain unknown functions $\%F(x,y)$ in the independent and dependent variable.



In [18]:

    
ode1FromSym([0,y],y,x);









    Out[18]:





\[\tag{${\it \%o}_{175}$}\frac{d}{d\,x}\,y={\it \%F}\left(x\right)\,y\]

function: (expr) differentialInvariants(X,y,x) Compute differential invariants from the symmetry



In [19]:

    
differentialInvariants([1,0],y,x);









    Out[19]:





\[\tag{${\it \%o}_{176}$}\left[ y={\it \%c} \right] \]

function: (expr) riccati_RNF(ode,y,x) try to put the input Riccati ODE into it's Rational Normal Form.



In [20]:

    
riccati_RNF('diff(y,x)=x*y^2+5*x+5,y,x);









    





Riccati ODE: y' = a(x) + b(x)*y + c(x)*y^2, with [a(x),b(x),c(x)] =  \(5\,\left(x+1\right)\) ,  \(0\) ,  \(x\)






    Out[20]:





\[\tag{${\it \%o}_{177}$}\frac{d}{d\,t}\,u=\frac{5\,t^3\,u^2}{t+1}+\frac{10\,t^2\,u^2}{t+1}+\frac{5\,t\,u^2}{t+1}-\frac{u}{t+1}+\frac{t}{t+1}+\frac{1}{t+1}\]

Bibliography

[1] E.S. Cheb-Terrab and T. Kolokolnikov, First-order ordinary differential equations, symmetries and linear transformations, Euro. J. of Applied Mathematics 14 (2003)

[2] E.S. Cheb-Terrab and A.D. Roche, Symmetries and first order ODE patterns, Computer Physics Communications 113 (1998)

[3] E. Kamke, Differentialgleichungen, L$\ddot o$sungsmethoden und L$\ddot o$sungen, Leipzig (1959)

[4] E.S. Cheb-Terrab, L.G.S. Duarte and L.A.C.P. da Mota, Computer algebra solving of first order ODEs using symmetry methods, Computer Physics Communications 101 (1997)

[5] F. Schwarz, Symmetry analysis of Abel's equation, Studies in applied mathematics 100 (1998)

[6] F. Schwarz, Algorithmic solution of Abel's equation, Computing 61 (1998)