I would like to calculate the Gâteaux 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];
(D[#, t] &)[functionalHead @@
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];
(D[#,t]&)[functionalHead@@
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[functionalHead]<>"@@Flatten["
<>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.