Helmholtz decomposition and potential functions for n-dimensional analytic vector fields

Erhard Glötzl, Oliver Richters
Supplementary Material for the article published in: Journal of Mathematical Analysis and Applications, 2023.

This worksheet allows to calculate a Helmholtz decomposition for various n-dimensional vector fields. First, the functions are defined, then Theorem 6.1 and several applications are implemented.
Version 2, February 2023.

Potential updates will be published at https://www.oliver-richters.de/helmholtz and https://doi.org/10.5281/zenodo.7512798

Functions

getAllVariables

headlist={Or,And,Equal,Unequal,Less,LessEqual,Greater,GreaterEqual,Inequality};
getAllVariables[f_?NumericQ]:=List[]
getAllVariables[{}]:=List[]
getAllVariables[t_]/;MemberQ[headlist,t]:=List[]
getAllVariables[ll_List]:=Flatten[Union[Map[getAllVariables[#]&,ll]]]
getAllVariables[Derivative[n_Integer][f_][arg__]]:=getAllVariables[{arg}]
getAllVariables[f_Symbol[arg__]]:=Module[{fvars},If[MemberQ[Attributes[f],NumericFunction]||MemberQ[headlist,f],fvars=getAllVariables[{arg}],(*else*)fvars=f[arg]];
fvars]
getAllVariables[other_]:=other
getUniqueVariables[all_] := Sort[DeleteDuplicates[Flatten[List[getAllVariables[all]]]]]
On[Assert];
Assert[getUniqueVariables[(Subscript[x, 2])(Subscript[x, 3])] == {(Subscript[x, 2]), (Subscript[x, 3])}];
Assert[getUniqueVariables[Exp[(Subscript[x, 2])(Subscript[x, 3])]] == {(Subscript[x, 2]), (Subscript[x, 3])}];
Assert[getUniqueVariables[Sin[(Subscript[x, 3])]Exp[(Subscript[x, 3])(Subscript[x, 2])]] == {(Subscript[x, 2]), (Subscript[x, 3])}];

from https://mathematica.stackexchange.com/questions/21257/extracting-variables-from-an-expression

Integrals and Incomplete Laplacians

(* Riemann-Liouville Integral, integrate func p-times with respect to (Subscript[x, m]) *)
IntegrateXmP[func_,p_,m_]:=Switch[p,0,func,_,Integrate[IntegrateXmP[func,p-1,m],{(Subscript[x, m]),0,(Subscript[x, m])}]];

(* Incomplete Laplacian: (Subscript[Δ, \m]) *)
LaplacianIncomplete[func_, m_]:=Laplacian[func, xvec] - D[func, {(Subscript[x, m]), 2}];

(* Incomplete Laplacian to the power of p: ((Subscript[Δ, \m])^p) *)
LaplacianIncompleteP[func_, p_, m_]:=Switch[p,0,func,_,LaplacianIncomplete[LaplacianIncompleteP[func, p-1, m], m]];

Calculate Potential Matrix (Subscript[F, ij])

(* === CalculateFij returns Fij with all the inputs manually given === *)
CalculateFij[f_,k_,m_,λ_,u_,W_]:=Module[{Fij},
(* Check whether conditions of Theorem 6.1 are fulfilled *)
On[Assert];
Assert[Length[f]==Length[xvec]==dim, "Length[f], Length[xvec] and dim must be identical"];
Assert[1 ≤ k ≤ dim, "k must be between 1 and dim"];
Assert[1 ≤ m ≤ dim, "m must be between 1 and dim"];
Assert[u≠0 || u==1-"a"^2/"w"^2 (*added to avoid assertion errors*), "u must be non-zero"];
If[k ≠ m, For[i=1,i≤dim,i++,Assert[ D[f[[k]],(Subscript[x, i])]===0 || i==k || i==m]]];
Assert[FullSimplify[D[W,{(Subscript[x, m]),2λ}] - f[[k]]] == 0,
    StringJoin["Differentiating W=,",W," with respect to (Subscript[x, m]) with m=",m," 2λ-times with λ=,",λ, " must yield (Subscript[f, k])=",f[[k]]]];
Assert[FullSimplify[(-1)^λ LaplacianIncompleteP[W,λ,m] - (1-u)f[[k]] == 0], "Terminal condition not fulfilled"];
(* Calculate Fij *)
Fij = Simplify[Table[KroneckerDelta[k,i]*Sum[ (-1)^p/u*D[ LaplacianIncompleteP[W, p, m],(Subscript[x, j]), {(Subscript[x, m]),2λ-2p-2} ],{p,0,λ-1}] ,{i,dim},{j,dim}]];
Return[Fij]
]; (* end of function CalculateFij *)

(* === AutocalculateSingleFij returns Fij given a simple expression (no sums!) M and the vector coordinate k === *)
(* === Don't call this function directly, use AutocalculateFij[f] === *)
AutocalculateSingleFij[S_, k_] := Module[{i,j,Singlef,xiWithoutxk,Q1,Q2,Q3,Q4,m,λ,u,W,v1,v2,Fij},
If[S === 0, Return[Table[0,{i,dim},{j,dim}]]];
Singlef = Table[S KroneckerDelta[i,k],{i,dim}];
xiWithoutxk = Intersection[xvec,DeleteCases[getUniqueVariables[S], (Subscript[x, k])]];
(* find suitable monomial to differentiate *)
Q1 = PolynomialQ[S,xiWithoutxk];
Q2 = PolynomialQ[S,xvec] && Length[xiWithoutxk] == 1
       && Exponent[S,First[xiWithoutxk]] > Exponent[S,(Subscript[x, k])];
Q3 = PolynomialQ[S,(Subscript[x, k])] && Length[xiWithoutxk] == 1;
If[Q1 || Q3,
 m = If[Q1 && ! Q2, k, (* i *) First[Flatten[Position[xvec, First[xiWithoutxk]]]]];
 u = 1;
 λ = Ceiling[(Total[Exponent[S,xvec]] - Exponent[S,(Subscript[x, m])] + 1)/2]; (* nochmal checken ob M ein polynom ist *)
 W = IntegrateXmP[S, 2λ, m];
 If[PrintInfo, Print[S," in f",k," is partly monomial, use m=",m,", λ=",λ]];
 Fij = CalculateFij[Singlef, k, m, λ, u, W];
 Return[Fij]
];
(* try to catch functions exp, sin, cos *)
v1 = FullSimplify[D[S,{(Subscript[x, k]),2}]/S];
v2 = FullSimplify[LaplacianIncomplete[S, k]/S];
xiInv1 = Intersection[xvec, getUniqueVariables[v1]];
xiInv2 = Intersection[xvec, getUniqueVariables[v2]];
Q4 = ! v1 === 0 && Length[xiInv1] == 0 && Length[xiInv2] == 0;
If[Q4,
 λ = 1;
 m = k;
 W = S/v1;
 u = 1 + v2/v1;
 If[PrintInfo,Print[S," in f",k," can be solved with m=",m,", λ=",λ,", W=",W,", u=",u]];
 Fij = CalculateFij[Singlef, k, m, λ, u, W];
 Return[Fij]
];
Print["No decomposition found for expression ", S, " in f", k,"."];
Abort[];
]; (* end of function AutocalculateSingleFij *)

(* === AutocalculateFij[f] returns Fij for a given vector field f, if possible === *)
AutocalculateFij[f_] := Module[{M,Fij,Poly,PrintInfo,pol,k,i,j},
PrintInfo=False; (* switch verbose printing on / off here *)
Fij = Table[0,{i,dim},{j,dim}]; (* initialize with 0 *)
Do[
  Poly = Expand[f[[k]]];
  Poly = If[Head[Poly]===Plus, List@@Poly, List[Poly]]; (* Make sure Poly is a proper list, even if f[[k]] has only one summand *)
  If[PrintInfo, Print["k=",k,": ",Poly]];
  Do[
      M = Poly[[pol]];
      Fij = Fij + AutocalculateSingleFij[M,k];
   , {pol,Length[Poly]}
   ];
, {k,dim}
];
Return[Fij]
]; (* end of function AutocalculateFij *)

ROTij[v_]:=FullSimplify[Table[ D[v[[i]],(Subscript[x, j])]-D[v[[j]],(Subscript[x, i])],{i,dim},{j,dim}]];
GF[Fij_]:=FullSimplify[Sum[Fij[[k,k]], {k,dim}]];
RijF[Fij_]:=FullSimplify[Fij-Transpose[Fij]];
gF[Fij_]:=FullSimplify[Grad[GF[Fij],xvec]];
rF[Fij_]:=FullSimplify[Table[Sum[D[RijF[Fij][[j,i]],(Subscript[x, i])],{i,dim}],{j,dim}]];

PrintHelmholtz

PrintHelmholtz[f_,Fij_,SF_]:=(
PrintCurl=False;
Print["f = ",MatrixForm[f]];
G=GF[Fij];
Rij=RijF[Fij];
g=gF[Fij];
r=rF[Fij];
Print["Div f = ",Div[f,xvec]];
If[dim<4 && PrintCurl,Print["Curl f = ",Curl[f,xvec]]];
Print[OverBar["ROT"]," f = ",If[MemberQ[SF,"ROTf"],ROTij[f],MatrixForm[ROTij[f], TableSpacing->{1,3}]]];
Print["Fij = ",If[MemberQ[SF,"Fij"],Fij,MatrixForm[Fij, TableSpacing->{1,4}]]];
Print["G = ",G];
Print["Rij = ",If[MemberQ[SF,"Rij"],Rij,MatrixForm[Rij, TableSpacing->{1,4}]]];
Print["g = ",If[MemberQ[SF,"g"],g,MatrixForm[g]]];
Print["r = ",If[MemberQ[SF,"r"],r,MatrixForm[r]]];
If[dim<4 && PrintCurl,Print["Curl g = ",Curl[g,xvec]]];
Print[OverBar["ROT"]," g = ",MatrixForm[ROTij[g]]];
Print["Div r = ",FullSimplify[Div[r,xvec]]];
Assert[FullSimplify[f - g-r] === Array[0 &, dim], "f not equal g + r"];
Print["g + r = ", MatrixForm[Simplify[g+r]]];
Print["f = g + r: ",FullSimplify[Expand[f]]==FullSimplify[Expand[g+r]]];
)
PrintHelmholtz[f_,Fij_]:=(AutocalculateFij[f]
If[SubsetQ[{"ROTf","Fij","Rij","g","r"},Fij],
 PrintHelmholtz[f,AutocalculateFij[f],Fij],
 PrintHelmholtz[f,Fij,{}]
];
)
PrintHelmholtz[f_]:=(
On[Assert];
Assert[Length[f]==Length[xvec]==dim];
PrintHelmholtz[f,AutocalculateFij[f]];
)

Theorem 6.1 and its applications

You can adapt the following formulas up to a certain extent. The script will tell you if it cannot solve the task.
PrintHelmholtz can be called
- with one argument: the vector field
- with two arguments: the vector field, the potential matrix Fij
- with two arguments: the vector field, a list {“ROTf”,”Fij”,”Rij”,”g”,”r”} that selects StandardForm as output for these variables
- with three arguments: vector field, potential matrix, a list where to use StandardForm

Theorem 6.1

dim=3;     (* <-- dimension of vector *)
k=1;         (* <-- index of non-zero component *)
m=3;         (* <-- coordinate to be integrated *)
xvec = Array[Subscript[x,#]&, dim]; (* define vector *)

fk=(Subscript[x, k])^2 (Subscript[x, m]);      (* <-- if m ≠ k, it must depend only on xk and xm *)
f=Table[fk KroneckerDelta[i,k],{i,dim}];
u=1;   (* non-zero real number *)
λ=2;   (* integer such that: (∂_xm)^(2λ) W = fk, and (-1)^λ (Subscript[Δ, \xm])^λ W == (1-u) fk *)
W=IntegrateXmP[fk,2λ,m]; (* you can also specify a different function here *)
Fij = CalculateFij[f,k,m,λ,u,W];
PrintHelmholtz[f,Fij]

Example: f = ((e^(Subscript[x, 1])) (e^(Subscript[x, 2])), 0)

dim=2; (* <-- dimension *)
xvec = Array[Subscript[x,#]&, dim]; (* define (Subscript[x, i]) *)
f={Exp[(Subscript[x, 1])]Exp[(Subscript[x, 2])],0};
PrintHelmholtz[f]

Example on linear vector field: f = L x

dim = 4;
xvec = Array[Subscript[x,#]&, dim];
Lij = Array[Subscript["L", #1, #2] &,{dim,dim}];
f = Table[Sum[Lij[[i,j]] (Subscript[x, j]),{j,1,dim}],{i,1,dim}];
PrintHelmholtz[f]

Example on 2-D linear vector field: f = L x

dim = 2;
xvec = Array[Subscript[x,#]&, dim];
Lij={{1,1},{-1,1}};
Print["Lij = ",MatrixForm[Lij,TableSpacing->{1,4}]];
f = Table[Sum[Lij[[i,j]] (Subscript[x, j]),{j,dim}],{i,dim}];
PrintHelmholtz[f]

Example 1: f = [((Subscript[x, 1])^2) ((Subscript[x, 2])^20), 0], short, integrating (Subscript[x, 2])

dim=2; (* <-- dimension *)
xvec = Array[Subscript[x,#]&, dim]; (* define (Subscript[x, i]) *)
f={(Subscript[x, 1])^2 (Subscript[x, 2])^100,0};
PrintHelmholtz[f]

Example 1: f = [((Subscript[x, 1])^2) ((Subscript[x, 2])^20), 0], long, integrating (Subscript[x, 1])

dim=2; (* <-- dimension *)
xvec = Array[Subscript[x,#]&, dim]; (* define (Subscript[x, i]) *)
f={(Subscript[x, 1])^2 (Subscript[x, 2])^20,0};
λ=Exponent[f[[1]],(Subscript[x, 2])]+1;
W=IntegrateXmP[f[[1]],2λ,1];
PrintHelmholtz[f,CalculateFij[f,1,1,λ,1,W],{"Fij","Rij","g","r"}]

Example 1: f = [((Subscript[x, 1])^2) ((Subscript[x, 2])^20) ((Subscript[x, 3])^2), 0], long, integrating (Subscript[x, 1])

dim=3; (* <-- dimension *)
xvec = Array[Subscript[x,#]&, dim]; (* define (Subscript[x, i]) *)
f={(Subscript[x, 1])^2 (Subscript[x, 2])^20 (Subscript[x, 3])^2,0,0};
PrintHelmholtz[f,{"Rij","Fij","g","r"}]

Example 2: f = [cos(w (Subscript[x, 1])) exp(a (Subscript[x, 2])), 0]

dim=2; (* <-- dimension *)
xvec = Array[Subscript[x,#]&, dim]; (* define (Subscript[x, i]) *)
f={Cos["w" (Subscript[x, 1])]Exp["a" (Subscript[x, 2])],0} ; (* assuming (a^2) ≠ (w^2) *)
PrintHelmholtz[f]

Complicated example with Sinus and Exp

dim=3; (* <-- dimension *)
xvec = Array[Subscript[x,#]&, dim]; (* AutocalculateFij[f]define (Subscript[x, i]) *)
f={Exp[2(Subscript[x, 1])+ (Subscript[x, 2])](Cos[(Subscript[x, 3])]Sin[(Subscript[x, 3])]+1)+((Subscript[x, 1])+(Subscript[x, 2]))(Subscript[x, 3]),(Subscript[x, 3])^4,(Subscript[x, 2])^3 (Subscript[x, 3])^4+(Subscript[x, 1])Sin[(Subscript[x, 3])]Cos[(Subscript[x, 2])]}; 
PrintHelmholtz[f,{"Rij","Fij","ROTf"}]

Example 3: Rössler attractor

dim=3; (* <-- dimension *)
xvec = Array[Subscript[x,#]&, dim]; (* define (Subscript[x, i]) *)
f={-(Subscript[x, 2])-(Subscript[x, 3]), (Subscript[x, 1])+"a" (Subscript[x, 2]), "b"+((Subscript[x, 1]) -"c")(Subscript[x, 3])};
PrintHelmholtz[f]

Example 4: Lorenz attractor

dim=3; (* <-- dimension *)
xvec = Array[Subscript[x,#]&, dim]; (* define (Subscript[x, i]) *)
f={"a"((Subscript[x, 2])-(Subscript[x, 1])), (Subscript[x, 1])("b"-(Subscript[x, 3]))-(Subscript[x, 2]),(Subscript[x, 1]) (Subscript[x, 2]) - "c" (Subscript[x, 3])};
PrintHelmholtz[f]

Example 5: Competitive Lotka–Volterra equations with n species

dim = 3; (* here: dim = n species *)
xvec = Array[Subscript[x,#]&, dim];
rvec = Array[Subscript["ρ",#]&, dim];
αij = Array[Subscript["α", #1, #2] &,{dim,dim}];
f = Table[rvec[[i]] (Subscript[x, i])(1-Sum[αij[[i,j]] (Subscript[x, j]),{j,dim}]),{i,dim}];
PrintHelmholtz[f,{"Fij"}]