How to make `D` apply chain rule on user-defined symbols as it does on `Dot` and `NonCommutativeMultiply`

Question

I have my own version of Dot called mydot, and it behaves quite similar to the built-in version in that it represents the non-commutative product of a string of objects. Also, mydot is supposed to support differentiation that behaves in the same way as Dot. But I am unable to get D to operate on mydot the same way D works on built-in Dot:

D[Dot[x, (x + a y), z], y]

But

D[mydot[x, (x + a y), z], y]

I tried:

Derivative[ndiff__][mydot][x__] := 
  ReplaceAll[Derivative[ndiff][Dot][x], Dot :> mydot];
D[mydot[x, (x + a y), z], y]

As you can see, the serious problem I'm having is that it is pulling the factor a outside of mydot, which is incorrect.

Trying to get this to behave properly by making definitions of Derivative fails, because the chain rule has already been applied:

Derivative[ndiff__][mydot][x__] := stuff;
D[mydot[x, (x + a y), z], y]

How do I program D and Derivative to apply chain rule to the symbol mydot correctly?

Answer 1

By default, Mathematica always differentiates function heads that contain the independent differentiation variables, unless the function head is listed in the "DifferentiationOptions"->"ExcludedFunctions" system options. So, we can get the behavior you want by adding mydot to this list:

SetSystemOptions[
    "DifferentiationOptions" -> 
    "ExcludedFunctions"-> DeleteDuplicates[
        Append[
            OptionValue[SystemOptions[], "DifferentiationOptions"->"ExcludedFunctions"],
            mydot
        ]
    ]
];

and then adding a custom differentiation rule for mydot:

mydot /: D[mydot[a__], y_] := D[Dot[a], y] /. Dot->mydot

Then, we have:

D[mydot[x, (x+a y), z], y]

mydot[x, a, z]

as desired. Here are the rest of @Jens' tests:

D[y^2 + 2 mydot[x + b y, (x + a y), z], y]
D[y^2 + 2 mydot[x, f[x + a y], z], y]
D[mydot[x, (x + a y^2), z], y, y]

2 y + 2 (mydot[b, x + a y, z] + mydot[x + b y, a, z])

2 y + 2 mydot[x, a f'[x + a y], z]

mydot[x, 2 a, z]

Answer 2

Maybe you can use the following:

ClearAll[dmydot, mydot]

Derivative[ndiff__][mydot][x__] := dmydot[ndiff][x];

dmydot /: Times[f__, dmydot[n1___, 1, n2___][x__]] := 
 mydot @@ MapAt[Times[f] &, {x}, Length[{n1}] + 1]

All I'm doing here is pull the inner derivative that came from the chain rule back into the mydot to replace the original entry at that position. This is done by defining the derivative as dmydot and giving it a rule that acts when something multiplies it. Now it can of course happen that this "something" also contains factors that were there even before the chain rule kicked in. But at least the operation of pulling such additional factors into the mydot product can never be incorrect. Whether you are willing to live with this possibility, I'm not sure...

Here are some tests:

D[y^2 + 2 mydot[x + b y, (x + a y), z], y]

(* ==> 2 y + 2 (mydot[b, x + a y, z] + mydot[x + b y, a, z]) *)

This now behaves exactly like Dot:

D[y^2 + 2 Dot[x + b y, (x + a y), z], y]

(* ==> 2 y + 2 (b.(x + a y).z + (x + b y).a.z) *)

Here is a more general function f inside mydot:

D[y^2 + 2 mydot[x, f[x + a y], z], y]

(* ==> 2 y + mydot[x, 2 a Derivative[1][f][x + a y], z] *)

I think this is all mathematically correct, and the main issue of unwanted inner derivatives in front of the product has been fixed. In addition, it also works with higher derivatives:

D[mydot[x, (x + a y^2), z], y, y]

(* ==> mydot[x, 2 a, z] *)