Derivative on uneven grid using interpolate

Question

I have a following problem: I have an uneven 2D grid of points with unknown function value. Let's take some nasty region:

region = ImplicitRegion[x^2 + y^2 - (1/2)^2 >= 0, {{x, -1, 1}, {y, -1, 1}}];

ToElementMesh will create an approximate mesh of region:

mesh = ToElementMesh[region, "MaxBoundaryCellMeasure" -> 0.05];

We can draw mesh:

MeshRegion[mesh]

mesh consists of some black magic box, but we can find a list of coordinates:

mesh[[1]]
{{0.5,-1.11022*10^-15}, {0.497508, 0.04986}, ..., {-0.701033, -0.460366}}

The next step is to feed this list with yet unknown functional values:

values = ToExpression[Table["f" <> ToString[k], {k, Length[mesh[[1]]]}]];

FEMvalues = Table[{mesh[[1, i]], values[[i]]}, {i, Length[mesh[[1]]]}]

The result is:

{{{0.5,-1.11022*10^-15}, f1}, {{0.497508, 0.04986},f2}, ..., {{-0.701033, -0.460366},f2240}

I'm interested in approximate first and second-order derivatives in grid points. Very popular and frequently advised approach would be to use Interpolation (Interpolation[FEMvalues]) and then ask for derivative in some point. The error message appears: "Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All. Order will be reduced to 1.".

Of course, first order interpolation produces piecewise linear function, so no derivative exists in grid points. Is there another way to automatically express derivative in some (bulk) grid point as expression of few neigbour functional values? I'm sure Mathematica has some built-in function to do this, but I could not find anything...thanks in advance!

Answer 1

As mentioned in the comment Mathematica has a function called ElementMeshInterpolation which come handy for FEM interpolation. However it only supports machine numbers as gird values. If you have numerical values you can use it like the following.

Needs["NDSolve`FEM`"];
region = ImplicitRegion[x^2 + y^2 - (1/2)^2 >= 0, {{x, -1, 1}, {y, -1, 1}}];
mesh = ToElementMesh[region, "MeshOrder" -> 2,MaxCellMeasure -> .0001,
"MaxBoundaryCellMeasure" -> 0.005];
solution = Sin[4 #1 + 5 #2] & @@@ mesh[[1]];
if = ElementMeshInterpolation[{mesh}, solution,"ExtrapolationHandler" ->
{Function[Indeterminate],"WarningMessage" -> False}];

However if you want to use symbolic grid values you can proceed with quadratic polynomial interpolation. To interpolate the value of the first and second order derivatives say at $p=(x,y)$, we choose five nearest neighbors of $p$. We know symbolic values of the unknown function on these six grid points (the point $p$ and five nearest ones to it). So we can get six linear equations with six unknowns. Use Solve and you are done with a basic approximation.

grid = mesh[[1]];
values = ToExpression[Table["f" <> ToString[k], {k, Length[grid]}]];
numvalues = Sin[4 #1 + 5 #2] & @@@ grid;
(* to define symbval  use numvalues if you want to test numerically *)
symbval = Dispatch[MapThread[Rule[#1, #2] &, {grid, values }]];
interpol[{x_, y_}, val_] := a6 x^2 + a5 y^2 + a4 x y + a3 x + a2 y + a1 == val;
InterpolationSym[x_, y_] := Block[{closest, pts, val, result},
   pts = grid;
   closest = Nearest[pts, {x, y}, 6];
   val = closest /. symbval;
   result = 
    First@Quiet@
      Solve[ReleaseHold@
        MapThread[Hold@interpol[#1, #2] &, {closest, val}], {a1, a2, 
        a3, a4, a5, a6}];
   {"FirstOrder" -> {a3 + 2 a6 x + a4 y, a2 + a4 x + 2 a5 y}, 
     "SecondOrder" -> {{2 a6, a4}, {a4, 2 a5}}} /. result
   ];

Test it with symbolic values.

InterpolationSym[-0.7, 0.8]

Of course it is not the best of approximation. Lets check it numerically. Second order derivatives!

D[if[x, y], {{x, y}, 2}] /. {x -> -0.7, y -> 0.8} (* FEM based *)

{{-14.9661, -18.7395}, {-18.7395, -23.4637}}

D[Sin[4 x + 5 y], {{x, y}, 2}] /. {x -> -0.7, y -> 0.8} (* Exact *)

{{-14.9126, -18.6408}, {-18.6408, -23.301}}

Our approximation will be

{"FirstOrder" -> {1.44952, 1.81188}, 
 "SecondOrder" -> {{-10.9205, -18.2583}, {-18.2583, -23.4065}}}

So you get a starting point here...