Exponential of a Differential Operator

Question

In Mathematica, is it possible to exponentiate a differential operator such that the operator will act on a function, $f(x,p)$? Specifically, I wondering if I can get Mathematica to do this:

$\exp(c\frac{d}{dp})f(x,p)=(1+c\frac{d}{dp}+\frac{c^2}{2}\frac{d^2}{dp^2}+\cdots)f(x,p)=(f(x,p)+c\frac{df}{dp}+\frac{c^2}{2}\frac{d^2f}{dp^2}+\cdots)$

where $c\in\Re$.

Thanks.

Answer 1

There are two possible approaches:

(1)

To map an operator exponential onto a MatrixExp, one would have to assume that there exists an orthogonal basis of the function space in which the target functions f[x] can be expanded. Then construct a matrix representing the operator in that basis. For example, the momentum operator $p$ is diagonal in the Fourier basis, and therefore we can do the exponential very easily in that basis. But we may not always be able to determine the expansion coefficients in a chosen basis in closed form. As a result, you end up with matrices corresponding to a truncated expansion of the function f[x]. Exact results can be obtained in this approach if your function space is finite-dimensional and the operator doesn't throw you out of that space.

(2)

Instead of truncating the expansion of f[x] in some basis, you can truncate the series for the exponential at some order of the expansion. Truncation in either approach genrally can't be avoided, but this approach is very common since the exponential series is absolutely convergent over the complex numbers.

However, the usual form of the series expansion for the exponential requires you to calculate factorials at every order. This is inefficient from a computational point of view.

Instead, the exponential of an operator can be conveniently (and efficiently) defined by using a recursive formula for the exponential. This is implemented here:

Clear[f, operatorExp];
operatorExp[dop_, n_] := Function[f, Fold[f + dop[#1]/#2 &, f, Reverse@Range[n]]]

This one-liner defines the exponential of dop up to a number n of terms. The result is itself an operator, taking as its argument a function of the form f[x]. Here I'll illustrate what this does to a general operator dop:

Clear[dop,x];
operatorExp[dop, 4][f[x]]

$\text{dop}\left(\frac{1}{2} \text{dop}\left(\frac{1}{3} \text{dop}\left(\frac{1}{4} \text{dop}(f(x))+f(x)\right)+f(x)\right)+f(x)\right)+f(x)$

Now try it for a well-known example (the Taylor series), and then for an operator with non-commuting terms. To do this, I use Function to define the operators dOp and dOp2. They assume that the independent variable of the differential operation is x.

dOp = Function[f, D[f, x]];
operatorExp[dOp, 3][f[x]]

$f(x) + f'(x)+\frac{1}{2} \left(\frac{1}{3} f^{(3)}(x)+f''(x)\right)$

dOp2 = Function[f, x + D[f, x]];
operatorExp[dOp2, 3][f[x]]

$x + f(x) + f'(x)+\frac{1}{2} \left(\frac{1}{3} f^{(3)}(x)+f''(x)+1\right)$

Applied to the specific problem in the question

Here I use the momentum as the differentiation variable in the operator:

Clear[c,p];
dOp3 = Function[f, c D[f, p]];
Expand[operatorExp[dOp3, 3][f[x, p]]]

$\frac{1}{6} c^3 f^{(0,3)}(x,p)+\frac{1}{2} c^2 f^{(0,2)}(x,p)+c f^{(0,1)}(x,p)+f(x,p)$

By using Expand, I turned the recursively factored form back into something that looks like the usual Taylor expansion.

Answer 2

may be something like this:

Plus @@ (c^#/#! D[f[x, p], {p, #}] & /@ Range[0, 5])

or

Sum[c^i/i! D[f[x, p], {p, i}], {i, 0, 5}]

Answer 3

Here is a trick for making a series expansion of a function of a single operator — i.e. the series expansion has no non-commutation issues — and then correctly manipulating powers of the operator, where the operator may have non-commuting internal parts.

Expand the exponential Exp[op] of an operator op — 3 terms are sufficient for illustrative purposes.

Series[Exp[op], {op, 0, 3}] // Normal

(* 1 + op + op^2/2 + op^3/6 *)

Convert the leading 1 into an identity operator id.

% // ReplacePart[#, 1 -> id] &

(* id + op + op^2/2 + op^3/6 *)

Operate on a state s.

(%)[s] // Through

(* id[s] + op[s] + ((op^2)/2)[s] + ((op^3)/6)[s] *)

Convert the powers of operators into the equivalent composition of operators. This trick sets things up so that, when op is expanded out in terms of its internal parts, any non-commuting parts in op are handled correctly.

% /. ((u_: 1) op^n_)[v_] :> u (Composition @@ Table[op, {n}])[v]

(* id[s] + op[s] + 1/2 op[op[s]] + 1/6 op[op[op[s]]] *)

Now use coordinate representation s[x] (i.e. a wave function) for the state s, and (for example) use (# + I D[#, x]) &) for the operator op .

% /. s -> s[x] /. {id -> (# &), op -> ((# + I D[#, x]) &)} // Expand

(* (8*s[x])/3 + (5/2)*I*Derivative[1][s][x] - Derivative[2][s][x] - (1/6)*I*Derivative[3][s][x] *)

This approach is very general — i.e. the (# + D[#, x])& is merely an example of the sort of operator you can use.

Answer 4

One way would be to get the series by treating the operator as a variable, then replace the variable with operator. Let say, you want the series to order 9.

n = 9;
op = Series[Exp[c dp], {dp, 0, n}] // Normal;
op /. {dp^m_ -> D[#, {p, m}]} &@f[x, p]

Here I first treat the operator as variable dp first and then replace it with the operator once I have my series. And you have your answer. You can use it for any function of operator (I mean in general, not those tricky function with singularity or even worse complexity!)

For any arbitrary form of operator

I can still use the same method, but the power will be taken care by Nest which will allow the repeated application of the operator.

opd[x_]=(D[#,x]#+x #)&; (*The operator*)
func[op_]:=Exp[c op]; (*The func. of operator*)
f[x_,p_]=y[x ,p]; (*Target function*)

n=1; (*order of series expansion*)
opr[op_]=op Series[func[op],{op,0,n}]//Normal; 
out[g_,x_,p_]=(opr[op])/.op^m_.->Nest[opd[x],g[x,p],m-1]//Quiet;

out[f,x,p]//Expand

The output is $y[x,p]+c x y[x,p]+c y[x,p] y^{(1,0)}[x,p]$.

Now you can use any form or function of operator.

Regarding this I like to mention one thing. While making the series opr[op] I multiply it with op so that the power series starts from op instead of 1. That is why I use the Nest m-1 times for op^m. Because 1 is not interpreted as op^0 and in final result it shows 1 instead of y[x,p]. I don't know if there is any better way to handle it. I post this question little ago. I would be glad to hear comments about that which may also improve this answer.

Better Approach

opd[x_]=(D[#,x]#+x #)&;(*The operator*)
func[op_]:=Exp[c op];(*The func.of operator*)
f[x_,p_]=y[x,p];(*Target function*)
n=1;(*order of series expansion*)
opr[op_]= Series[func[op],op,0,n}]//Normal;
out[g_,x_,p_]=Sum[Coefficient[opr[op],op,m]Nest[opd[x],g[x,p],m],{m,0,n}];

out[f,x,p]//Expand