findSym.nb

(*
This notbook computes the determining equations for the lie \
symmetries of an evolutionary PDE $u_t=K[u,u_x,u_{xx},...]$.
Equation numbers in documation refer to \
http://www.math.umn.edu/~olver/sm_
*)


ClearAll["Global`*"]
SetDirectory[NotebookDirectory[]]

"/Users/E/Dropbox/RESEARCH/misc"


(*K is the RHS of the PDE*)
(*cond is the corresponding infinitesimal \
symmerty condition*)

(*Example: the heat equation*)
cond[Q_] := D[Q, {x, 2}];
K[Q_] := D[Q, {x, 2}];


(*Example: the heat equation*)

cond[Q_] := D[Q, {x, 2}] + 2 D[u[x, t], x] D[Q, x];
K[Q_] := D[Q, {x, 2}] + D[u[x, t], x]^2;

(*Example: aggregation equation*)

cond[Q_] := D[Q, {x, 2}] - 4 a u[x, t] Q;
K[Q_] := D[Q, {x, 2}] - 2 a Q^2;


(*The characteristic (3.8)*)

Q = \[Phi][x, t, 
    u[x, t]] - \[Tau][x, t, u[x, t]] D[u[x, t], t] - \[Xi][x, t, 
     u[x, t]] D[u[x, t], x];
(*the maximum derivatives to compute*)
maxderivs = 4;
J = Flatten[Table[{#, i - #} & /@ Range[0, i], {i, 0, maxderivs}], 
   1];
Jderivs = {{x, #[[1]]}, {t, #[[2]]}} & /@ J;

(*defining Dx allows us to treat the derivative of Q as a mathematica \
function*)
 Dx[Q_] := D[Q, x];
(*we need to mathematica treat the derivatives D[u[x,t],t] as \
indepndent variables*)
(*10 should be big enough*)

uxReplace = (D[u[x, t], {x, #}] -> 
      ToExpression["u" <> ConstantArray["x", #]]) & /@ 
   Range[0, maxderivs];
(*now we use the original equation to obtain only partial derivates \
with respect to x*)
(*so replace replaces all partial derivatives \
D[u[x,t],?] with the correspoding ux..x type terms*)
(*MAKE SURE \
#[[,]] MATCH ORDER K, Dx ORDER (????)*)

replace = (D[u[x, t], #[[1]], #[[2]]] -> 
      Simplify[
       Nest[K, Nest[Dx, u[x, t], #[[1, 2]]], #[[2, 2]]] /. 
        uxReplace]) & /@ Jderivs;

(*this line actually outputs the terms resulting from enacting the \
infintesimal generator in evolutionary form on the PDE according to \
(3.31)*)
(*the rest is just making a nice looking table*)
termsList = ({#} /. {Times -> 
         List} & /@ (Expand[(D[Q, t] - cond[Q]) /. 
         replace] /. {Plus -> List})) //. {x_List} -> x;


(* make sure every list starts with an integer *)

padded = If[IntegerQ[#[[1]]], #, Prepend[#, 1]] & /@ termsList;
(* pad each list on the right so they are all the same length*)

expanded = 
  PadRight[Delete[Prepend[#, Last[#]], -1], 10, 1] & /@ padded;
(*get the coefficients*)

coeffs = (# /. {List -> Times}) & /@ 
   Transpose[Transpose[expanded][[1 ;; 2]]];
(*get the monomials*)

monos = (# /. {List -> Times}) & /@ 
   Transpose[Delete[Transpose[expanded], {{1}, {2}}]];
(*put back together*)
out = Transpose[{monos, coeffs}];
(*generate determining equations *)

reduced2 = {#, Total[Transpose[Cases[out, {#, _}]][[2]]]} & /@ 
   DeleteDuplicates[monos];
TableForm[reduced2]