# Computing with matrix differential operators

I am looking the achieve a working set of Mathematica functions that would allow me to take products, commutators, etc. of (2x2) matrix differential operators and apply them to vector-valued functions of the form {f[x], g[x]}.

Here is what I have now and the problem I am encountering. The function

CircleTimes[a_,b_]:=Inner[Composition[#1,#2]&,a,b]


does the right job of computing the matrix operator resulting from the product of two matrix differential operators a and b. With the function

CenterDot[a_,v_]:=Inner[Through[#1[#2]&, Plus],a, v]


one can compute the action of the matrix differential operator a on the vector-valued function v, and it seems to be working fine.

Problems arise when I want to compute the products of more than 2 operators or their commutator on vector-valued functions.

For example, if one defines

d=CircleTimes[a,CircleTimes[b,c]]


which is the matrix operator corresponding to a x b x c, then

v=CenterDot[d,{f[x],g[x]}]


does not produce anything useful. The same happens if we define for example the anticommutator

e=Plus[CircleTimes[a,b], CircleTimes[b,a]]


then

w=CenterDot[e,{f[x],g[x]}]


does not apply properly. It seems that Mathematica does not complete the process of applying the various operators to the functions f[x] and g[x].

You can try with the following operators to see the problem

I1:={{0 &, (x # + D[#,x] &)}, {0 &, - D[#,x] &}
I2:= {{0 &, x^2# &}, {0 &, D[#,x] &}


We see that with

u={f[x],g[x]}


the following are perfectly computed

CenterDot[I1,u]
CenterDot[I2,u]
CenterDot[CircleTimes[I1,I2],u]


but trying to evaluate thing like the anticommutator in the following way

I3:=Plus[CircleTimes[I1,I2], CircleTimes[I2,I1]]

CenterDot[I3,u]


one finds that Mathematica does not complete the evaluation and we are left with stuff.

My guess is that I should be telling Mathematica to apply the "Through" function as long as there are things to evaluate, but I am not sure how this can be done.

Hope this is clear :)

Many thanks !

I add the following for the sake of information; see the first answer for a better solution (I think).

As a way to have Mathematica evaluate everything correctly, I came up with the following modification of the CenterDot[a_,b_] function. It reads

CenterDot[a_,b_]:=FixedPoint[Through[#,Plus]&/@# &,#]&/@Inner[#1[#2]&,a,b]


This arranges for Mathematica to perform all evaluations, distributing over the sums. It is also suitable for combination with the CirclePlus[a_,b_] defined in the first answer.

Here is a modification of your code that works. The problem was that you weren't combining the Compositions for each matrix entry correctly in CircleTimes: using Plus with Composition doesn't lead to a function which can be applied subsequently to another expression. You have to make each sum of compositions into a function itself.

First I'll copy the definitions that I didn't modify:

Clear[f, g, x];

I1 := {{0 &, (x # + D[#, x] &)}, {0 &, -D[#, x] &} }

I2 := {{0 &, x^2 # &}, {0 &, D[#, x] &}}

u = {f[x], g[x]};

CenterDot[a_, v_] := Inner[Through[#1[#2] &, Plus], a, v]


Now my modifications:

Clear[CircleTimes];
CircleTimes[a_, b_] :=
Inner[Composition[#1, #2] &, a, b,
Function[f, Total@Through[{##}[f]]] &]

CirclePlus[a_, b_] :=
MapThread[Function[f, Total@Through[{##}[f]]] &, {a, b}, 2]


The modification in CircleTimes is that I specified the replacement of the Plus operation in Inner as Function[f, Total@Through[{##}[f]]] & which is now itself a function that adds the individual operators after applying them to whatever argument f you provide.

Next, the same issue with Plus also arises in I3. So I defined an addition of operators CirclePlus in which the same trick is used as above, by combining the matrices of operators element-wise using MapThread.

Now the tests:

CenterDot[I1, u]

(* ==> {x g[x] + Derivative[g][x], -Derivative[g][x]} *)

CenterDot[I2, u]

(* ==> {x^2 g[x], Derivative[g][x]} *)

CenterDot[CircleTimes[I1, I2], u]

(*
==> {x Derivative[g][x] + (g^\[Prime]\[Prime])[x], -(
g^\[Prime]\[Prime])[x]}
*)

I3 := CirclePlus[CircleTimes[I1, I2], CircleTimes[I2, I1]]

CenterDot[I3, u]

(*
==> {x Derivative[g][x] -
x^2 Derivative[g][x] + (g^\[Prime]\[Prime])[x], -2 (
g^\[Prime]\[Prime])[x]}
*)

Id := {{1 &, 0 &}, {0 &, 1 &}}

CenterDot[CirclePlus[Id, Id], u]

(* ==> {2, 2} *)


As you can see, the tests produce outcomes that look correct. Edit: the last example is from the comment, and I changed the code to accommodate such constant operators. Note that this is not the same as the identity which multiplies the test function by 1. My understanding of the comment was that the desired operator was meant to replace the test function by 1.

• Thank you very much ! This explains a lot. The only problem remaining is that something strange happens when the operators are ordinary matrices. For example if we want to add the Identity to the Identity with CirclePlus[a_,b_] as defined above, one gets thins like 2[f[x]] like "2" is an operator waiting to be evaluated. (Even though the identity was defined as a true operator) Dec 17, 2014 at 0:16
• Ah yes, that's because Plus inside Thread is not held unevaluated when its arguments have values that are known at the time of evaluation. To make it works for such cases it's probably best to avoid Plus - I'll edit the answer to deal with the example you gave.
– Jens
Dec 17, 2014 at 2:46

This is how I would define your functions, similar but I think simpler than @Jens answer:

ClearAll[CirclePlus, CircleTimes, CenterDot]

SetAttributes[CirclePlus, {Flat, Listable}]

CirclePlus[a__][x_] := Activate @* Through @* Inactive[Plus][a] @ x
CircleTimes[a_, b_] := Inner[Composition, a, b, CirclePlus]
CenterDot[a_, b_] := Inner[Compose, a, b]


For the OP example:

I1 := {{0&, (x # + D[#,x]&)}, {0&, -D[#,x]&}}
I2 := {{0&, x^2 #&}, {0&, D[#,x]&}}
u = {f[x], g[x]};


The same results as @Jens:

CenterDot[I1, u] //TeXForm


$\left\{g'(x)+x g(x),-g'(x)\right\}$

CenterDot[I2, u] //TeXForm


$\left\{x^2 g(x),g'(x)\right\}$

CenterDot[CircleTimes[I1, I2], u] //TeXForm


$\left\{g''(x)+x g'(x),-g''(x)\right\}$

I3 := CirclePlus[CircleTimes[I1, I2], CircleTimes[I2, I1]]
CenterDot[I3, u] //TeXForm


$\left\{g''(x)+x^2 \left(-g'(x)\right)+x g'(x),-2 g''(x)\right\}$

Id := {{1&, 0&}, {0&, 1&}}
CenterDot[CirclePlus[Id, Id], u] //TeXForm


$\{2,2\}$

• Very elegant (+1). It seems Compose was once a valid command...
– Jens
Jun 1, 2017 at 18:45
• It was deprecated in M2 (and replaced with Composition), but never removed. I wish it were still documented. Jun 1, 2017 at 18:48
• Interesting - there isn't really any consistency between use of nouns versus verbs for operations. E.g., Add doesn't exist, but Subtract does.
– Jens
Jun 1, 2017 at 18:52

You can make use of my DifferentialOperator paclet to do this. Install with:

PacletInstall["https://github.com/carlwoll/DifferentialOperator/releases/download/0.1/DifferentialOperator-0.0.3.paclet"]


<<DifferentialOperator


Here is an animation defining I1 and I2 The other definition is:

u = {f[x], g[x]};


Now, let's check:

I1 . u //TeXForm


$\left\{g'(x)+x g(x),-g'(x)\right\}$

I2 . u //TeXForm


$\left\{x^2 g(x),g'(x)\right\}$

I1 . I2 . u //TeXForm


$\left\{g''(x)+x g'(x),-g''(x)\right\}$

(I1 . I2 + I2 . I1) . u //TeXForm
`

$\left\{g''(x)+x^2 \left(-g'(x)\right)+x g'(x),-2 g''(x)\right\}$