In the course of a long investigation I had been doing on the Adomian decomposition (and its implementation in Mathematica), I encountered an interesting paper by Duan showing how to recursively generate the underlying polynomials. That long note will be posted some other day, but for now, I would like to demonstrate how the recursive method can be used to compute the partial Bell polynomials.
The partial Bell polynomials
Y
n,k
(
x
1
,…,
x
n-k+1
)
are already built-in as BellY, but the following method can surely be adapted to other polynomial sets with similar structure. In particular, consider the following defining relation for the partial Bell polynomials:
In[]:=
With{n=7,k=4},​​BellYn,k,Array[x.,n],Sumn!ApplyTimes,
m
Array[x.,n]
Range[n]!
1
m!
,{m,FrobeniusSolve[Range[n],n]⋂FrobeniusSolve[ConstantArray[1,n],k]}
Out[]=
105x.[1]
3
x.[2]
+210
2
x.[1]
x.[2]x.[3]+35
3
x.[1]
x.[4],105x.[1]
3
x.[2]
+210
2
x.[1]
x.[2]x.[3]+35
3
x.[1]
x.[4]
The sticky part here is in the generation of the summation indices through two calls to FrobeniusSolve; that is, the nonnegative indices
j
1
,…,
j
n-k+1
that satisfy the conditions
j
1
+…+
j
n-k+1
=k
and
j
1
+2
j
2
+…+(n-k+1)
j
n-k+1
=n
. Even for modestly-sized
n
, the computation can already take quite a while:
In[]:=
AbsoluteTiming[res1=With[{n=13},Table[FrobeniusSolve[Range[n],n]⋂FrobeniusSolve[ConstantArray[1,n],k],{k,n}]];]
Out[]=
{23.833,Null}
Theorem 1 in Duan’s paper shows the recursive approach for generating these summation indices; as usual, we use memoization in its implementation:
In[]:=
idx[n_Integer?Positive,1]:={UnitVector[n,n]}
In[]:=
idx[n_Integer?Positive,k_Integer?Positive]:=idx[n,k]=​​Module[{res},​​res=Table[PadRight[v,n]+UnitVector[n,1],{v,idx[n-1,k-1]}];​​If[k<=Quotient[n,2],​​res=Union[res,Table[ArrayPad[v,{1,k-1}],{v,idx[n-k,k]}]]];​​res]
Compared to the FrobeniusSolve approach, generating a whole pile of these indices recursively is quite fast:
In[]:=
AbsoluteTiming[res2=With[{n=13},Table[idx[n,k],{k,n}]];]
Out[]=
{0.0044075,Null}
In[]:=
res1===res2
Out[]=
True
Even larger values of
n
are nearly unmanageable with the Frobenius approach:
In[]:=
TimeConstrained[With[{n=19,k=11},FrobeniusSolve[Range[n],n]⋂FrobeniusSolve[ConstantArray[1,n],k]],10]
Out[]=
$Aborted
The recursive approach still works like a charm:
In[]:=
With{n=19,k=11},​​BellYn,k,Array[x.,n],Sumn!ApplyTimes,
m
Array[x.,n]
Range[n]!
1
m!
,{m,idx[n,k]}
Out[]=
1964187225
3
x.[1]
8
x.[2]
+18332414100
4
x.[1]
6
x.[2]
x.[3]+36664828200
5
x.[1]
4
x.[2]
2
x.[3]
+16295479200
6
x.[1]
2
x.[2]
3
x.[3]
+775975200
7
x.[1]
4
x.[3]
+10999448460
5
x.[1]
5
x.[2]
x.[4]+24443218800
6
x.[1]
3
x.[2]
x.[3]x.[4]+6983776800
7
x.[1]
x.[2]
2
x.[3]
x.[4]+2618916300
7
x.[1]
2
x.[2]
2
x.[4]
+436486050
8
x.[1]
x.[3]
2
x.[4]
+3666482820
6
x.[1]
4
x.[2]
x.[5]+4190266080
7
x.[1]
2
x.[2]
x.[3]x.[5]+349188840
8
x.[1]
2
x.[3]
x.[5]+523783260
8
x.[1]
x.[2]x.[4]x.[5]+11639628
9
x.[1]
2
x.[5]
+698377680
7
x.[1]
3
x.[2]
x.[6]+349188840
8
x.[1]
x.[2]x.[3]x.[6]+19399380
9
x.[1]
x.[4]x.[6]+74826180
8
x.[1]
2
x.[2]
x.[7]+11085360
9
x.[1]
x.[3]x.[7]+4157010
9
x.[1]
x.[2]x.[8]+92378
10
x.[1]
x.[9],1964187225
3
x.[1]
8
x.[2]
+18332414100
4
x.[1]
6
x.[2]
x.[3]+36664828200
5
x.[1]
4
x.[2]
2
x.[3]
+16295479200
6
x.[1]
2
x.[2]
3
x.[3]
+775975200
7
x.[1]
4
x.[3]
+10999448460
5
x.[1]
5
x.[2]
x.[4]+24443218800
6
x.[1]
3
x.[2]
x.[3]x.[4]+6983776800
7
x.[1]
x.[2]
2
x.[3]
x.[4]+2618916300
7
x.[1]
2
x.[2]
2
x.[4]
+436486050
8
x.[1]
x.[3]
2
x.[4]
+3666482820
6
x.[1]
4
x.[2]
x.[5]+4190266080
7
x.[1]
2
x.[2]
x.[3]x.[5]+349188840
8
x.[1]
2
x.[3]
x.[5]+523783260
8
x.[1]
x.[2]x.[4]x.[5]+11639628
9
x.[1]
2
x.[5]
+698377680
7
x.[1]
3
x.[2]
x.[6]+349188840
8
x.[1]
x.[2]x.[3]x.[6]+19399380
9
x.[1]
x.[4]x.[6]+74826180
8
x.[1]
2
x.[2]
x.[7]+11085360
9
x.[1]
x.[3]x.[7]+4157010
9
x.[1]
x.[2]x.[8]+92378
10
x.[1]
x.[9]

CITE THIS NOTEBOOK

A recursive method for partial Bell polynomials​
by J. M.​
Wolfram Community, STAFF PICKS, April 25, 2025
​https://community.wolfram.com/groups/-/m/t/3450233