# Derivative of an interpolated function on an unstructured grid

I'm writing my own FEM Poisson equation solver in C++ and I'm using Mathematica as a benchmark - currently it seems that my program outputs a correct solution, however I want to check whether the laplacian of the result is equal to the right-hand side. Unfortunately, the interpolation itself raises a warning "Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All. Order will be reduced to 1". FEM deals with mostly unstructured grids, so as a result, the order of interpolation is now 1 (piecewise linear functions). It is no surprise that when I try to compute the laplacian of the result, I get flat zero function.

For the sake of simplicty, let's use the following code

ClearAll[x, y]
c = {}; n = 10000;
For[i = 0, i <= n, i++,
p = {x -> Random[], y -> Random[]};
AppendTo[c, {x, y, Sin[5 x] Sinh[5 y] + 1/6 x^3 - 1/12 y^4} /. p];
]


The following

f = Interpolation@c;
lapf = D[f[x, y], {x, 2}] + D[f[x, y], {y, 2}];
Plot3D[lapf, {x, 0, 1}, {y, 0, 1}]


gives flat zero function, although the laplacian of Sin[5 x] Sinh[5 y] + 1/6 x^3 - 1/12 y^4 is x - y^2. What can be done in case I have only numerical data on an unstructured grid and I need to evaluate laplacian at some point?

• What about choosing a nice, smooth function $u$, compute the Laplacian $f = \Delta u$ of it analytically, sample it (so that you obtain a discrete function $f_h$) and compare the degrees of freedom of the FEM solution $u_h$ (i.e., $A_h \, u_h = M_h \, f_h$ with stiffness matrix $A_h$ and mass matrix $M_h$) to $u$ sampled on the grid? That's what I usually do to debug my solvers... Jan 7, 2018 at 22:41
• I see what you're suggesting...it seems like a feasible approach I might want to try. Jan 7, 2018 at 23:04
• Don't forget to pick the right interpolation order before taking derivatives. Mar 15, 2018 at 14:21
• Have you tried the Obtuse package?: familydahl.se/mathematica/#Obtuse Jul 13, 2018 at 15:35
• If noise is a concern you may want to have a look at mathematica.stackexchange.com/a/10997/1089 and mathematica.stackexchange.com/a/176375/1089 Jul 13, 2018 at 19:33

You could import your higher order mesh (2nd order is possible) and construct an ElementMesh from the data like so:

Needs["NDSolveFEM"]
mesh = ToElementMesh[
"Coordinates" -> {{0., 0.}, {1., 0.}, {1., 1.}, {0., 1.}, {0.5,
0.}, {1., 0.5}, {0.5, 0.5}, {0.5, 1.}, {0., 0.5}},
"MeshElements" -> {TriangleElement[{{1, 2, 3, 5, 6, 7}, {3, 4, 1,
8, 9, 7}}]}];


Construct some example values.

values = Function[{x, y}, x*y] @@@ mesh["Coordinates"];


Create the interpolating function.

if = ElementMeshInterpolation[{mesh}, values];


Plot the result.

Plot3D[if[x, y], {x, y} \[Element] mesh]


However, be careful with computing second order derivatives for that; this is more of a test of the interpolating function that it is for the PDE solver you write.

To test the solver you use an somewhat arbitrary expression and insert that in the differential equation. That will give you a right hand side. Solve the PDE with this new right hand side and check how close the solution is to the somewhat arbitrary expression you used to create the right hand side. This process is called manufactured solution.

First, instead of procedural code to produce your sample data, how about using functional programming on a smaller sample of 200 points:

data = Apply[
Function[{x, y},
{x, y, Sin[5 x] Sinh[5 y] + (1/6) x^3 - (1/12) y^4}
], RandomReal[{0, 1}, {200, 2}], 1]


Then, as the warning message suggests, try

f = Interpolation[data, InterpolationOrder -> All]


A plot of f shows that it agrees well with your original function. Next compute the Laplacian

lapf = Function[{x, y}, Evaluate[Laplacian[f[x, y], {x, y}]]]


which also agrees very well with the Laplacian of f.

• G'day, (1) f["InterpolationOrder"] gives {9223372036854775806, 9223372036854775806} which seems remarkable; is Mathematica trying to interpolate over the whole domain? (2) In Version 11.0.1, I get a blank plot for f and lapf, please see here. This does not happen if InterpolationOrder -> All is omitted. (3) If I try it with 2000 points, and attempt to evaluate say lapf[0.5,0.5], it gives the warning "The interpolation points {<>} are not poised, so an interpolating polynomial of total degree 62 could not be found". What is going on? Cheers!
– user40265
Apr 11, 2018 at 4:32