Notebook contains functions for weight-homogenous decomposition of the system of first order ordinary differential equations (ODEs).
LICENCE: GNU GPL v3 (http://www.gnu.org/copyleft/gpl.html)
NOTE: The programs were written by means of high standards, however, author is not responsible for any damages that result from the use of these programs. Use them at your own risk.

Functions:

NOTE: Execute this section before you do any calculations.

Helper functions:

Helper function. Check if given set is a subset of other set.

SubsetQ[set_,subset_]:=Intersection[set,subset]==subset;

Helper function. Add distinct components of vectors a and b (union of components and then sum them) with respect to the variables

AddDistinctComponents[a_,b_,variables_]:=Module[{monomialsA,monomialsB,union,vect },
monomialsA=Map[MonomialList[#,variables]&,a];
monomialsB=Map[MonomialList[#,variables]&,b];
union=MapThread[Union[#1,#2]&,{monomialsA,monomialsB}];
vect=Map[Total,union];
Return[vect];
];

Helper function. Join components of vectors of functions of given variables with the same weights. sol={ {vector1, weight1},{vector2,weight2},...}

JoinSameScaleinvariantTerms[sol_,variables_]:=Module[{tmp,positions},
(*triangular iteration and join the elements with the same weights*)
tmp=sol;
Do[
Do[ 
If[sol[[i,2]]==sol[[j,2]],
tmp[[i,1]]=AddDistinctComponents[tmp[[i,1]],tmp[[j,1]],variables];
] , 
{j,i+1,Length[sol]}
],
{i,1,Length[sol]}
];
(*mark element that have the same weights*)
positions={};
Do[
Do[ 
If[tmp[[i,2]]==tmp[[j,2]],
positions =Union[positions,{{j}}];
] , 
{j,i+1,Length[sol]}
],
{i,1,Length[sol]}
];
(*Remove repeated parts*)
tmp=Delete[tmp,positions];
Return[tmp];
];

Create terms which are special substitutions of other terms, e.g., (w1->1+w2) /. (w2->1)

MakeSpecialExpCases[allTerms_,variables_]:=Module[{tmp,term},
tmp=allTerms;
reduce[a_,b_]:={Map[#1[[1]]->(#1[[2]]//.b[[1]])&,a[[1]]]};
Do[
If[ Length[allTerms[[i,2,1]]] < Length[variables],
Do[
If[i==j,Continue[]];
term=allTerms[[i]];
term={#[[1]],reduce[#[[2]],allTerms[[j,2]]]}&@term;
tmp=Append[tmp,term];
,{j,1,Length[allTerms]}
];
];
,{i,1,Length[allTerms]}
];
Return[tmp];
];

Join terms with arbitrary exponents that fits,e.g.,(w1->1,w2->2) and (w3->3):

JoinExpTerms[allTerms_,variables_]:=Module[{all,tmp,wterm1,wterm2, w1, w2,wintersect},
all=allTerms;
tmp=all;
(*Iterate until all w's are joined*)
Do[
(*Iterate over all pairs of terms-triangular iteration*)
Do[
Do[
(*Check if exponents are orthogonal*)
wterm1=Flatten[all[[i,2]]];
wterm2=Flatten[all[[j,2]]];
w1=Map[#[[1]]&,wterm1];
w2=Map[#[[1]]&,wterm2];
wintersect=Map[#[[1]]&,Intersection[wterm1,wterm2]];
(*Join terms which have "orthogonal" terms*)If[Intersection[w1,w2]==wintersect,
tmp=Append[tmp,{AddDistinctComponents[all[[i,1]],all[[j,1]],variables],{Union[all[[i,2,1]],all[[j,2,1]]]}}]],
{j,i,Length[all]}],
{i,1,Length[all]}];
(*At the end join the same scale ivariant terms*)
all=JoinSameScaleinvariantTerms[tmp,variables];,
{p,Min[Map[Length[#[[2,1]]]&,allTerms]],Length[variables]}
];
Return[all];];

Jacobians:

JacobianMatrix[f_List?VectorQ,x_List]:=Outer[D,f,x]/;Equal@@(Dimensions/@{f,x});

JacobianDeterminant[f_List?VectorQ,x_List]:=Det[JacobianMatrix[f,x]]/;Equal@@(Dimensions/@{f,x});

Creates base vectors:

e[i_,n_]:=Module[{tmp},tmp=ConstantArray[0,n];tmp[[i]]=1;Return[tmp]];

Analoous to GroupBy[ sol,First] avaiable in Mathematica 10:

JoinSameWeightTerms[sol_]:=Module[{tmp,positions},
(*triangular iteration and join the elements with the same weights*)
tmp=sol;
Do[
Do[ 
If[sol[[i,1]]==sol[[j,1]],
tmp[[i,2]]=tmp[[i,2]]+tmp[[j,2]];
] , 
{j,i+1,Length[sol]}
],
{i,1,Length[sol]}
];
(*mark element that have the same weights*)
positions={};
Do[
Do[ 
If[tmp[[i,1]]==tmp[[j,1]],
positions =Union[positions,{{j}}];
] , 
{j,i+1,Length[sol]}
],
{i,1,Length[sol]}
];
(*Remove repeated parts*)
tmp=Delete[tmp,positions];
Return[tmp];
];

Core functions:

Print scale invariant terms:

PrintScaleInvariantTerms[result_]:=Module[{},
Print["('anything'-arbitrary value. Note that some special values can be given in the list of solutions.)"];
Map[Print["Term:", #[[1]]//MatrixForm, "; exponent:" , Flatten[#[[2]]]//MatrixForm]&,result];]

Print scale-invariant solutions:

PrintScaleInvariantSolutions[result_]:=Module[{},
Print["Every scale-invariant term provides solution of the form y(x)=((αx)^w). ('anything'-arbitrary value. Note that some special values can be given in the list of solutions.):"];
Map[Print["Term:", #[[1]]//MatrixForm, "; exponent:" , Flatten[#[[2]]]//MatrixForm,"; coefficients:",#[[3]]//MatrixForm]&,result];]

Perform searching of all scale-invariant  terms grouping them and then finds scale-invariant solutions. ‘anything’ means that the value can be arbitrary.

AnalyzeScaleInvariantTerms[system_,variables_]:=Module[{eps,monomials,monomials0,tup,combinations,powers, powSubst,subst,tsubst,exponents,eq,sol,allTerms,wterms,remainingw,addw,reduce},
If[ MemberQ[variables,t], Print["Use other variable than",t];Return[{}]];
If[Dimensions[system]≠ Dimensions[variables],Print["No. equations not the same as a number of varibales!"];Return[{}];];
monomials=Map[MonomialList[#,variables]&,system];
monomials0 = Map[Union[{0},#]&,monomials];
tup=Tuples[monomials0];
combinations=Complement[tup,{Table[0,{Length[monomials]}]}];
powers=Table[(Subscript[w, i]),{i,1,Length[variables]}];
powSubst=Map[(t^(#))&,powers];
subst=Flatten[MapThread[{#1->#2}&,{variables,powSubst}]];
tsubst=combinations/.subst;
exponents=Exponent[tsubst,t]/.{-∞->flaggedToDel};
eq=Map[MapThread[#1-1==#2&,{powers,#}]&,exponents];
sol=Map[Solve[#,powers]&,eq];
(*Combine terms with solutions*)
allTerms=MapThread[{#1,#2}&,{combinations,sol}];
allTerms=Select[allTerms, FreeQ[Flatten[#[[2]]],flaggedToDel]&];
(*Remove empty solutions*)
allTerms=Select[allTerms,#≠{_,{}}&];
(*Check if there is any solution*)
If[Min[Map[Length[#[[2,1]]]&,allTerms]]==0,
Print["No scale-invariant term found. Is the field linear?"];
Return[{}];
];
(*Simplify unevaluated terms by substitution*)
reduce[a_]:={Map[#1[[1]]->(#1[[2]]//.a[[1]])&,a[[1]]]};
Do[allTerms[[i,2]]=reduce[allTerms[[i,2]]],{i,1,Length[allTerms]}];
allTerms=JoinSameScaleinvariantTerms[allTerms,variables];
(*Make terms that are substitutions of other terms when some exponents are unevaluated*)
allTerms=MakeSpecialExpCases[allTerms,variables];
(*Join components that have arbitray w's *)
allTerms=JoinExpTerms[allTerms,variables];
(*Once again simplify unevaluated terms by substitution*)
(*reduce[a_]:={Map[#1[[1]]->(#1[[2]]//.a[[1]])&,a[[1]]]};*)
Do[allTerms[[i,2]]=reduce[allTerms[[i,2]]],{i,1,Length[allTerms]}];
allTerms=JoinSameScaleinvariantTerms[allTerms,variables];
(*Add exponents that are missing, i.e., arbitrary*)
Do[
wterms={};
wterms=Map[#[[1]]&,Flatten[allTerms[[i,2]] ]];
remainingw= Complement[powers,wterms];
addw=Map[#->anything&,remainingw];
allTerms[[i,2,1]]=Union[allTerms[[i,2,1]],addw],
{i,1,Length[allTerms]}
];
(*PrintScaleInvariantTerms[allTerms];*)
Return[allTerms];
]

Calculate scale-invariant Solutions:

CalculateScaleInvariantSolution[solution_,variables_]:=Module[{powers,allTerms,alphaSubst,alphas,alphaVectSubst,alphaEq,alphaSol,alphaSolSubst,nonemptyTerms, powersubst},
allTerms=sol;
alphas=Table[(Subscript[α, i]),{i,1,Length[variables]}];
alphaSubst= Flatten[MapThread[{#1->#2}&,{variables,alphas}]];
alphaVectSubst=allTerms//.alphaSubst;
powers=Table[(Subscript[w, i]),{i,1,Length[variables]}];
alphaEq=Map[MapThread[#1*#2==#3&,{alphas,powers,#[[1]]}]&,alphaVectSubst];
alphaEq=MapThread[#1/.#2[[2]]&,{alphaEq,allTerms}];
alphaSol=Map[Solve[Flatten[#],alphas]&,alphaEq];
allTerms=MapThread[{#1[[1]],#2[[2]],{#3}}&,{allTerms,allTerms,alphaSol}];
nonemptyTerms=Select[allTerms,#≠{_,{},_}&];
(*PrintScaleInvariantSolutions[allTerms];*)
Return[allTerms];
];

Removes irrelevant soluions(not describing singularities and too arbitrary):

RemoveNonsingular[solution_,variables_]:=Module[{sol,tmp, powers,alphas, subst,coeff},
sol=solution;
(*Delete exponets that do not describe singular solutions - positive*)
tmp={};
powers=Table[(Subscript[w, i]),{i,1,Length[variables]}];
Do[
subst=powers/.sol[[i,2]];
 If[
TrueQ[Fold[And,True,Flatten[Map[Positive,  #]]]&subst], 
tmp=Append[tmp,i]
],
{i,1,Length[sol]}
];
sol=Delete[sol,tmp];
(*Delete zero solutions*)
alphas=Table[(Subscript[α, i]),{i,1,Length[variables]}];
Do[
tmp={};
coeff=sol[[i,3]];
Do[
subst=Map[alphas/.#&,coeff];
If[  
Fold[And,True,Map[#==0&,subst[[1,j]]]],
tmp=Append[tmp,j]],
{j,1,Length[coeff[[1]]]}
];
sol[[i,3,1]]=Delete[sol[[i,3,1]],tmp],
{i,1,Length[sol]}
];
sol=Select[sol,#≠{_,_,{{}}}&];
Return[sol];
];

Remove the solutions which have arbitrary values of exponets

RemoveArbitrary[solution_,variables_]:=Module[{sol},
sol=solution;
(*Remove exponents with arbitrary values*)
sol=Select[sol, FreeQ[Flatten[#[[2]]],anything]&];
Return[sol];
];

Print the solutions obtained by KovalevskayaMatrix[]:

PrintScaleInvariantKovalevskaya[result_]:=Module[{},
Print["Every scale-invariant term provides solution of the form y(x)=((αx)^w). ('anything'-arbitrary value. Note that some special values can be given in the list of solutions.):"];
Do[
Print["Term:", #[[1,1]]//MatrixForm, "; exponents:" , Flatten[#[[1,2]]]//MatrixForm]&@  result[[i]];
MapThread[Print["-> coefficients:",#1//MatrixForm,", Kovalevskaya matrix:", #2//MatrixForm,", eigenvalues: ",#3//MatrixForm ]&,  {result[[i,1,3,1]],result[[i,2]],result[[i,3]] }];,
{i,1,Length[result]}];
];

Compute Kovalevskaya matrices for given scale invariant solutions and corresponding eigenvaues:

KovalevskayaMatrix[sol_,variables_]:=Module[{K,tmpK,powers,JacMat,alphas,alphaSubst,tmp,tmpEigen,eigenK},
K={};
eigenK={};
alphas=Table[(Subscript[α, i]),{i,1,Length[variables]}];
alphaSubst= Flatten[MapThread[{#1->#2}&,{variables,alphas}]];
Do[
tmpK={};
tmpEigen={};
powers=Table[(Subscript[w, i]),{i,1,Length[variables]}];
powers=Flatten[powers//.sol[[i,2]]];
JacMat=JacobianMatrix[sol[[i,1]],variables];
JacMat=JacMat//.alphaSubst;
tmpK =Map[(JacMat/.#)&,sol[[i,3]]];
tmpK=Map[#-DiagonalMatrix[powers]&,tmpK[[1]]];
tmpEigen=Map[Eigenvalues[#]&,tmpK];
K=Append[K,tmpK];
eigenK=Append[eigenK,tmpEigen];
,{i,1,Length[sol]}];
tmp = MapThread[{#1,#2,#3}&,{sol,K,eigenK}];
(*PrintScaleInvariantKovalevskaya[tmp];*)
Return[tmp];
];

Analyze equation according to given weights:

Decompose system of first order ODEs with variables according to self-similar solution number solNumber:

WeightHomogenousDecompose[solution_,system_,variables_,solNumber_]:=Module[{selfsimmonomials,monomials,powers,powSubst,powValues,subst,powersubst,monomialsSubst,exponents,qs,tmp,weightVectors},
If[solNumber>Length[solution],Print["No Solution"];Return[{}]];
(*get self-similar monomials*)
selfsimmonomials=Map[MonomialList[#,variables]&,sol[[solNumber,1]]];
(*get all monomials from the system*)
monomials=Map[MonomialList[#,variables]&,system];
monomials=MapThread[Complement[#1,#2]&,{monomials,selfsimmonomials}];
(*create power exponents*)
powers=Table[(Subscript[w, i]),{i,1,Length[variables]}];
powSubst=Map[(t^(#))&,powers];
(*weights*)
powValues=Flatten[powers/.sol[[solNumber,2]]];
(*Substitute x->(t^w) into equation and find weights of all terms*)
subst=Flatten[MapThread[{#1->#2}&,{variables,powSubst}]];
powersubst=subst/.sol[[solNumber,2]];
monomialsSubst=monomials//.powersubst;
exponents=Exponent[monomialsSubst,t];
(*calculate q's*)
qs=MapThread[#1-#2+1&,{exponents[[1]],powValues}];
tmp={};
Do[
tmp=Append[tmp,MapThread[{#1,#2}&,{qs[[i]],monomials[[i]]}]],
{i,1,Length[monomials]}];
(*Construct vectors*)
weightVectors={};
Do[
weightVectors=Append[weightVectors,Map[{#[[1]],#[[2]]*e[i,Length[tmp]]}&,tmp[[i]]]],
{i,1,Length[tmp]}];
weightVectors=Flatten[weightVectors,1];
(*Add terms wit the same weigts*)
weightVectors=JoinSameWeightTerms[weightVectors];
Print["Weight=",Flatten[#[[2]]]//MatrixForm,". Scale-invariant term =",#[[1]]//MatrixForm]&@solution[[solNumber]];
Map[Print["-> q = ",#[[1]],", term:",#[[2]]//MatrixForm]&,weightVectors];
Return[weightVectors];
];

Perform full analysis:

WeightHomogenousAnalysis[system_,variables_]:=Module[{sol,km},
sol=AnalyzeScaleInvariantTerms[system,variables];
sol=CalculateScaleInvariantSolution[sol,variables];
sol=RemoveNonsingular[sol,variables];
km=KovalevskayaMatrix[sol,variables];
PrintScaleInvariantKovalevskaya[km];
Print["Weight-homogenous decomposition:"];
Do[
WeightHomogenousDecompose[sol,system,variables,i],
{i,1,Length[sol]}];
];

Tests:

NOTE: The examples are based on the book: Alain Goriely “Integrability and Nonintegrability of Dynamical Systems“, World Scientific Publishing Company (2001)

Example 2.11:

equation={(Subscript[x, 3]),(Subscript[x, 4]),2*(Subscript[x, 1])*((Subscript[x, 2])^3)+(3/4)*((Subscript[x, 1])^3)*(Subscript[x, 2]),5*((Subscript[x, 2])^4)+3*((Subscript[x, 1])^2)*((Subscript[x, 2])^2)+(3/16)*((Subscript[x, 1])^4)};
variables={(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4])};
sol=AnalyzeScaleInvariantTerms[equation,variables];
PrintScaleInvariantTerms[sol]

sol=CalculateScaleInvariantSolution[sol,variables];
PrintScaleInvariantSolutions[sol];

sol=RemoveNonsingular[sol,variables];
PrintScaleInvariantSolutions[sol];

WeightHomogenousAnalysis[equation,variables]

Example 2.12:

equation={4y-(x^2)/2,x*y+(x^3)};
variables={x,y};
sol=AnalyzeScaleInvariantTerms[equation,variables];
PrintScaleInvariantTerms[sol]

sol=CalculateScaleInvariantSolution[sol,variables]
PrintScaleInvariantSolutions[sol];

Remove irrelevant soltuions:

sol=RemoveNonsingular[sol,variables];
PrintScaleInvariantSolutions[sol]

Kowalewskaya Matrix and eigenvalues:

km=KovalevskayaMatrix[sol,variables];
PrintScaleInvariantKovalevskaya[km]

Weight-homogenous decomposition:

WeightHomogenousDecompose[sol,equation,variables,1]

WeightHomogenousDecompose[sol,equation,variables,2]

Example 2.13:

equation={(Subscript[x, 2])*(((((Subscript[x, 1])^2)+((Subscript[x, 2])^2)))^2),-(Subscript[x, 1])*(((((Subscript[x, 1])^2)+((Subscript[x, 2])^2)))^2)};
variables={(Subscript[x, 1]),(Subscript[x, 2])};
sol=AnalyzeScaleInvariantTerms[equation,variables];
PrintScaleInvariantTerms[sol]

sol=CalculateScaleInvariantSolution[sol,variables];
PrintScaleInvariantSolutions[sol];

Remove irrelevant soltuions:

sol=RemoveNonsingular[sol,variables];
PrintScaleInvariantSolutions[sol]

Example 2.14(a cubic planar system):

equation={(Subscript[λ, 1])*(Subscript[x, 1])+(Subscript[x, 1])*((Subscript[a, 1])*(Subscript[x, 1])+(Subscript[a, 2])*(Subscript[x, 1])*(Subscript[x, 2])+(Subscript[a, 3])*(Subscript[x, 2])),(Subscript[λ, 2])*(Subscript[x, 2])+(Subscript[x, 2])*((Subscript[b, 1])*(Subscript[x, 1])+(Subscript[b, 2])*(Subscript[x, 1])*(Subscript[x, 2])+(Subscript[b, 3])*(Subscript[x, 2]))};
variables={(Subscript[x, 1]),(Subscript[x, 2])};
sol=AnalyzeScaleInvariantTerms[equation,variables];
PrintScaleInvariantTerms[sol];

Example 2.15(the Klapper-Rado-Tabor model):

equation={(Subscript[x, 4]),(Subscript[x, 5]),(Subscript[x, 6]),(1/3)*(Subscript[x, 1])*(Subscript[x, 3])-(β/3)*(Subscript[x, 1]),(1/3)*(Subscript[x, 2])*(Subscript[x, 3])+(β/3)*(Subscript[x, 2])-(γ/3)*(Subscript[x, 1]),((Subscript[x, 3])^2)+8*(Subscript[x, 5])-(β^2)};
variables=Table[(Subscript[x, i]),{i,1,6}];
sol=AnalyzeScaleInvariantTerms[equation,variables];
PrintScaleInvariantTerms[sol];

WeightHomogenousDecompose[sol,equation,variables,2]

Example 3.11(the planar Lotka-Volterra system):

equation={x(a*y-x),y*(b*x-y)};
variables={x,y};
sol=AnalyzeScaleInvariantTerms[equation,variables];
PrintScaleInvariantTerms[sol];

Calculate scale-invariant terms:

sol=CalculateScaleInvariantSolution[sol,variables];
PrintScaleInvariantSolutions[sol];

Remove irrelevant:

sol=RemoveNonsingular[sol,variables];
PrintScaleInvariantSolutions[sol]

PrintScaleInvariantSolutions[RemoveArbitrary[sol,variables]];

km=KovalevskayaMatrix[sol,variables];
PrintScaleInvariantKovalevskaya[km]

WeightHomogenousDecompose[sol,equation,variables,1]

WeightHomogenousDecompose[sol,equation,variables,2]

WeightHomogenousDecompose[sol,equation,variables,3]

Example 3.7  (the three wave interactions):

equation={-2*(y^2)+γ*x+z-(δ^2),2*x*y+γ*y+((γ*δ)/2),-2*z*x-2*z};
variables={x,y,z};
sol=AnalyzeScaleInvariantTerms[equation,variables];
PrintScaleInvariantTerms[sol]

Calculate scale-invariant terms:

sol=CalculateScaleInvariantSolution[sol,variables];
PrintScaleInvariantSolutions[sol];

sol=RemoveNonsingular[sol,variables];
PrintScaleInvariantSolutions[sol]

km=KovalevskayaMatrix[sol,variables];
PrintScaleInvariantKovalevskaya[km]

WeightHomogenousDecompose[sol,equation,variables,1]

WeightHomogenousDecompose[sol,equation,variables,2]

WeightHomogenousDecompose[sol,equation,variables,3]

WeightHomogenousDecompose[sol,equation,variables,4]

WeightHomogenousAnalysis[equation,variables];

Example 3.12(the Duffing oscillator):

equation={y,3*a*y+2*b*x-2*(x^3)};
variables={x,y};
sol=AnalyzeScaleInvariantTerms[equation,variables];
PrintScaleInvariantTerms[sol]

WeightHomogenousDecompose[sol,equation,variables,1]

Example 3.15(the Henon-Heiles system):

equation={(Subscript[x, 3]),(Subscript[x, 4]),-a*(Subscript[x, 1])-2*b*(Subscript[x, 1])*(Subscript[x, 2]),-(Subscript[x, 2])-b*((Subscript[x, 1])^2)+((Subscript[x, 2])^2)};
variables=Table[(Subscript[x, i]),{i,1,4}];
sol=AnalyzeScaleInvariantTerms[equation,variables];
PrintScaleInvariantTerms[sol]

Calculate scale-invariant terms:

sol=CalculateScaleInvariantSolution[sol,variables];
PrintScaleInvariantSolutions[sol];

sol=RemoveNonsingular[sol,variables];
sol=RemoveArbitrary[sol,variables];
PrintScaleInvariantSolutions[sol];

Remaining terms of the decomposition:

WeightHomogenousDecompose[sol,equation,variables,1]

Kovalevskaya matrix:

km=KovalevskayaMatrix[sol,variables];
PrintScaleInvariantKovalevskaya[km]