Functional finite difference: low accuracy

Question

Finite difference approximation of functionals seems to have some unusual behavior. Usually the rule of thumb is that you can take functional finite differences down to around the square root of machine precision ($\approx 10^{-7}$) or often lower without a problem, but in my tests I observe that in Fenics, functional finite differences break down way earlier, around a stepsize of $10^{-2}$(!)

I've seen this in several situations, but heres a simple "minimal example" demonstrating the issue. One defines
$$f(u):=\int_{\partial \Omega} \nabla u \cdot n ds,$$
and computes the finite difference,
$$f'(u) \delta u \approx (f(u+s\delta u) - f(u))/s.$$
The dolfin code is as follows:

from dolfin import *

mesh = UnitSquareMesh(30,30)
V = FunctionSpace(mesh,'Lagrange',1)
u = Function(V)

n = FacetNormal(mesh)
f = dot(nabla_grad(u), n)*ds

u0 = project(Expression('sin(2*pi*x[0]+1)+cos(2*pi*x[1]-1)+x[0]*x[1]+x[1]'),V) #Some random function I made up
f0 = replace(f,{u:u0})

delta_u = project(Expression('1.3-tan(x[0])+sin(2*pi*x[0])*exp(x[1])'),V)
set_log_level(50)
print 's', 'form','diff'
for k in range(9):
    s = pow(10,-k)
    u1 = project(u0 + s*delta_u)
    f1 = replace(f,{u:u1})
    df_diff = (assemble(f1) - assemble(f0))/s

    df = derivative(f,u,delta_u)
    df_form = assemble(df)

    print s, df_form, df_diff

When run it outputs the following:

Solving linear system of size 961 x 961 (PETSc Krylov solver).
Solving linear system of size 961 x 961 (PETSc Krylov solver).
s form diff
1 -2.61447044141 -2.61445236177
0.1 -2.61447044141 -2.61439465518
0.01 -2.61447044141 -2.61379277866
0.001 -2.61447044141 -2.60777870884
0.0001 -2.61447044141 -2.54763865842
1e-05 -2.61447044141 -1.94623822025
1e-06 -2.61447044141 4.06776613726
1e-07 -2.61447044141 64.2078098556
1e-08 -2.61447044141 665.608249162

In the right column you can see the accuracy of the finite difference deteriorating as the stepsize is made smaller, far before one would expect it to.

What is going on here?

Thanks for any help.

Nick Alger · Answer 1 · Apr 4, 2015

Figured out the answer to my own question.

It appears that project() is only linear up to single-precision floating point accuracy. That is,

project(u + project(v-u)) - project(v) $\approx$ 1e-8 (rather than 1e-16).

For example, the following code,

from dolfin import *
mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh,'Lagrange',1)

u0 = project(Expression('x[0]*x[1]'),V)
delta_u = project(Expression('x[0] + 2*x[1]'),V)

print assemble(project(project(u0 + delta_u,V) - u0)*dx) - assemble(delta_u*dx)

outputs:

5.79858969729e-8

Small stepsizes in finite differences make this loss of accuracy even worse, leading to the issues seen in the original question.

To remedy the problem, essentially you have to

be very careful about when and where you call "project" and "interpolate", either not calling it at all if possible, or delaying it until the last possible moment.
add functions in the same space by adding their vectors, rather than adding the functions then projecting. Ie.,

u1 = Function(V)
u1.vector()[:] = u0.vector() + s*delta_u.vector() # <--- GOOD

u1 = project(u0 + s*delta_u, V) # <-- BAD

Making these changes to remove calls to project() in the original code fixes the problem.

Functional finite difference: low accuracy

1 Answer