Computing with matrix differential operators

Question

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.

Answer 1

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[1][g][x], -Derivative[1][g][x]} *)

CenterDot[I2, u]

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

CenterDot[CircleTimes[I1, I2], u]

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

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

CenterDot[I3, u]

(*
==> {x Derivative[1][g][x] - 
  x^2 Derivative[1][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.

Answer 2

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\}$

Answer 3

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"]

and load with:

<<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\}$