How to get accurate derivative of a exponential decaying function?

Question

Hi,

I am trying to find the first derivative to an exponential decaying function exp(-|r|). I am only looking for the scalar field to this function, not the vector field. I have modified the example in FEniCS tutorial 1.1.9 Computing derivatives to find the scalar field. I also added a Dirichlet boundary condition to handle the cusp at origo.

The code:

from dolfin import *

# mesh
point1 = Point(-2.5,-2.5)
point2 = Point(2.5,2.5)
mesh = RectangleMesh(point1,point2,16,16, 'crossed')
V_g = FunctionSpace(mesh, 'CG', 1)

# function to derivative
f0 = Expression('exp(-sqrt(x[0]*x[0] + x[1]*x[1]))')
u_k = interpolate(f0, V_g)

# exact derivative
f_exact = Expression('exp(-2*(sqrt(x[0]*x[0] + x[1]*x[1])))')
grad_exact = interpolate(f_exact, V_g)

# plotting function to derivative
plot(u_k, title="u_k")

# boundary conditions, setting cusp point to 1
def u0_boundary(x):
    tol = 1E-15
    return abs(x[0]) < tol and abs(x[1]) < tol

u0 = Expression("1")
bc = DirichletBC(V_g,u0,u0_boundary,method='pointwise')

# finding derivative 
w = TrialFunction(V_g)
v = TestFunction(V_g)
a = inner(w, v)*dx
L = inner(u_k.dx(0)*u_k.dx(0)+u_k.dx(1)*u_k.dx(1), v)*dx
grad_u = Function(V_g)
solve(a == L, grad_u,bc)

# Plot derivative
plot(grad_u, title="grad(u)")

# Dump solution to file in VTK format
file = File("grad.pvd")
file << grad_u  

# Hold plot
interactive()

# error plot
error = abs(grad_exact - grad_u)
plot(error, title="error")  
interactive()

I am not sure that the error plot is working as I expected. I am trying to see where the biggest error occurs and not just the total error for the entire field. But If I am reading this right, the projection gets better with higher order element like CG3 and beyond. Is correct or have I done some errors in my code?

Is this the correct way of finding the most accurate derivative as possible or is there a better way?

Answer 1

You can use the following approach.

x = SpatialCoordinate(mesh)
r = (x[0]**2 + x[1]**2)**.5
f = exp(-r)
g = grad(f)**2

This results in a function g that can be used in variational forms.

To obtain a Function, use projection.

g_h = project(g, V)

Comments

Thank Magne for our answer, I do not fully understand the difference between the two solutions. It seems like we are doing the same thing but using different function, or am I wrong.

Could this also be used to find the second derivate, i.e., the laplacian of exp(-r)?

---

It's not the same thing.

The approach you've outlined first approximates $\exp(-\lvert r\rvert)$ with a finite element function, then evaluates the derivative of this FE function.

The approach in the answer provides a ufl object that can be evaluated at quadrature points without any interpolation into a FE space. You can also evaluate higher order derivatives this way.