In this example we use the package IntroToSymmetry.m to work out the point group of the 3-D incompressible Navier-Stokes equations,

Ut+U*Ux+V*Uy+W*Uz+Px-nu*Uxx-nu*Uyy-nu*Uzz=0,
Vt+U*Vx+V*Vy+W*Vz+Py-nu*Vxx-nu*Vyy-nu*Vzz=0,
Wt+U*Wx+V*Wy+W*Wz+Py-nu*Wxx-nu*Wyy-nu*Wzz=0
Ux+Vy+Wz=0.
 
P is the kinematic pressure, pressure/density and nu is the kinematic viscosity.

This next command turns off spurious spelling error warnings.

Off[General::spell]

Clear all symbols in the current context.

ClearAll[Evaluate[Context[]<>"*"]]

First read in the package which is located in User Home Folder/Library/Mathematica/Applications/SymmetryAnalysis.

Needs["SymmetryAnalysis`IntroToSymmetry`"]

Enter the list of independent variables.

independentvariables={"x","y","z","t"};

Enter the list of dependent variables.

dependentvariables={"U","V","W","P"};

Enter the list of function and/or constant names that need to be preserved when the equation is converted to generic variables.

frozennames={"nu"};

Enter the maximum derivative order of the input equation(s).

porder=2;

The maximum derivative order that the infinitesimals are assumed to depend on is specified by the input parameter rorder. This parameter is only nonzero when the user is looking for Lie contact groups or Lie-Backlund groups. For the usual case where one is searching for point groups set rorder=0.

rorder=0;

When searching for Lie-Backlund groups (rorder=1 or greater) one can, without loss of generality, leave the independent variables untransformed. The corresponding infinitesimals (the xse's) are set to zero by setting xseon=0. If one is searching for point groups then set xseon=1. The choice xseon=1 is also an option when looking for Lie-Backlund groups and this can be useful when looking for contact symmetries.

xseon=1;

When searching for Lie-Backlund groups it is necessary to differentiate the input equations with respect to all the independent variables producing derivatives of order p+r. These higher order differential consequences are appended to the set of rules applied to the invariance condition. This process is carried out automatically when internalrules=1. For point groups the equation or equation system is the only rule or set of rules needed and one sets internalrules=0.

internalrules=0;

Construct the array of rules to be applied to the invariance condition.

rulesarray={
"D[U[x,y,z,t],t]->-
U[x,y,z,t]*D[U[x,y,z,t],x]-
V[x,y,z,t]*D[U[x,y,z,t],y]-
W[x,y,z,t]*D[U[x,y,z,t],z]-
D[P[x,y,z,t],x]+
nu*D[U[x,y,z,t],x,x]+
nu*D[U[x,y,z,t],y,y]+
nu*D[U[x,y,z,t],z,z]",
"D[V[x,y,z,t],t]->-
U[x,y,z,t]*D[V[x,y,z,t],x]-
V[x,y,z,t]*D[V[x,y,z,t],y]-
W[x,y,z,t]*D[V[x,y,z,t],z]-
D[P[x,y,z,t],y]+
nu*D[V[x,y,z,t],x,x]+
nu*D[V[x,y,z,t],y,y]+
nu*D[V[x,y,z,t],z,z]",
"D[W[x,y,z,t],t]->-
U[x,y,z,t]*D[W[x,y,z,t],x]-
V[x,y,z,t]*D[W[x,y,z,t],y]-
W[x,y,z,t]*D[W[x,y,z,t],z]-
D[P[x,y,z,t],z]+
nu*D[W[x,y,z,t],x,x]+
nu*D[W[x,y,z,t],y,y]+
nu*D[W[x,y,z,t],z,z]",
"D[U[x,y,z,t],x]->-D[V[x,y,z,t],y]-D[W[x,y,z,t],z]"};

Enter the x-momentum equation as a string.

inputequation1=
"D[U[x,y,z,t],t]+
        U[x,y,z,t]*D[U[x,y,z,t],x]+
    V[x,y,z,t]*D[U[x,y,z,t],y]+
    W[x,y,z,t]*D[U[x,y,z,t],z]+
    D[P[x,y,z,t],x]-
      nu*D[U[x,y,z,t],x,x]-
  nu*D[U[x,y,z,t],y,y]-
  nu*D[U[x,y,z,t],z,z]";

Enter the y-momentum equation as a string.

inputequation2=
"D[V[x,y,z,t],t]+
        U[x,y,z,t]*D[V[x,y,z,t],x]+
    V[x,y,z,t]*D[V[x,y,z,t],y]+
    W[x,y,z,t]*D[V[x,y,z,t],z]+
    D[P[x,y,z,t],y]-
      nu*D[V[x,y,z,t],x,x]-
  nu*D[V[x,y,z,t],y,y]-
  nu*D[V[x,y,z,t],z,z]";

Enter the z-momentum equation as a string.

inputequation3=
"D[W[x,y,z,t],t]+
        U[x,y,z,t]*D[W[x,y,z,t],x]+
    V[x,y,z,t]*D[W[x,y,z,t],y]+
    W[x,y,z,t]*D[W[x,y,z,t],z]+
    D[P[x,y,z,t],z]-
      nu*D[W[x,y,z,t],x,x]-
  nu*D[W[x,y,z,t],y,y]-
  nu*D[W[x,y,z,t],z,z]";

Enter the continuity equation as a string.

inputequation4=
"D[U[x,y,z,t],x]+D[V[x,y,z,t],y]+D[W[x,y,z,t],z]";

Now work out the determining equations of the Lie point group which leaves equation1 invariant.

Timing[FindDeterminingEquations[
independentvariables,dependentvariables,frozennames,porder,rorder,xseon,inputequation1,rulesarray,internalrules]]

The determining equations of equation1.

zdeterminingequations1=zdeterminingequations;

Length[zdeterminingequations1]

Now work out the determining equations of the Lie point group which leaves equation2 invariant.

Timing[FindDeterminingEquations[
independentvariables,dependentvariables,frozennames,porder,rorder,xseon,
inputequation2,rulesarray,internalrules]]

The determiningequations of equation2.

zdeterminingequations2=zdeterminingequations;

Length[zdeterminingequations2]

Now work out the determining equations of the Lie point group which leaves equation3 invariant.

Timing[FindDeterminingEquations[
independentvariables,dependentvariables,frozennames,porder,rorder,xseon,
inputequation3,rulesarray,internalrules]]

The determiningequations of equation3.

zdeterminingequations3=zdeterminingequations;

Length[zdeterminingequations3]

Now work out the determining equations of the Lie point group which leaves equation3 invariant.

Timing[FindDeterminingEquations[
independentvariables,dependentvariables,frozennames,porder,rorder,xseon,
inputequation4,rulesarray,internalrules]]

The determiningequations of equation3.

zdeterminingequations4=zdeterminingequations;

Length[zdeterminingequations4]

zdeterminingequations5=
Join[zdeterminingequations1,
zdeterminingequations2,zdeterminingequations3,zdeterminingequations4];

Length[zdeterminingequations5]

Here is the correspondence between z-variables and conventional variables.

ztableofrules

Now solve the determining equations in terms of multivariable polynomials of some selected order. The Mathematica function Solve uses Gaussian elimination to solve a large number of linear equations for the polynomial coefficients. The time roughly follows

time/timeref=((number of equations)/(number of equationsref))^n

where the exponent is between 2.4 and 2.7. The Mathematica function Timing outputs the time required for the SolveDeterminingEquations function to execute. Running the notebook for, say, order 1and order 2 can provide a reasonably accurate estimate for the time required to run for third and higher order.

Timing[SolveDeterminingEquations[
independentvariables,dependentvariables,rorder,xseon,zdeterminingequations5,order=3]]

Here are the infinitesimal transformation functions for the independent variables.

xsefunctions

and for the dependent variables.

etafunctions

infinitesimalgroups1=infinitesimalgroups/.
{z1->x,z2->y,z3->z,z4->t,z5->U,z6->V,z7->W,z8->P};

Check the groups. xse's are on the left, eta's are on the right.

Column[infinitesimalgroups1]

Check the infinitesimals in the determining equations.

Simplify[ToExpression[zdeterminingequations]]

The equation admits one dilation group. Note the increasing powers of time that are evidence of arbitrary time dependent translations of the spatial coordinates. In addition there is an arbitrary function of time that can be added to the pressure. Lets join these groups with the rest and then form the commutator table of the remaining groups.

translationgroups=
{{{f1[t],0,0,0},{D[f1[t],t],0,0,-x*D[f1[t],t,t]}},
{{0,f2[t],0,0},{0,D[f2[t],t],0,-y*D[f2[t],t,t]}},
{{0,0,f3[t],0},{0,0,D[f3[t],t],-z*D[f3[t],t,t]}},
{{0,0,0,0},{0,0,0,g[t]}}}

infinitesimalgroups2=Drop[infinitesimalgroups1,{8,21}]

infinitesimalgroups3=Join[infinitesimalgroups2,translationgroups]

Column[infinitesimalgroups3]

Create the commutator table.

MakeCommutatorTable[
independentvariables,dependentvariables,infinitesimalgroups3]

MatrixForm[commutatortable]

MaxMemoryUsed[ ]