A function which maps $(ax + by)(cx + dy) \mapsto (a \partial_x + b \partial_y) (c \partial_x + d \partial_y)$

I am new to Mathematica and, as the title says, looking for a way of mapping (for example) the polynomial $$(ax + by)(cx + dy) \mapsto \left(a \frac{\partial}{\partial x} + b \frac{\partial}{\partial y}\right) \left(c \frac{\partial}{\partial x} + d \frac{\partial}{\partial y}\right) \\ = ac \frac{\partial^2}{\partial_x^2} + (ac + bd) \frac{\partial}{\partial_x}\frac{\partial}{\partial_y} + bd \frac{\partial^2}{\partial_y^2},$$ where the latter operator operates on some predefined function. The $\{a,b,c,d\}$ are constants.

Of course this is simple to do by hand in the case shown above, but I will need to be doing this for rather large polynomials where the manual way is not possible. A systematic approach of some kind would be appreciated.

Thanks to anyone who can help me.

• Closely related, but just for a single operator: Polynomial expansion of operator. You can use makeOperator from my answer here, too, because your example doesn't contain any operators that don't commute. If that were the case, one would have to do some additional work. – Jens May 16 '17 at 22:21
• The parameters {a,b,c,d} are intended as constants (that is, independent of {x,y})? – Daniel Lichtblau May 16 '17 at 23:16
• @DanielLichtblau Yes. – oscarafone May 16 '17 at 23:17

Maybe something like the following:

L /: L[x___]L[y___] := Sort @ L[x,y]
L /: Power[L[x__], n_Integer?Positive] := L @@ Flatten @ ConstantArray[{x}, n]

toOperator[expr_, var_] := With[
{op = Collect[expr /. v:var->L[v], _L] /. L[x__] -> Inactive[D][#, x]},
Activate @ Function @ op
]

For the OP example:

toOperator[(a x + b y)(c x + d y), x|y] //TeXForm

$(a d+b c) \frac{\partial ^2\text{$\#$1}}{\partial x\, \partial y}+a c \frac{\partial ^2\text{$\#$1}}{\partial x\, \partial x}+b d \frac{\partial ^2\text{$\#$1}}{\partial y\, \partial y}\&$

• Excellent! I have spent some time trying to understand this solution (I'm a novice wrt Mathematica) and I cannot quite get it. But I am extremely happy that you've found a solution. Thanks so much! – oscarafone May 30 '17 at 22:21

I have one solution.

ClearAll[PolyMap];
PolyMap[p1_ + p2_, vars_] /; (PolynomialQ[p1, vars] && PolynomialQ[p2, vars]) :=  PolyMap[p1, vars] + PolyMap[p2, vars];
PolyMap[a_, vars_] /; VectorQ[vars, ! MemberQ[a, #, Infinity] &] := a*# &;
PolyMap[a_*x_, vars_] /; (VectorQ[vars, ! MemberQ[a, #, Infinity] &] && MemberQ[vars, x]) := a D[#, x] &;
PolyMap[a_*x_^n_, vars_] /; (IntegerQ@n && n > 0 && MemberQ[vars, x] && VectorQ[vars, ! MemberQ[a, #, Infinity] &]) := a*D[#, {x, n}] &;
PolyMap[a_*x_,vars_] /; (PolynomialQ[a x, vars] && MemberQ[vars, x]) := D[#, x] &@*PolyMap[a, vars];
PolyMap[a_*x_^n_, vars_] /; (PolynomialQ[a x^n, vars] && IntegerQ@n && n > 0&& MemberQ[vars, x]) := D[#, {x, n}] &@*PolyMap[a, vars];

Here is an Example:

PolyMap[a x^2 + b x y^4 + c + c/d x y, {x, y}] It looks a little of ugly, but it works:

PolyMap[a x^2 + b x y^4 + c + c/d x y, {x, y}][x^3 y^3 z^3] // Through
(*(9 c x^2 y^2 z^3)/d + 6 a x y^3 z^3 + c x^3 y^3 z^3*)