Define operator algebra

Question

My objective is to work out commutators like this $$[f(x,y)\partial_x^2+g(x,y)\partial_x+h(x,y),a(x,y)\partial_x^2+b(x,y)\partial_x+c(x,y)]$$ or $$[f(x,y)\partial_x^2+g(x,y)\partial_x+h(x,y),a(x,y)\partial_y^2+b(x,y)\partial_y+c(x,y)]$$ by starting from simpler commutator identities, for example, using $$[\partial_x,f(x,z)]=f_x(x,z)$$ it is easy to derive $$[\partial^2_x,f(x,z)]=2f_x(x,z)\partial_x+f_{xx}(x,z)$$ and so on, until I get to the first two commutators that I want. But how can I do that with Mathematica?

My difficulty is that I don't know how to teach Mathematica to use identities like the third equation to derive the fourth and up to the first two.

Answer 1

Define the 2 operators whose commutator you wish to derive. I have ignored the dependence on $y$ to simplify the expressions without losing any essential content.

op1 = (f[x] Dt[#, {x, 2}] + g[x] Dt[#, x] + h[x] #) &;
op2 = (a[x] Dt[#, {x, 2}] + b[x] Dt[#, x] + c[x] #) &;

Apply the commutator to a test function $u(x)$, and then separately collect the $u(x), u^\prime (x), u^{\prime\prime}(x), u^{\prime\prime\prime}(x)$ terms.

op1[op2[u[x]]] - op2[op1[u[x]]] // Collect[#, {u[x], u'[x], u''[x], u'''[x]}] & 

(* the result is a bit too long to quote here *)

You can then read off the operator-valued version of the commutator by simply removing the test function $u(x)$ from this result.

Alternatively, a fully operator-valued derivation of the commutator could follow the lines I outline below, which was a quick bit of interactive operator manipulation that I concocted.

Define the 2 operators, using op[...] as a container to hold an "operator product", which allows us to control the order in which the various operators (and c-numbers) occur.

o1 = op[f, dx, dx] + op[g, dx] + op[h];
o2 = op[a, dx, dx] + op[b, dx] + op[c];

Form the commutator.

op[o1, o2] - op[o2, o1]

Distribute the op over Plus.

% //. op[u_ + v_, w_] :> op[u, w] + op[v, w]
% //. op[w_, u_ + v_] :> op[w, u] + op[w, v]

Flatten the nested op.

% /. op[op[u__], op[v__]] :> op[u, v]

Move the dx operators to the right.

% //. op[u___, dx, v_?(# =!= dx &), w___] :> op[u, v, dx, w] + op[u, d[v, x], w]

Move the c-numbers outside the op.

% //. op[u_?(# =!= dx &), v___] :> u op[v]

Tidy up the notation.

% /. d[d[u_, v_], v_] :> d[u, {v, 2}]
% /. op[] -> 1

Collect terms by derivative.

% // Collect[#, op[__]] &

This gives the same result as the quick derivation earlier.