Computation of second derivative fails on boundaries

Question

Hi,

In order to compute the solutions of an elasticity problem on a 1D domain thanks to a manually implemented Newton's method I need to compute a function and its second derivative prior to passing it in the weak form as an initial guess. Although I have no problem computing the function, I have some issues trying to compute its second derivative near the boundaries of the domain. I'll try to reproduce what's happening in a short example:

Having a unit interval mesh, I have defined R, a vector valued function and projected it on V which is a vector function space of Lagrange finite elements. Then I compute its first derivative calling Dx(R,0). To this point everything goes smoothly and the output of the function looks as follow:

Then I basically try to compute the second derivative of this function by first projecting Dx(R,0) on V and calling again Dx():

dR=project(dR,V)
ddR=Dx(dR,0)

Which returns this function:

As you can see, there are each time wiggles near the domain's boundaries and I think it may impede the convergence of my Newton method. So far, I tried to:

• Refine the mesh near the boundaries, which led to the same result plus another set of wiggles on the interface between the refined and non-refined parts of the mesh.
• Use higher order Lagrange elements for the computation, which led to the same results
• Compute the derivative on high order elements then projecting on lower order to get some averaging but it didn't work out as well.

As none of these brought satisfying results I begin to think that computation of derivative is not as clear as I thought it was, if anyone has already encountered such problem and/or have ideas of workarounds it would be much appreciated. As this function is the initial guess of a Newton's method I think that anything that could smoothen the second derivative enough to get convergence will do the trick.

Many thanks!
Regards,

MFV.

Answer 1

Hi

This is a difficult problem because of the regularity of the piecewise solutions, see, e.g., this previous question. However, if you have your vector valued function R in CG2, then at least you can get the solution nicely in DG0

mehs = UnitIntervalMesh(50)
CG2 = FunctionSpace(mesh, "CG", 2)
DG0 = FunctionSpace(mesh, "DG", 0)
f = interpolate(Expression("x[0]*(1-x[0])"), CG2)
d2f = project(f.dx(0).dx(0), DG0)
assert numpy.allclose(d2f.vector().array(), -2) #  True

Otherwise, I think you need som boundary conditions on grad(R) to avoid the wiggles.

Comments

Hi,

First of all thank you for your reply, I understood not long after this post what was my error concerning the use of higher order elements. Before, when I used CG1, I could compute first derivative nicely but couldn't call again Dx because, of course, the result of a derivative in a CG0 space is zero. So I had to project the derivative back to CG1 before calling Dx again, which was very silly considering that it could only result in loss of information. Removing the projection and using CG2 space the second derivative computes much more nicely.

Still, I'm very interested in your suggestion of projecting in DG0 space, I'll try it in order to see if I get any improvement in convergence of my Newton method.

Anyways thank you very much for your time.