# Gateaux (directional) derivatives and higher order differentials of a functional

I would like to calculate the Gateaux derivative of a functional (i.e. a function depending on functions). A simple example for the Dirichlet functional:

$L(u(x))=\int_0^1 \frac{1}{2} (u'(x))^2 dx$

Its derivative in the direction of $v(x)$ is then calculated by:

$dL[u](v)=\left.\frac{\partial}{\partial t}L(u(x)+t v(x))\right|_{t=0}=\int_0^1 u'(x) v'(x) dx$

In Mathematica I define the functional under consideration without arguments, like that:

L[u_] := Integrate[(1/2)*D[u, x]^2, {x, 0, 1}]


and in order to compute the derivative I use the following recipe:

d[functionalHead_, functions_, testfunctions_, independentvar_] :=
Block[{t},
myop[x_, y_] := Operate[x, y, 0];
Flatten[{
Inner[myop, Flatten[{functions}], Flatten[{independentvar}], List] +
t*Inner[myop, Flatten[{testfunctions}], Flatten[{independentvar}],
List]}]] /. t -> 0]

d[L, u, v, x]


$\int_0^1 u'[x] v'[x] \, dx$

The funny list constructions in the definition of d is because I need to allow for multiple arguments that do not exclusively depend on the same space variables. Notice the the output functional is in a slightly different (more correct) form than the input functional. Now I would like to calculate higher order derivatives, here's the second one:

$d(dL[u](v))[w]=\left.\frac{\partial}{\partial t}dL[u+t w](v)\right|_{t=0}=\int_0^1 w'(x) v'(x) dx$

and higher ones are then defined recursively (in this case the Dirichlet functional is quadratic and hence third and higher order derivatives are zero), and problems appear: If I want to use my function I have to redefine the output of the first derivative so that it looks like

dL[u_] := Integrate[D[u, x] D[v[x], x], {x, 0, 1}]


then we get

dL[u]

0


but

d[dL, u, w, x]


$\int_0^1 v'[x] w'[x] \, dx$

Since I need to calculate a bulk of higher (third, fourth, fifth etc.) derivatives of some nonlinear functionals the workaround of redefining without an argument is very ineffective. Does anybody have any idea how to make the process more systematic?

EDIT: I came up with a solution. I redefined my derivative operator (now called d1) to use pure functions instead of symbols and produce a correct output which can then be used for an iteration. Here's the code:

d1[functionalHead_, functions_, independentvar_, i_: 1] :=
Block[{t,
testfunctions :=
Flatten[ToExpression[StringReplace[ToString[{functions}],
s : Except["," | "}" | "{" | "[" | "]"] .. :> s <>
ToString[i]]]]},
myop[x_,y_]:=Operate[x,y,0];
Flatten[
Inner[Function,
Flatten[{independentvar}],
Flatten[
Inner[myop, Flatten[{functions}],
Flatten[{independentvar}],List]
+t*Inner[myop,Flatten[{testfunctions}],
Flatten[{independentvar}], List]
],List]]]/.t->0];

d[functionalHead_, functions_, independentvar_, n_: 1] :=
Block[{i, tempfunc, tempfunctions, temp, dtemp},
tempfunctions = functions;
ToExpression[
ToString[temp]<>"["<>
StringReplace[ToString[{tempfunctions}],{s:"{"|
"}"..:>"",s:Except[","|"}"|"{"|"["|"]"]..:>s<>"_"}]<>"]:="
<>ToString[{tempfunctions}]<>"]"];
For[i = 1, i <= n, i++,
ToExpression[ToString[dtemp]<>"["
<>StringReplace[ToString[{tempfunctions}],{s:"{"|
"}"..:>"",s:Except[","|"}"|"{" |"["|"]"]..:>s<>"_"}]<>"]=d1["
<>ToString[temp]<>","<>ToString[tempfunctions]<>","
<>ToString[independentvar]<>","<>ToString[i]<>"];"];
ToExpression[ToString[temp]<>"["<>StringReplace[ToString[
{tempfunctions}],{s:"{"|"}"..:>"",s:Except[","|"}"|"{"|
"["|"]"]..:>s<>"_"}]<>"]="<>ToString[dtemp]<>"@@Flatten["
<>ToString[{tempfunctions}]<>"]"];];
temp @@ Flatten[{tempfunctions}]]


The recursion step is sadly very messy: I wanted that the test functions are automatically named u1, v1 for the first derivative or u2 for the second and so on. This led me to a confusing manipulation of variable names which I was unable to do in a Mathematica-way so I transformed everything to strings, did the manipulation and pulled everything back as expressions. Note the Except command is being used because I also wanted to include greek and special characters as function names.

The whole thing now works like so:

L[u_,v_,w_]:=Integrate[(1/2)*D[u[x],x]^2+(1/2)*D[v[x],x]^2,{x,0,1}]+(1/2)*w[s]^2
d[L, {u, v, w}, {x, x, s}]


$\int_0^1 \left(u'[x] \text{u1}'[x]+v'[x] \text{v1}'[x]\right) \, dx+w[s] \text{w1}[s]$

d[L, {u, v, w}, {x, x, s},2]


$\int_0^1 \left(\text{u1}'[x] \text{u2}'[x]+\text{v1}'[x] \text{v2}'[x]\right) \, dx+\text{w1}[s] \text{w2}[s]$

EDIT2: I've changed the title since it seems to be some discrepancies with respect to the terminology used in Wikipedia and in other references: to be concise, I am following the book from Ambrosetti and Prodi on nonlinear analysis.

-
I think your question will be answered by the VariationalMethods package. An example of how to use it in practice are the the Euler equation for a classical Lagrangian. Maybe your main point are higher variational derivatives. Then perhaps one can point to the question Can we teach Mathematica about functional differentiation? –  Jens Jul 3 '14 at 20:22
google.com/… –  belisarius Jul 3 '14 at 20:35
@Jens: the variational methods package takes me one step further than I would like to: the Euler equation actually come from equating the derivative to zero and then integrating by parts - but this last step of integration by parts is exactly what I must avoid... –  tks Jul 3 '14 at 20:41
@belisarius cakemusic.com –  tks Jul 3 '14 at 20:41
@A.S. I see - so the question is really about defining the higher-order derivative operators in a more automatic way, if I understand correctly. –  Jens Jul 3 '14 at 20:45

I think your definition of d is not properly generalizable because the list dimensions don't match when doing higher derivatives. So I instead use a simpler definition of the Gateaux derivative from Wikipedia which does exactly the same thing as what you're trying to do. I call it gatD, and it takes the operator, the function u and a List of test functions. The length of the latter list uniquely determines the order of the derivative. The independent variables don't need to be specified if, as you do in the question, the operator definition itself contains the assumptions about what the names of these independent variables are. Here I take your example operator L which assumes that the independent variable is x.

L[u_] := Integrate[(1/2)*D[u, x]^2, {x, 0, 1}]

gatD[oper_, func_, tFunc_List] :=
Module[{t = Array["t", Length[tFunc]]},
D[oper[func + t.tFunc], Sequence @@ t] /. Thread[t -> 0]]


Now some tests to verify that I get the same results as in your question, but without manual intervention to re-define any operators:

gatD[L, u[x], {v[x]}]


$$\int_0^1 u'(x) v'(x) \, dx$$

gatD[L, u[x], {v[x], w[x]}]


$$\int_0^1 v'(x) w'(x) \, dx$$

-
Very nice approach indeed. It seems however that is is not possible to handle more than one arguments. I worked out a solution based on my previous approach: it does not look nice at all but seems to do the trick (see edit)! Maybe there are some ideas of how to make it more readable and/or faster...(I would also upvote your answer but my rep is sadly not sufficient) –  tks Jul 6 '14 at 12:56