# Derivative on uneven grid using interpolate

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!

• Related: mathematica.stackexchange.com/q/82679/4999. See also the tutorial FiniteElementProgramming for more examples. I'm not sure what the highest order derivative that is supported, but probably at least second. – Michael E2 Sep 17 '15 at 2:23
• The element mesh interpolation function returns the same error as interpolation function - "Interpolation on unstructured grids is currently only supported for machine numbers." so I cannot fill functional values with symbols like f[1], f[2] etc. I need to approximate derivatives somehow and then further work with them (make set of equations, solve them for f[1], f[2], ...). – user16320 Sep 17 '15 at 6:35

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["NDSolveFEM"];
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...

• This is great, thanks! I thought of something like this but you brought that to perfection. It can be even generalized to higher order, but I don't know if the error would go to zero... – user16320 Sep 17 '15 at 21:37
• Now I wanted to calculate InterpolationSym[3.061616997868383*^-17, 0.5], which is at point grid[[17]] and there are errors {} has a length of zero and no first element. >>, {First[{}]} is neither a list of replacement rules nor a valid \ dispatch table, and so cannot be used for replacing. >> and the result is {"FirstOrder" -> {a3 + 0.5 a4 + 6.12323*10^-17 a6, a2 + 3.06162*10^-17 a4 + 1. a5}, "SecondOrder" -> {{2 a6, a4}, {a4, 2 a5}}} /. First[{}] what is wrong with grid[[17]]? – user16320 Sep 20 '15 at 15:38