Can we teach Mathematica about functional differentiation?

The key relation for functional differentiation is $$\frac{\delta}{\delta f(y)}f(x)=\delta(x-y),$$ where $\delta(x-y)$ is the Dirac delta function, and the usual properties of differentiation (e.g. linearity, chain rule) still hold.

Another example: $$\frac{\delta}{\delta f(y)}(f(x)f(x))=2 \delta(x-y) f(x).$$

A possible Mathematica version could look like this

FunctionalD[f[x],f[y]]=DiracDelta[x-y]


How would one go about implementing this?

EDIT:

There's also the case of multiple arguments e.g. $$\frac{\delta}{\delta f(y)}(f(x_1)f(x_2))= \delta(x_1-y) f(x_2)+\delta(x_2-y) f(x_1).$$

-
short answer is yes; I did it in the past (with version 2.2. of mathematica too!) but more on a problem by problem basis than what you are asking. – chris Nov 13 '12 at 18:00

FunctionalD[expr_, f_[y_]] :=
With[{xs = DeleteDuplicates[Cases[expr, f[x_ /; FreeQ[x, _f]] :> x, {0, Infinity}]]},
Total[D[expr, f[#]] DiracDelta[# - y] & /@ xs]
]


edit 1: the FreeQ is needed to get the chain rule correct for f[f[x]]

edit 2: the DeleteDuplicates is needed per Stephen Luttrell's comment. I still need to think about Chris' comment.

-
 I think you need to use Cases[expr, f[x_] :> x, {0, Infinity}] // Union in order to retain only the distinct instances of f[x_]. For instance, FunctionalD[f[x] + g[f[x]], f[y]] gives twice the expected result doing it your way. – Stephen Luttrell Nov 13 '12 at 21:22 another issue may be: FunctionalD[Integrate[K[x - y] f[y], {y, 0, Infinity}], f[x]] returns the DiracDelta outside the integral. – chris Nov 13 '12 at 21:23 The other thing I completely forgot -- it's been a while since I've had to functional differentiate -- is that one should have e.g. FunctionalD[f'[y],f[x]] = -DiracDelta'[y-x]. That could probably be fixed by supplementing the pattern in the Cases to include derivatives of f, and include an corresponding appropriate derivative of the DiracDelta in the Map. – Jeremy Michelson Nov 14 '12 at 2:03 I had to downvote this because it returns incorrect results as you notice yourself in the previous comment. – Jens Nov 15 '12 at 4:45 @Jeremy: As far as I'm concerned the rule for the derivative is not overly important (in many cases one would simply get rid of the derivative through a Fourier transform). – Emerson Nov 15 '12 at 21:46

Mathematica already knows quite a lot about functional derivatives. In particular, you can do variational derivatives. That is, you have to give it the functional and the function (I would strongly suspect that your problem can be written so as to use the VariationalD function). To get started, have a look at the tutorial for the Variational Methods package.

Edit

In case someone really needs the translation from VariationalD to the functional derivative in the question, here it is. Of course, you have to load the package I mentioned above before doing anything else.

Needs["VariationalMethods"];
FunctionalD[expr_, fn_, var1_, var2_] :=
VariationalD[(expr /. var1 -> var2)  DiracDelta[var1 - var2], fn, var2]


The first argument is the expression expr to be differentiated, the second is the function fn with respect to which the derivative is taken. As the third and fourth argument, you have to provide the names of the independent variables: var1 refers to expr and var2 to fn.

An example corresponding to the discussion of Jeremy's answer :

FunctionalD[f'[x], f[y], x, y]

(* ==> Derivative[1][DiracDelta][x - y] *)


The only thing I did in FunctionalD is to multiply the given expression expr by a delta function such that the resulting integral underlying VariationalD leaves only expr evaluated at var1. The rest is done by VariationalD, no need to re-invent those manipulations.

Edit 2

I just came across another approach to defining the functional derivative without recourse to the VariationalMethods package. It's directly from the documentation for DiracDelta (under "Applications"):

FunctionalD[functional_, f_[y_]] :=
Assuming[ Element[y, Reals],
Limit[((functional /. f :>  Function[x, f[x] +
ε DiracDelta[x - y]]) -
functional)/ε, ε -> 0]]

FunctionalD[f[x]^2, f[y]]

(* ==> 2 DiracDelta[x - y] f[x] *)

FunctionalD[f[x1] f[x2], f[y]]

(* ==> DiracDelta[x2 - y] f[x1] + DiracDelta[x1 - y] f[x2] *)


This satisfies all the requirements of the question without loading a package.

-
AFAICT this doesn't handle the case of multiple arguments (see edit in question). – Emerson Nov 15 '12 at 21:42
@Emerson We can easily handle multiple arguments as soon as you understand what they are supposed to mean in terms of the variational approach. The new rule you added in your edit simply means you're varying f at y in a product of two functionals, one with x1 and the other with x2 as the integration variable. After applying the regular product rule to that, you're back to the FunctionalD in my answer. Of course one could try to automate the regular differentiation parts too - but the core of the calculation is always what I described already. – Jens Nov 16 '12 at 4:30

When I did this in FeynCalc some 17 years ago I did not care about DeltaFunction, since it is not really needed, at least not for deriving Feynman rules. You can find a couple of examples here : http://www.feyncalc.org/FeynCalcBook/FunctionalD/

So

FunctionalD[QuantumField[f]^2, QuantumField[f]]


results in 2 QuantumField[f]` The code is free, so you can just download the package and look at FunctionalD.m and friends (like RightPartialD.m, ExpandPartialD.m , etc. )

It is all a bit rusty and should be rewritten, but that takes time.

-