Polynomial expansion of operator

I am new to Mathematica, I am trying to generate the polynomial function of a operator. So for example, the operator $L$ is $\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}$, and I want to generate the polynomial $\sum_{n=0}^{n=k} (L/2)^n$ and apply that polynomial operator to a function. I was trying to use Nest, any help??

Define

L = (1/2) (D[#, x] + D[#, y]) &

We see that L works as desired. For instance:

Simplify[Nest[L, f[x, y], 3]]
(* (Derivative[0, 3][f][x, y] + 3*Derivative[1, 2][f][x, y] + 3*Derivative[2, 1][f][x, y] +
Derivative[3, 0][f][x, y])/8 *)

And the Sum can be constructed in a similar manner. For instance:

Simplify[Sum[Nest[L, f[x, y], n], {n, 0, 3}]]
(* (8*f[x, y] + 4*Derivative[0, 1][f][x, y] + 2*Derivative[0, 2][f][x, y] +
Derivative[0, 3][f][x, y] + 4*Derivative[1, 0][f][x, y] +
4*Derivative[1, 1][f][x, y] + 3*Derivative[1, 2][f][x, y] + 2*Derivative[2, 0][f][x, y] +
3*Derivative[2, 1][f][x, y] + Derivative[3, 0][f][x, y])/8 *)

Update

As kindly pointed out by @b.gatessucks in the Comment below, computation of the final result can be simplified with NestList. (Thanks!)

Simplify[Total[NestList[L, f[x, y], 3]]]
• You can make it more efficient by using NestList as in Simplify[Total[NestList[L, f[x, y], 3]]]. – b.gates.you.know.what Jan 25 '15 at 8:40

Here is a function makeOperator that takes any polynomial together with a replacement rule that maps the desired variable onto the desired operator. It outputs the result as a new operator:

Clear[makeOperator];
makeOperator[poly_, Rule[x_, op_]] /; PolynomialQ[poly, x] :=
Module[{f},
Function[#1, #2] & @@ {f, Expand[poly]} /.
Power[x, n_: 1] :> Nest[op, f, n]]

I define operators using Function. Since that has attribute HoldAll, the necessary replacements in the polynomial have to be done outside the Function body, and are injected afterwards using Apply (@@). The pattern Power[x, n_: 1] detects powers of the variable (including first powers) and replaces them by Nest. Expand makes sure that the polynomial is in a canonical form before doing the replacements, in particular it eliminates parentheses like $x(x+c)$.

Here is a test with the operator L in the question, and a Hermite polynomial:

Clear[x, y];
L = Function[f, D[f, x] + D[f, y]];

hp = HermiteH[5, x]

(* ==> 120 x - 160 x^3 + 32 x^5 *)

hpOp = makeOperator[hp, x -> L];

hpOp[ψ[x, y]]

$$120 \left(\frac{\partial \psi }{\partial x}+\frac{\partial \psi }{\partial y}\right)\\ -160 \left(3 \frac{\partial ^3\psi }{\partial x^2\, \partial y}+3 \frac{\partial ^3\psi }{\partial x\, \partial y^2}+\frac{\partial ^3\psi }{\partial x^3}+\frac{\partial ^3\psi }{\partial y^3}\right)\\+32 \left(5 \frac{\partial ^5\psi }{\partial x^4\, \partial y}+10 \frac{\partial ^5\psi }{\partial x^3\, \partial y^2}+10 \frac{\partial ^5\psi }{\partial x^2\, \partial y^3}+5 \frac{\partial ^5\psi }{\partial x\, \partial y^4}+\frac{\partial ^5\psi }{\partial x^5}+\frac{\partial ^5\psi }{\partial y^5}\right)$$