(* ::Package:: *)

(* :Name: DESymmetries` *)

(* :Title: Functions for Finding and Using Symmetries of Differential Equations *)

(* :Author: Karel Houfek *)

(* :Summary:
 This package provides several functions for finding and using symmetries of ODE or PDE or of their systems.
*)

(* :Package Version: 1.0 *)

(* :History:
 Version 1.0 - March 2011.
   - contains ChangeOfVariables
              CheckPointSymmetryOfDEs
              GetAnsatzPolynomial
 Version 1.1 - March 2013
   - added    GetConditionsForPointSymmetries
              ShowPointSymmetriesAndCommutationRelations
 *)

(* :Keywords:
 change of variables, ordinary differential equations, partial differential equations, point symmetries.
*)

(* :Mathematica Version: 7.0 *)

(* :Limitation:
*)

(* :Discussion:
*)

BeginPackage["DESymmetries`"];

ChangeOfVariables::usage = "ChangeOfVariables[pdes_, substitutions_, oldIndVar_, oldDepVar_, newIndVar_, newDepVar_] returns partial differential equations expressed in new variables."
CheckPointSymmetryOfDE::usage = "CheckPointSymmetryOfDE[pdes_, substitutions_, IndepVar_, DepVar_, xis_, etas_] checks whether a point transformation given by infinitesimals xis and etas is a point symmetry of a given system of PDEs."
GetConditionsForPointSymmetries::usage = "GetConditionsForPointSymmetries[zeroes_, colvars_] gets conditions from the infinitesimal criterium by returning coefficients of variables colvars in zeroes."
GetAnsatzPolynomial::usage = "GetAnsatzPolynomial[variables_, orders_, coef_, n_] constructs a general polynomial which can be of different orders in each variable but each monomial can be of maximal order n, const is the name of coefficients."
ShowPointSymmetriesAndCommutationRelations::usage = "ShowPointSymmetriesAndCommutationRelations[X_, f_, \[Epsilon]_, IndepVar_, DepVar_, xis_, etas_, par_, npar_] prints all independent infinitesimal generators and corresponding global transformations given by xis and etas, also prints a commutator table."



Unprotect[
  ChangeOfVariables,
  CheckPointSymmetryOfDE,
  GetConditionsForPointSymmetries,
  GetAnsatzPolynomial,
  ShowPointSymmetriesAndCommutationRelations
];

Begin["`Private`"];

Clear[
  ChangeOfVariables,
  CheckPointSymmetryOfDE,
  GetConditionsForPointSymmetries,
  GetAnsatzPolynomial,
  ShowPointSymmetriesAndCommutationRelations
]

ChangeOfVariables[pdes_, substitutions_, oldIndVar_, oldDepVar_, newIndVar_, newDepVar_, showAll_] :=
  Module[{order, nivar, ndvar, ivar, jac, invJac, subst, trans, allder, varDer},

  (* order of the system of the PDEs *)
  order = Max[Total /@ Cases[pdes, Derivative[d___][_][___]->{d}, Infinity]];
  If[order<1,
    Print["Warning: No derivatives found in the PDEs !"];
    Return[pdes /. subst // Simplify]
  ];
  If[showAll==1,Print["The highest order of derivative is ", order]];

  (* auxiliary variables - abbreviations *)
  nivar = Length[oldIndVar];
  ndvar = Length[oldDepVar];
  ivar = Delete[oldIndVar,0];

  (* calculate Jacobian and its inverse *)
  jac = Table[D[oldIndVar[[j]] /. substitutions, newIndVar[[i]]], {i,nivar}, {j,nivar}];
  invJac = Inverse[jac] // Simplify;
  If[showAll==1,Print["Jacobian (matrix A) is ", MatrixForm[jac], " and its inverse is ", MatrixForm[invJac]]];

  ndvar=Length[oldDepVar];
  (* trans is auxiliary array with dependent variables
    and their derivatives expressed in new coordinates *)
  Do[
    trans[oldDepVar[[j]],{}]=oldDepVar[[j]][ivar]/.substitutions, {j,ndvar}
  ];

  (* array of substitutions for derivatives *)
  subst={};
  Do[ (* loop over all orders of derivatives *)
    allder=Union[Subsets[Flatten[Table[oldIndVar,{i,1,p}]],{p}]]; (* all derivatives of a given order*)
    Do[ (* loop over all derivatives *)
      Do[ (* loop over dependent variables *)
        Do[ (* loop over independent variables *)
          varDer=Append[f,oldIndVar[[i]]];
          trans[oldDepVar[[j]],varDer] = Sum[invJac[[i,k]]D[trans[oldDepVar[[j]],f],newIndVar[[k]]],{k,nivar}]//Simplify;
          If[showAll==1,
            Print["Derivative of ", oldDepVar[[j]], " with respect to ", varDer, ":"];
            Print["   ", trans[oldDepVar[[j]],varDer]]
          ];
          AppendTo[subst,D[oldDepVar[[j]][ivar],Delete[varDer,0]]->trans[oldDepVar[[j]],varDer]],
          {i,nivar}
        ],
        {j,ndvar}
      ],
      {f,allder}
    ];
    (* get rid of repetitions which appear for higher derivatives *)
    Union[subst],
    {p,0,order-1}
  ];
  Return[pdes /. subst /. substitutions // Simplify]
]

CheckPointSymmetryOfDE[PDEs_, substitutions_, IndepVar_, DepVar_, xis_, etas_] :=
  Module[{order, solsubs, nivar, ndvar, ivar,
          der, subder, varder,
          usubst, backusubst,
          PDEs2, subs2, sol2,
          i, j, ii, jj, var, ni, zero,
          \[Eta], \[Xi]},

  (* order of the system of the PDEs *)
  order = Max[Total /@ Cases[PDEs, Derivative[d___][_][___]->{d}, Infinity]];
  If[order < 1,
    Print["Warning: No derivatives found in the PDEs !"];
    Return[{}]
  ];

  (* determine substitutions for given derivatives *)
  solsubs = Solve[PDEs == 0, substitutions];
  
  (* auxiliary variables - abbreviations *)
  nivar = Length[IndepVar];
  ndvar = Length[DepVar];
  ivar = Delete[IndepVar,0]; (* get rid of {} *)

  (* all combinations of independent variables with respect to which one has to differentiate,
     including zeroes derivative as the empty set *)
  der = Union[Subsets[Sort[Flatten[Table[IndepVar, {i,1,order}]]], order]];
  
  (* determination of infinitesimals of the extension of the point group *)
  Do[ (* loop over independent variables *)
    \[Xi][var] = xis[var],
    {var, IndepVar}
  ];
  Do[ (* loop over dependent variables *)
    \[Eta][var][{}] = etas[var],
    {var, DepVar}
  ];
  Do[ (* loop over dependent variables *)
    Do[ (* loop over all derivatives *)
      subder = j[[1;;Length[j]-1]];
      varder = Last[j];
      \[Eta][var][j]=D[\[Eta][var][subder],varder]-
      Sum[D[\[Xi][ii],varder] D[var[ivar],Delete[subder,0],ii],{ii,IndepVar}](*;
      Print["   \[Eta][",var,"][",j,"] = ",\[Eta][var][j]]*), 
      {j, der[[2;;]]} (* skip zeroth derivative, element {} *)
    ],
    {var, DepVar}
  ];

(* substitute auxiliary arrays like u[i] for derivatives in both PDEs and infinitesimals *)
  usubst = {};
  backusubst = {};
  Do[
    i = 0;
    Do[
      If[Length[j]==0,
         AppendTo[usubst,var[ivar]-> var[i]];AppendTo[backusubst,var[i]-> var[ivar]],
         AppendTo[usubst,D[var[ivar],Delete[j,0]]->var[i]];AppendTo[backusubst,var[i]->D[var[ivar],Delete[j,0]]]
      ];
      i++;,
      {j,der}
    ];
    ni[var]=i-1,
    {var,DepVar}
  ];
  PDEs2 = PDEs /. usubst;
  subs2 = substitutions /. usubst;
  Do[\[Xi][var] = \[Xi][var] /. usubst, {var, IndepVar}];
  Do[
    Do[
      \[Eta][var][jj] = \[Eta][var][jj] /. usubst,
      {jj, der}
    ],
    {var, DepVar}
  ];
(* Printing assignments of auxiliary arrays *)
  Do[
    ii = 0;
    Do[
      If[Length[j]==0,Print["   ",var[ii]," = ",var[ivar]],
        Print["   ",var[ii]," = ",D[var[ivar],Delete[j,0]]]
      ];
      ii++,
      {j,der}
    ],
    {var,DepVar}
  ];

  sol2 = Solve[PDEs2==0, subs2];

(* Infinitesimal criterion of invariance:
   Apply the prolongation of X to PDEs and substitute *)
  Do[
    jj = 0;
    zero[i] = Sum[\[Xi][ii] D[PDEs2[[i]],ii], {ii, IndepVar}];
    Do[
      zero[i] += Sum[\[Eta][var][j] D[PDEs2[[i]],var[jj]],{var,DepVar}];
      jj++,
      {j, der}
    ];
    zero[i] = zero[i] /. sol2 // Simplify,
    {i, Length[PDEs2]}
  ];

  Return[Table[Delete[zero[i],0],{i,1,Length[PDEs2]}]]
]

GetConditionsForPointSymmetries[zeroes_, colvars_] :=
  Module[{eqs},
  eqs = Flatten[CoefficientList[Table[
          Numerator[Cancel[Together[Expand[zeroes[[i]]]]]] // Simplify,
          {i,1,Length[zeroes]}], colvars]];
  Return[DeleteCases[Union[Table[Factor[eqs[[i]]],{i,1,Length[eqs]}]],0]]
];

GetAnsatzPolynomial[variables_, orders_, coef_, n_] :=
  Module[{auxarray, allterms, p, v, i},
  auxarray={};
  Do[
    AppendTo[auxarray,
      Table[variables[[v]],{i,1,orders[[v]]}]
    ],
    {v, 1, Length[variables]}
  ];
  allterms = Union[Subsets[Sort[Flatten[auxarray]], n]];
  p=0;
  Do[
    p += coef[i] Times[Delete[allterms[[i]], 0]],
    {i, 1, Length[allterms]}
  ];
  Return[p]
];

ShowPointSymmetriesAndCommutationRelations[X_, f_, \[Epsilon]_, IndepVar_, DepVar_, xis_, etas_, par_, npar_, addsub_] :=
  Module[{subst, subvar, var, allvar, \[Epsilon]subst, globtrans, eqs, com, lhs, sol, comXX, XX, ff, a, i, j, U, V},
  subvar = Table[var[Delete[IndepVar,0]] -> var, {var, DepVar}];
  allvar = Delete[Flatten[{IndepVar,DepVar}],0];
  subst = ConstantArray[0, npar];
  Do[
    subst[[j]] = {};
    Do[If[i===j,AppendTo[subst[[j]], par[i] -> 1],AppendTo[subst[[j]], par[i] -> 0]],{i,1,npar}],
    {j,1,npar}
  ];
  XX = Function[k, Function[y, (Total[Table[xis[var] D[y,var],{var,IndepVar}]] + 
                               Total[Table[etas[var] D[y,var],{var,DepVar}]])/.subst[[k]]/.addsub/.subvar]];
  Print["Infinitesimal operators:"];
  Do[Print[X[i],f[allvar]," = ",XX[i][f[allvar]]], {i,1,npar}];

  \[Epsilon]subst = Table[var -> var[\[Epsilon]],{var,{allvar}}];
  globtrans = {};
  Do[
    eqs = Flatten[{Table[D[var[\[Epsilon]],\[Epsilon]] == (xis[var]/.subst[[i]]/.addsub/.subvar/.\[Epsilon]subst), {var,IndepVar}],
                   Table[D[var[\[Epsilon]],\[Epsilon]] == (etas[var]/.subst[[i]]/.addsub/.subvar/.\[Epsilon]subst), {var,DepVar}],
                   Table[var[0] == ToString[var], {var,{allvar}}]}];
    (* Print[eqs]; *)
    AppendTo[globtrans, Flatten[Expand[DSolve[eqs,Table[var[\[Epsilon]],{var,{allvar}}],{\[Epsilon]}]]//Simplify]],
    {i,1,npar}
  ];
  Print["Corresponding global transformations:"];
  Do[Print[X[i]," gives ",globtrans[[i]]], {i,1,npar}];
  (* Print[Column[globtrans]]; *)

  com = Function[{U,V},Function[ff,U[V[ff]]-V[U[ff]]//Simplify]];
  Do[
    Do[
      lhs = DeleteCases[Flatten[CoefficientList[Expand[com[XX[i],XX[j]][f[allvar]]-Sum[a[k] XX[k][f[allvar]],{k,1,npar}]],
              Flatten[{Table[D[f[allvar],var],{var,{allvar}}],allvar}]]],0];
      sol = Solve[lhs==0];
      comXX[i,j] = Sum[a[k] X[k],{k,1,npar}]/.sol,
      {j,1,npar}
    ],
    {i,1,npar}
  ];
  Print["Commutator table:"];
  Print[TableForm[Array[comXX,{npar,npar}],TableHeadings->{Automatic,Automatic}]];
  Return[Null];
];


End[ ]; (* End `Private` Context. *)

SetAttributes[
{ ChangeOfVariables,
  CheckPointSymmetryOfDE,
  GetConditionsForPointSymmetries,
  GetAnsatzPolynomial,
  ShowPointSymmetriesAndCommutationRelations },
{ Protected, ReadProtected }
];

EndPackage[ ]; (* End package Context. *)

