An Integration–Annihilator method for analytical solutions of Partial Differential Equations

Oliver Richters, Erhard Glötzl
Supplementary material for: arXiv:???.???.

This worksheet implements the Integration–Annihilator method and applies it to the Poisson equation, the Polyharmonic equation and the Helmholtz equation, each with unweighted and weighted Laplacians. Also, we calculate a Helmholtz decomposition for various n-dimensional vector fields.

We first define some helper functions and then the DeriveQ function that implements Theorems 4.1–4.3 if λ and W are given by the reader. The Autocalculate function allows to call various HelperAuto functions that automatically determine λ and W for different equations.

Version v1. May 2025. Updates will be published at https://oliver-richters.de/annihilator and doi:10.5281/zenodo.15363375.

Integration–Annihilator method

Integrals and Incomplete Laplacians

(* a vector with only ones irrespective of length *)
OneVector /: Part[OneVector, _] := 1

(* 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],{xvec[[m]],0,xvec[[m]]}]];

(* Weighted Laplacian *)
WLaplacian[func_, w_]:=Sum[w[[i]] D[func, {xvec[[i]], 2}], {i,1,dim}];

(* Incomplete Laplacian: (Subscript[Δ, \m]) *)
LaplacianIncomplete[func_, m_]:=Laplacian[func, xvec] - D[func, {xvec[[m]], 2}];
WLaplacianIncomplete[func_, m_, w_]:=WLaplacian[func, w] - w[[m]] D[func, {xvec[[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]];
WLaplacianIncompleteP[func_, p_, m_, w_]:=Switch[p,0,func,_,WLaplacianIncomplete[WLaplacianIncompleteP[func, p-1, m, w], m, w]];

(* Laplacian to the power of p: (Δ^p) *)
LaplacianP[func_, p_]:=Simplify[Switch[p,0,func,_,Laplacian[LaplacianP[func, p-1], xvec]]];
WLaplacianP[func_, p_, w_]:=Simplify[Switch[p,0,func,_,WLaplacian[WLaplacianP[func, p-1, w], w]]];

(* Generalized Helmholtz operator (Δ^j + ν)^k *)
Helmholtz[func_, ν_, j_, k_, w_] := Simplify[Switch[k,0,func,_,Helmholtz[WLaplacianP[func,j,w] + ν*func, ν, j, k-1, w]]]

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

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])}];

Integration–Annihilator method: Theorems 4.1–4.3

(* solution of Theorems 4.1 and 4.2 *)
DeriveQ[q_,opD_,opA_,λ_,W_,k_:1,u_:1] := Module[{opB},
  opB[g_]:=opD[g]-opA[g];
  On[Assert];
  Assert[Not[VectorQ[q]], "q must be a scalar function"];
  Assert[1 ≤ k, "k must be at least 1"];
  Assert[FullSimplify[Nest[opA, W, λ+k] - q] == 0,
    StringJoin["Differentiating W=,",W," with respect to A=",opA," λ+k-times with λ=,",λ, " must yield q=",q]];
  Assert[FullSimplify[Nest[opB, W, λ+1] - (-1)^λ (u-1)q] == 0, "Terminal condition not fulfilled"];
  Assert[k == 1 || u == 1, "u=1 is required as k>1 does not support Theorem 4.2."];
  Assert[u≠0 || u==1-"a"^2/"w"^2 (*added to avoid assertion errors*), "u must be non-zero"];
  Return[Simplify[Sum[1/u Binomial[k+p-1,p] (-1)^p Nest[opA, Nest[opB, W, p], λ-p], {p,0,λ}]]];
]; (* end of function SolutionD *)

(*AutocalculateSum expands input into list and applies solver elementwise *)
(*thesolver are defined as HelperAuto functions below for each equation specifically *)
AutocalculateSum[f_,k_,thesolver_,args__] := Module[{M,F,Poly,PrintInfo,pol,v,xiInv,opD},
PrintInfo=False; (* switch verbose printing on / off here *)
F = 0;
Poly = Expand[f];
(* Make sure Poly is a proper list, even if f has only one summand *)
Poly = If[Head[Poly]===Plus, List@@Poly, List[Poly]]; 
If[PrintInfo, Print[Poly]];
opD = thesolver["returnD", k, args];
Do[
  M = Poly[[pol]];
  If[M === 0, Continue[]];
  (* catch case where D Q = v q with v a real number / Theorem 4.3. *)
  v = FullSimplify[opD[M]/M];
  xiInv = Intersection[xvec, getUniqueVariables[v]];
  If[! v === 0 && Length[xiInv] == 0,
   F = F + M/v^k,
   F = F + thesolver[M,k,args];
  ];
, {pol,Length[Poly]}
];
Return[F]
]; (* end of function AutocalculateSum *)

(* Autocalculate determines whether q is a vector or not and either maps or not *)
Autocalculate[q_,args__] := Module[{},
If[VectorQ[q],
  Return[Map[Autocalculate[#,args] &, q, {Length[Dimensions[q]]}]],
  Return[AutocalculateSum[q,args]]
]
]; (* end of function Autocalculate *)

Poisson and Polyharmonic Equation / (Δ^k)Q = q

Poisson and polyharmonic equation with weighted Laplacian, see Corollary 5.1–5.3

The function call is SolvePolyharmonic[q,k,w] with w an n-dimensional weight vector for the Laplacian.
You can also drop w and call the function with two arguments only, defaulting to an unweighted Laplacian.
If called with one single argument, k defaults to 1, which is equivalent to SolvePoisson[q].

Helper functions

(* These HelperManual functions returns a solution if m, λ, W (and u) are given *)

(* calculate solution for Poisson equation given m, λ, u, W *)
HelperManualPoisson[q_,m_,λ_,u_,W_]:=Module[{g,opA,opD},
  On[Assert];
  Assert[1 ≤ m ≤ dim, "m must be between 1 and dim"];
  opA[g_]:=D[g,{xvec[[m]],2}];
  opD[g_]:=Laplacian[g,xvec];
  Return[DeriveQ[q,opD,opA,λ,W,1,u]];
]; (* end of function SinglePoisson *)

(* calculate solution for polyharmonic equation with weighted Laplacian given m, λ, W, k and weight w *)
HelperManualPolyharmonic[q_,m_,λ_,W_,k_:1,w_:OneVector]:=Module[{g,opA,opD},
  On[Assert];
  Assert[1 ≤ m ≤ dim, "m must be between 1 and dim"];
  Assert[Length[w]==dim || w == OneVector, "w must have length dim"];
  opA[g_]:=w[[m]] D[g,{xvec[[m]],2}];
  opD[g_]:=WLaplacian[g,w];
  Return[DeriveQ[q,opD,opA,λ,W,k]];
]; (* end of function DerivePolyharmonic *)

(* This HelperAuto function tries to automatically find a solution for a simple expression (no sums!) === *)
(* Don't call this function directly, use the Solve functions below *)

(* Polyharmonic equation with weighted Laplacian, see Corollary 5.1–5.3 *)
HelperAutoPolyharmonic[S_,k_,w_] := Module[{Q1,Q2,m,λ,W,i},
  opD[g_]:=WLaplacian[g,w];
  If [S==="returnD", Return[opD]];
  PrintInfo=False; (* switch verbose printing on / off here *)
  isPolyWithoutxi = Table[PolynomialQ[S, Intersection[xvec, DeleteCases[getUniqueVariables[S], xvec[[i]]]]], {i,dim}];
  Q1 = AnyTrue[isPolyWithoutxi, TrueQ];
  If [Q1,
    Q2 = PolynomialQ[S,xvec];
    If[Q2,
      m = First[Position[Exponent[S,xvec],Max[Exponent[S,xvec]]]][[1]];
      , (* else *)
      m = First[Position[isPolyWithoutxi, True]][[1]];
    ];
    λ = Ceiling[(Total[Exponent[S,xvec]] - Exponent[S,xvec[[m]]] - 1)/2];
    W = IntegrateXmP[S, 2λ+2k, m] w[[m]]^(-λ-k);
    If[PrintInfo, Print[S," in q is partly monomial, use m=",m,", λ=",λ,", W=",W]];
    Return[HelperManualPolyharmonic[S,m,λ,W,k,w]]
  ];
  Print["No decomposition found for expression ", S, " in q."];
  Abort[];
];
(* end of function HelperPolyharmonic *)

(* These functions should be called *)
SolvePolyharmonic[q_,k_,w_] := Autocalculate[q, k, HelperAutoPolyharmonic, w];
SolvePolyharmonic[q_,k_:1] := SolvePolyharmonic[q,k,OneVector];
SolvePoisson[q_] := SolvePolyharmonic[q]

PrintLaplacianP[q_, Q_, k_:1, w_:OneVector] := Module[{qnew},
  Print["q = ", Simplify[q]];
  Print["Q = ", Q];
  qnew = FullSimplify[WLaplacianP[Q, k, w]];
  Print["Δ^", k, " Q = ", qnew];
  Print["Is Δ^", k, " Q = q? ", TrueQ[FullSimplify[q - qnew == 0]]];
];

Poisson Equation example: monomial

dim=4;
xvec = Array[Subscript[x,#]&, dim];
q = (Subscript[x, 1])^7 (Subscript[x, 2])^7;
λ = 9;
Q = HelperManualPoisson[q, m=1, λ, u=1, W=IntegrateXmP[q,2λ+2,m=1]];
PrintLaplacianP[q, Q, 1];

Poisson Equation example: sin, cos, exp and monomial

q = Cos[5(Subscript[x, 3])] Sin[3(Subscript[x, 1])] + Exp[2(Subscript[x, 1])] (Subscript[x, 1])^3 (Subscript[x, 2])^4;
Q = SolvePoisson[q];
PrintLaplacianP[q, Q];

Polyharmonic Equation example: product of sines

q = Sin[(Subscript[x, 1])]Sin[(Subscript[x, 2])];
k = 4;
Q = SolvePolyharmonic[q, k];
PrintLaplacianP[q, Q, k];

Weighted Polyharmonic Equation example: sin, cos, exp and monomial

q = Sin[3(Subscript[x, 1])] Cos[5(Subscript[x, 3])] + Exp[2(Subscript[x, 1])] (Subscript[x, 1])^3 (Subscript[x, 2])^4;
wnew = Join[{a}, ConstantArray[1, Length[xvec]-1]];
k=2;
Q = SolvePolyharmonic[q, k, wnew];
PrintLaplacianP[q, Q, k, wnew];

Wave equation

dim = 4;
xvec = Join[{t}, Array[Subscript[x,#]&, dim  - 1]];
Print["xvec = ", xvec];
wvec = Join[{1}, ConstantArray[-c^2, dim  - 1]];
Print["wvec = ", wvec];
q = t Sin[t] (Subscript[x, 1])^2 (Subscript[x, 2]);
k=1;
Print["PDE: ", WLaplacianP[A[xvec],k,wvec]==q];
Q=SolvePolyharmonic[q, k, wvec];
PrintLaplacianP[q, Q, k, wvec];
Qmanual = (-2 Cos[t]-t Sin[t])(Subscript[x, 1])^2 (Subscript[x, 2]) + 2(Subscript[x, 2])c^2(4Cos[t] + t Sin[t]);
PrintLaplacianP[q, Qmanual,k,wvec];

Helmholtz Equation ((Δ^j)-ν)kQ = q

Helmholtz equation with weighted Laplacian, see Corollary 5.4

The function call is SolveHelmholtz[q,ν,j,k,w] with w an n-dimensional weight vector for the Laplacian.
You can also drop w and call the function with four arguments only, defaulting to an unweighted Laplacian.
If called with even less arguments, j and k default to 1.

Helper functions

(* This HelperAuto function tries to automatically find a solution for a simple expression (no sums!) === *)
(* Don't call this function directly, use the Solve functions below *)
HelperAutoHelmholtz[S_,k_,ν_,j_,w_] := Module[{λ,g,opA,opD,W},
 On[Assert];
 opA[g_]:=ν*g;
 opD[g_]:=Helmholtz[g, ν, j, 1, w];
 Assert[Length[w]==dim || w == OneVector, "w must have length dim"];
 If [S==="returnD", Return[opD]];
 If[PolynomialQ[S,xvec],
  λ = Ceiling[(Total[Exponent[S,xvec]]-1)/(2j)];
  W = S*ν^(-λ-k);
  Return[DeriveQ[S,opD,opA,λ,W,k]];
 ];
 Print["No decomposition found for expression ", S, " in q."];
 Abort[];
];

(* These functions should be called *)
SolveHelmholtz[q_,ν_,j_,k_,w_] := Autocalculate[q, k, HelperAutoHelmholtz, ν, j, w];
SolveHelmholtz[q_,ν_,j_:1,k_:1] := SolveHelmholtz[q, ν, j, k, OneVector];
SolveHelmholtz[q_,ν_] := SolveHelmholtz[q, ν, 1, 1];

Example with unweighted Laplacian and k=j=1.

dim=4;
xvec = Array[Subscript[x,#]&, dim];

q =(Subscript[x, 1])^4 (Subscript[x, 2])^3;
Print["PDE: ", Helmholtz[A[xvec], v, 1, 1, OneVector] ==0];
Q = Expand[SolveHelmholtz[q,v]];
Print["q = ", q]
Print["Q = ", Q]
Print["(Δ - v)Q = ", Simplify[Laplacian[Q, xvec] + v*Q]];

Example with weighted Laplacian and k=5 and j=3.

wnew = ConstantArray[1, dim];
wnew[[1]] = 2;
j = 3;
k = 5;
ν = 3;
Print["Weighted Helmholtz Equation with ν=", ν, ", k=", k, ", j=", j, ":"];
Q = SolveHelmholtz[q,ν,j,k,wnew];
Print["q = ", q]
Print["Q = ", Q]
Print["((((Δ^j) - ν))^k)Q = ", Simplify[Helmholtz[Q,ν,j,k,wnew]]];

Helmholtz Decomposition

For a detailed introduction, see: Erhard Glötzl, Oliver Richters: Helmholtz decomposition and potential functions for n-dimensional analytic vector fields. In: Journal of Mathematical Analysis and Applications 525(2), 127138, September 2023,  doi:10.1016/j.jmaa.2023.127138 (arXiv:2102.09556).

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

Definitions

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}]];

(* === 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 *)

(* === AutocalculateFij[f] returns Fij for a given vector field f, if possible === *)
AutocalculateFij[f_] := Module[{M,Fij,Poly,PrintInfo,pol,k,i,j},
Return[Grad[SolvePoisson[f], xvec]]
];

PrintHelmholtz Helper function

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]];
)

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: f = [((Subscript[x, k])^2) (Subscript[x, m]), 0, 0]

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 s.t.: ((∂_(Subscript[x, m]))^(2λ)) W = fk, (((-1))^λ) ((Subscript[Δ, \m])^λ) W == (1-u) fk *)
W=IntegrateXmP[fk,2λ,m]; (* you may specify different functions here *)
Fij = CalculateFij[f,k,m,λ,u,W];
PrintHelmholtz[f,Fij]

Example with 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 with 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: f = [((Subscript[x, 1])^2) ((Subscript[x, 2])^20), 0]

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: 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) ≠ 0 *)
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: 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: 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: 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"}]