# Iterated differential operator

I would like to implement the following differential operator in Mathematica,

$$\prod_{j=1}^n \left({\mathrm d \over \mathrm d x}-j\right)$$

On a case by case basis, I can just expand it out and write things in terms of D[f, {x, j}] for $$j = 1, 2, \dots, n$$, but I would like to know how to define a function that works for general $$n$$.

• Is the product actually a product? You mention D[f,{x,j}] which takes the $j$th derivative of f w.r.t x, but (for univariate f) you could also mean op[f_, x_, n_] := Product[D[f[x], x] - j, {j, 1, n}] – flinty Jul 22 '20 at 0:29
• Sorry for poor wording. The product is supposed to define an $n$-th order differential operator. So for example, when applied to a test function $f[x]$ the $n=2$ case would give $f''[x]-3 f'[x]+2f[x]$. – B. Erenu Jul 22 '20 at 0:33
• Please explain how you get that result - it's not obvious to me. – flinty Jul 22 '20 at 0:39
• For $n=2$ the differential operator of interest is $(d/dx-1)(d/dx-2)= d^2/dx^2 - 3 d/dx + 2$. Then just apply this to the test function. – B. Erenu Jul 22 '20 at 0:42
• Related: mathematica.stackexchange.com/q/71643/1871 mathematica.stackexchange.com/q/125277/1871 Also, you may want to try this package: mathematica.stackexchange.com/a/162590/1871 – xzczd Jul 22 '20 at 3:48

Clear["Global*"]


Defining the operator recursively,

dOp[func_, x_Symbol, 1] := dOp[func, x, 1] = D[func, x] - func;

dOp[func_, x_Symbol, n_Integer?Positive] := dOp[func, x, n] =
D[dOp[func, x, n - 1], x] - n*dOp[func, x, n - 1];


For example,

dOp[f[x], x, 2] // Expand

(* 2 f[x] - 3 f'[x] + f''[x] *)


Looking at the first several,

Table[{n, dOp[f[x], x, n] // Expand}, {n, 1, 6}] //
Grid[#, Alignment -> Left, Dividers -> All] &


The coefficients are Stirling numbers of the first kind, StirlingS1

Table[StirlingS1[n, m], {n, 2, 7}, {m, 1, n}] // Grid


Consequently, the operator can alternatively be written as the sum

dOp2[func_, x_Symbol, n_Integer?Positive] :=
Sum[StirlingS1[n + 1, m + 1] D[func, {x, m}], {m, 0, n}]


Verifying the equivalence of the definitions,

And @@ Table[dOp[f[x], x, n] == dOp2[f[x], x, n] // Simplify, {n, 1, 15}]

(* True *)

• +1 for StirlingS1 :) – Αλέξανδρος Ζεγγ Jul 22 '20 at 6:08
• Thanks for the thorough response! – B. Erenu Jul 23 '20 at 18:10

One can use Composition (which should be the true meaning of the product of operators) with (pure) Function:

ClearAll[op, ff]
ff[j_] := f \[Function] D[f, x] - j f
op[n_] := Composition @@ Table[ff[j], {j, n}]

op[2][f[x]] // Simplify


2 f[x] - 3 f'[x] + f''[x]

Update

The version without Apply (@@) (actually Composition can be also replaced by RightComposition, because the order seems not to matter much in this case), by virtue of a syntax of Array able to replace heads:

Clear[op2]
op2[n_] := Array[ff, n, 1, Composition]

• Thanks for your help! – B. Erenu Jul 23 '20 at 18:10
op[f_, x_, n_] := Block[{d},
Total[MapIndexed[#1*If[#2[[1]] == 1, 1, D[f[x], {x, #2[[1]] - 1}]] &,
CoefficientList[Product[d - j, {j, 1, n}], d]]]]


op[g, y, 2] gives 2 - 3 g'[y] + g''[y]

op[f_, x_, n_] :=
Block[{d},
Total[MapIndexed[#1*D[f[x], {x, #2[[1]] - 1}] &,
CoefficientList[Product[d - j, {j, 1, n}], d]]]]


op[g, y, 2] gives 2*g[y] - 3*g'[y] + g''[y]`

• Excellent, thanks! I actually want the lowest term to be 2f[x] instead of 2, but this is easily fixed by changing the 1 after #2[[1]] == 1 to f[x]. – B. Erenu Jul 22 '20 at 1:05