# Evaluate InterpolatingFunction of two variables on a 2D grid

I am using NDSolve to solve a system of pde's over a 2D domain: 2 equations and 2 unknowns which I call 'u' and 'v'. NDSolve returns an InterpolatingFunction for both 'u' and for 'v' and these are both functions of two variables. I can plot and visualise these solutions easily and everything looks good. Now I am simply trying to evaluate these InterpolatingFunctions at points on a 2D grid that I have generated. If I type:

vFun = v /. sol ;
vFun[1,2]


Mathematica does not evaluate the InterpolatingFunction but instead returns:

{InterpolatingFunction[Domain: {{0.,2.5},{0.,2.5}} Output: scalar]}[1,2]

I also tried:

vFun = v[x,y] /. sol ;
vFun[1,2]


What do I need to do to get it to return the numerical value?

HERE IS AN EXCERPT FROM MY CODE:

rzMax = 2.5;
Rabl = 0.5; Rnoz = 1;
zBas = 1.5;

R = ImplicitRegion[{(0 <= x <= rzMax &&
0 <= y <= zBas) || (Rabl <= x <= rzMax &&
zBas < y <= rzMax)}, {x, y}];

RegionPlot[R];
ClearAll[x, y];

sol = NDSolve[{
D[u[x, y], {x, 2}] +
D[u[x, y], {y, 2}] + (1/x) D[u[x, y], {x, 1}] ==
NeumannValue[0, x == 0] + NeumannValue[0, y == rzMax],
D[v[x, y], {x, 2}] +
D[v[x, y], {y, 2}] + (1/x) D[v[x, y], {x, 1}] ==
NeumannValue[0, x == 0] + NeumannValue[0, y == rzMax],
DirichletCondition[u[x, y] == 0, x == rzMax || y == 0],
DirichletCondition[
v[x, y] - Piecewise[{{1 - x^2, x <= Rnoz}, {0, x > Rnoz}}] == 0,
y == 0],
DirichletCondition[v[x, y] == 0, x == rzMax]},{u, v}, {x, y} \[Element] R];

Plot[v[x, y = 0] /. sol, {x, 0, 5}, PlotRange -> {0, 1}]

ContourPlot[v[x, y] /. sol, {x, 0, rzMax}, {y, 0, rzMax},PlotRange ->
{0, 1}, PlotLegends -> Automatic]

ClearAll[x, y];
vFun = v /. sol ;
vFun[1, 2];

• Could you post the main equations (or something similar to work on)? It is hard to guess the source of the problem. Dec 22 '17 at 11:18
• I updated with an excerpt from my code. Thanks for taking a look! Dec 22 '17 at 14:40
• Try vFun = v /. First[sol]. Note the extra braces {..} in what you report is returned. (Note also the use of ... /. First[s] and ... /. First[%] in the examples in the doc page for NDSolve.) Dec 23 '17 at 1:39

Here is what I get in MMA 11.1, Linux

{u1, v1} = {u, v} /. sol[[1]]
vfun[x_, y_] = v1[x, y]
vfun[0.470079, 1.49863]


0.0941901

• Thank you for your response. That worked perfect. Dec 23 '17 at 20:39

Instead of answering the main question, which arises because NDSolve, like other *Solve functions, returns a list of solutions rather than just the solution, I would rather address the more interesting objective stated in the OP:

Now I am simply trying to evaluate these InterpolatingFunctions at points on a 2D grid that I have generated.

So let's say I'm particularly interested in the solution values at a list of points, gridpts:

R = ImplicitRegion[{(0 <= x <= rzMax && 0 <= y <= zBas) ||
(Rabl <= x <= rzMax && zBas < y <= rzMax)}, {x, y}];

gridpts = Select[Flatten[Array[List, {9, 9}, {{0., 2.5}, {0., 2.5}}], 1],
RegionMember[R]];

RegionPlot[R, Epilog -> {Red, PointSize@Medium, Point@gridpts}]


If we use an NDSolve call like the OP's, then probably, the values will be interpolated, which adds an extra error to the discretization error in the NDSolve solution. Here we see the relationship of the mesh of the OP's vFun and the grid:

Show[
First[vFun]["ElementMesh"]["Wireframe"],
Graphics[{Red, PointSize@Medium, Point@gridpts}]
]


Wouldn't it be nicer to include the points in the mesh constructed by NDSolve, so that NDSolve will calculate the values I'm interested in? This can be done by manually including the points in NDSolveFEMToBoundaryMesh:

Needs["NDSolveFEM"];
bmesh = ToBoundaryMesh[
R,
"IncludePoints" -> gridpts
];

emesh = ToElementMesh[bmesh, MaxCellMeasure -> 0.01]

Show[
emesh["Wireframe"],
Graphics[{Red, PointSize@Medium, Point@gridpts}]
]


{uIF, vIF} =
NDSolveValue[{D[u[x, y], {x, 2}] +
D[u[x, y], {y, 2}] + (1/x) D[u[x, y], {x, 1}] ==
NeumannValue[0, x == 0] + NeumannValue[0, y == rzMax],
D[v[x, y], {x, 2}] +
D[v[x, y], {y, 2}] + (1/x) D[v[x, y], {x, 1}] ==
NeumannValue[0, x == 0] + NeumannValue[0, y == rzMax],
DirichletCondition[u[x, y] == 0, x == rzMax || y == 0],
DirichletCondition[
v[x, y] - Piecewise[{{1 - x^2, x <= Rnoz}, {0, x > Rnoz}}] == 0,
y == 0], DirichletCondition[v[x, y] == 0, x == rzMax]}, {u,
v}, {x, y} ∈ emesh];


Now we can extract the values of interest as follows. The list of vIF["ValuesOnGrid"] is a list of solution values at each point in emesh["Coordinates"] in the same order. We can find the indices of the grid points in the coordinate list with Nearest. It's often faster than repeatedly using Position.

myVals = Extract[vIF["ValuesOnGrid"],
Nearest[emesh["Coordinates"] -> "Index", gridpts]
];

opVals = First[vFun] @@@ gridpts;   (* values of OP's solution *)


The difference between myVals and the OP's is quite small, on the order of 10^-4 or less.

ListPlot[myVals - opVals // RealExponent]


• That's a very good idea! However, there are several use cases, where this would not be desireable, for example in multi grid methods: Think about having two resolutions of the mesh: a fine one and a coarse one and you would like to use the solution on the coarse mesh as an initial guess for interative methods on the fine one. For this, you would have to interpolate the coarse grid solution to the fine one and the point is having considerably less vertices in the coarse mesh so that the coarse mesh solution is not as expensive to compute. Dec 23 '17 at 9:41
• Thank you for your response. As far as the direct question which was to evaluate the interpolating function at the coordinates [1,2], your comment above answered it: that is, when I use vFun = v /. First[sol], it worked. Thanks also Michael E2 for your in depth response. Ultimately, I am looking to do something like that and so it was very informative. Dec 23 '17 at 20:40
• @HenrikSchumacher Thanks! There is something I don't understand in the iterative procedure: Why use gridpts for the coarse mesh, if the solution is ultimately discarded? Just use them in the fine mesh. Or in a iterative process, wait until the mesh gets fine enough that adding gridpts does not increase the computational burden appreciably. Including gridpts is only absolutely necessary in the last iteration, and only if the error from interpolating is unacceptable. Dec 23 '17 at 23:40
• @MichaelE2 Okay, that was only part of the story. For multigrid methods, you have to interpolat several times back and forth between different resolutions of the grid. The diagram on the following site shows bit of what's going on: en.wikipedia.org/wiki/Multigrid_method. Thinking of it, in order to make this procedure fast, one should not use interpolation functions at all as they have too much overhead for lookup (I know that because a collegue of mine tried ;o). Instead, one would prefer hierarchical meshes as you did (and some fine-mesh-vertex to coarse-mesh-triangle lookup). Dec 23 '17 at 23:54
• @HenrikSchumacher It might be obvious I have no real practical experience with multigrid methods, so thanks for the further comment. My main purpose was to point out a tool in Mathematica's FEM relevant to the OP's general purpose, even though I expect the answer will languish here, because the OP's question is really about something else. Dec 24 '17 at 0:17

This should be the way to fix it:

ContourPlot[v[x, y] /. sol, {x, y} ∈ R, PlotRange -> {0, 1},
PlotLegends -> Automatic]


But there still remain similar error messages, e.g., for the point {0.470079,1.49863} which lies definitely in the region R. I think, this is really a bug in NDSolve (the solution does not get interpolated up to the boundary of the domain R. Anyway, the resulting plot looks okay to me. You can suppress the error messages with Quiet if you like.

The rest can be fixed by observing that the objects returned by Solve-like functions is always a list of rules of solutions, even if there is a unique solution. The code below simply takes the first rule so that vFun becomes really a function:

ClearAll[x, y];
vFun = v /. sol[[1]];
vFun[1, 2]

• Thank you for your responses. As far as the direct question which was to evaluate the interpolating function at the coordinates [1,2], the response of Michael E2 worked perfectly: that is, when I use vFun = v /. First[sol], it worked. As they point out, there was an extra pair of brackets that I hadn't noticed. Sumit and Henrik Schumacher's suggestions also worked. Thanks Michael E2 for your in depth response as well. Ultimately, I am looking to do something like that as well and so it was very informative. Dec 23 '17 at 9:16
• @seanD If you would like to adress the other posters, you should do so directly under their post. Otherwise they won't get a notification. You can also ping one user per comment with @. Dec 23 '17 at 9:36
• +1 for fixing the real problems. This Q might fall in the "found in the documentation" close reason, but the OP has been nice. (I hope you realize he can't vote until he has 20 or 25 rep.) I must say I was confused when I first saw your answer. It was clear the first code wouldn't work, and I stopped reading. Then I reread it and realized code at the bottom fixed the list of lists problem, but I thought, "The first code still doesn't work." Later I reread again, and realized you had written it that way. But I read fast here. It seems better, imo, if the code were in the order to be executed. Dec 23 '17 at 23:49
• @MichaelE2 Admittedly, the first part of my answer was more about how things should be in a perfect world... and thus less helpful... Dec 23 '17 at 23:59