Sorry for the relatively simple question, but the doc is a bit short, and I couldn't find a corresponding QA.
Say I have the functional
$\mathcal F[u] = \int_\Omega q(u(x))\, w(x)\; dx$
where $q$ is a known nonlinear function and $w$ is a known weight. I would like to compute the derivative of $\mathcal F$ w.r.t. $u$, evaluated at $u=u_0$, for some $u_0$.
If the form is defined as
F = q(u)*w*dx,
how could this be done using derivative? Also I would say that u and w need to be defined as Function, since the form is nonlinear, is that right?
EDIT:
It seems to work fairly well in 1D, modulo normalisation, see here. But in 2D things go haywire (normalisation problem, not a big deal, but a lot of oscillations!)
Key to the plots: [theoretical gradient] [computed gradient] [ratio of the two]
# coding: utf-8
# This file is released under the WTFPL (wtfpl.net)
from __future__ import division
from dolfin import *
try:
from mshr import *
except:
pass
normalise = False
# consider the unit ball,
R = 1
N = 100
# in dimension 1 or 2.
dim = 2
if dim == 1:
mesh = IntervalMesh(N, -R, R)
else:
domain = Circle(Point(0.0, 0.0), R)
if dolfin.__version__ > '1.4.0+':
mesh = generate_mesh(domain, N)
else:
mesh = Mesh(domain, N)
V = FunctionSpace(mesh, "DG", 1)
u = Function(V)
# we will consider ∫ q(u(x)) w(x) dx,
# where w(x) = w(r) = 1+r²
w = project(Expression('x[0]*x[0]+x[1]*x[1]+1'), V)
# and q(y) = y²
q = lambda x:x*x
energy = w*q(u)*dx
# differentiate, should coefficient_derivatives only be specified if w depends on u?
deriv = derivative(energy, u) #, coefficient_derivatives={w:0})
# evaluate at some u0
u0 = project(Expression('1+cos(3*pi*x[0])/4'), V)
deriv = replace(deriv, {u:u0})
gradient_v = assemble(deriv).array()
# normalise
if normalise and dim == 1:
gradient_v = N/R * gradient_v
if normalise and dim == 2:
pass
# build the dolfin function
gradient = Function(V)
gradient.vector()[:] = gradient_v
# plot
plot(2*w*u0, title="Theoretical gradient", mode="color")
plot(gradient, title="Numerical gradient", mode="color")
plot(2*w*u0/gradient, title="Ratio", mode="color")
interactive()
