3.11 Obtain ODE in canonical coordinates (Lie symmetry)

How to obtain the ODE in canonical coordinates \(S(R)\) for an ode using Lie symmetry?

Here is a function I wrote with the help of Maple docs, that does that. For now, I am using first order ode. It takes the ode and the variables (coordinates) to use (\(R,S\)) is what normally used, and returns back the ODE in canonical coordinates. The ode should always be a quadrature for first order.

get_ODE_in_canonical:=proc(ode::`=`,y::symbol,x::symbol,S::symbol,R::symbol) 
  local infinitesimals,tr,itr,ODE; 
  infinitesimals := [DEtools:-symgen(ode)]; 
  if nops(infinitesimals)=0 then 
     RETURN(FAIL); 
  fi; 
  infinitesimals:=infinitesimals[1]; 
  tr  := DEtools:-canoni(infinitesimals,y(x),S(R)); 
  itr := op(1,[solve(tr,{x,y(x)})]); 
  ODE:= PDEtools:-dchange(itr,ode,[R,S(R)],simplify); 
  ODE:= op(solve(ODE,{diff(S(R),R)})); 
  RETURN(ODE); 
end proc:
 

Example

ODE := diff(y(x),x) = 1/x^2+y(x); 
get_ODE_in_canonical(ODE,y,x,S,R) 
 
    #diff(S(R), R) = 1/(exp(R)*R^2) 
 
ODE := diff(y(x),x) = x^2-2*x*y(x)+y(x)^2; 
get_ODE_in_canonical(ODE,y,x,S,R) 
   #diff(S(R), R) = 1
 

The following function is small change of the above. It accepts your choice of \(\xi ,\eta \) to use instead of the first one obtained by Maple. This makes it possible to experiment with different infinitesimals. Note that different infinitesimals will/can result in different canonical ODE, but the final solution for \(y(x)\) of course will always be the same after applying the transformation back to \(x,y\) coordinates.

restart; 
get_ODE_in_canonical_v2:=proc(ode::`=`,y::symbol,x::symbol,S::symbol,R::symbol,xi,eta) 
  local infinitesimals:=[xi,eta],tr,itr,ODE; 
  tr  := DEtools:-canoni(infinitesimals,y(x),S(R)); 
  itr := op(1,[solve(tr,{x,y(x)})]); 
  ODE:= PDEtools:-dchange(itr,ode,[R,S(R)],simplify); 
  ODE:= op(solve(ODE,{diff(S(R),R)})); 
  RETURN(simplify(ODE,symbolic)); 
end proc:
 

And now

ODE:=diff(y(x),x) = x^2-2*x*y(x)+y(x)^2; 
get_ODE_in_canonical_v2(ODE,y,x,S,R,1,1) 
 
    #diff(S(R), R) = 1/(R^2 - 1)