7
$\begingroup$

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 $.

$\endgroup$
5
  • $\begingroup$ 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}] $\endgroup$
    – flinty
    Jul 22, 2020 at 0:29
  • $\begingroup$ 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]$. $\endgroup$
    – B. Erenu
    Jul 22, 2020 at 0:33
  • $\begingroup$ Please explain how you get that result - it's not obvious to me. $\endgroup$
    – flinty
    Jul 22, 2020 at 0:39
  • $\begingroup$ 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. $\endgroup$
    – B. Erenu
    Jul 22, 2020 at 0:42
  • 1
    $\begingroup$ 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 $\endgroup$
    – xzczd
    Jul 22, 2020 at 3:48

3 Answers 3

10
$\begingroup$
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] &

enter image description here

The coefficients are Stirling numbers of the first kind, StirlingS1

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

enter image description here

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 *)
$\endgroup$
2
  • 2
    $\begingroup$ +1 for StirlingS1 :) $\endgroup$ Jul 22, 2020 at 6:08
  • $\begingroup$ Thanks for the thorough response! $\endgroup$
    – B. Erenu
    Jul 23, 2020 at 18:10
8
$\begingroup$

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]
$\endgroup$
1
  • $\begingroup$ Thanks for your help! $\endgroup$
    – B. Erenu
    Jul 23, 2020 at 18:10
4
$\begingroup$
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]


Update based on your comment:

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]

$\endgroup$
1
  • $\begingroup$ 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]. $\endgroup$
    – B. Erenu
    Jul 22, 2020 at 1:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.