# Non-rectilinear integration of InterpolatingFunction

Bug introduced in 7 or earlier and persisting through 11.0.1 or later

[CASE:3487737]

I have a 2-d array that I would like to resample into a different coordinate system and integrate along one of the new coordinates. Here is an example array:

arr = GaussianMatrix[50]*ArrayPad[DiamondMatrix[5], {45, 45}, DiamondMatrix[5]];
ArrayPlot[arr, PixelConstrained -> {1, 1}]


To accomplish the resampling, I interpolate:

if = ListInterpolation[arr, InterpolationOrder -> 1, Method -> "Hermite"];
DensityPlot[if[x, y], {x, 1, 101}, {y, 1, 101},
PlotPoints -> 101, MaxRecursion -> 0
]


(The interpolation is good, even if the DensityPlot looks a bit ragged.)

It is easy to integrate the InterpolationFunction along one of the "natural" coordinates (i.e., without resampling):

int = Integrate[if[x, y], {x, 40, 60}]
(* -> InterpolatingFunction[{{1., 101.}}, <>][y] *)
Plot[int, {y, 1, 101}]


However, any change (even a trivial one) to the coordinate system causes Integrate to produce a nonsensical answer:

Integrate[if[(x + y)/2, y], {x, 40, 60}]
(* -> 20 InterpolatingFunction[{{1., 101.}, {1., 101.}}, <>][(x + y)/2, y] *)

Integrate[if[x + 1, y], {x, 40, 60}]
(* -> 20 InterpolatingFunction[{{1., 101.}, {1., 101.}}, <>][1 + x, y] *)


As the answers to this related question propose, we could work around this by re-interpolating the resampled function using FunctionInterpolation before integrating. This would work:

if2 = FunctionInterpolation[
if[(x + y)/2, y], {x, 1, 101}, {y, 1, 101},
InterpolationOrder -> 1
];
int2 = Integrate[if2[x, y], {x, 40, 60}]
(* -> InterpolatingFunction[{{1., 101.}}, <>][y] *)


The problem is that FunctionInterpolation does a very poor job of sampling the function unless InterpolationPoints is far above twice the "resolution" of the underlying interpolation. In this case over 500(!) InterpolationPoints are needed to get a reasonable reproduction of the resampled function:

• Default (11) InterpolationPoints

• 101 InterpolationPoints

• 301 InterpolationPoints

• and so on...

As you can see, the result is aliased to the point of uselessness until InterpolationPoints gets really high. With so many sampling points required, this is not an intelligent way to proceed, and there is no benefit to using FunctionInterpolation and Integrate over tabulating the function values at high resolution and summing over the resulting array.

So:

• Is it possible to construct the resampled InterpolatingFunction from the original ListInterpolation more intelligently, to obtain better results than those given by FunctionInterpolation? Or,
• Can Integrate be made to behave more reasonably when given the integration directly?
• Are the changes in coordinates always linear/affine? – Michael E2 Dec 11 '15 at 3:19
• @MichaelE2 no, the coordinate change I'm really interested in isn't. However it isn't too far from it: it transforms a skewed trapezoid into a rectangle. But if you have a solution that works well only for affine transformations, I'd still be interested to see it. – Oleksandr R. Dec 11 '15 at 3:21

This was not intended to be a self-answered question, but I found some methods that work well enough.

### Method 1: NDSolve

I had been focused on doing the resampling before the integration, and could not get NDSolve to do that (because it objects that the input is not a differential equation, which is indeed true). But it can do both at the same time:

res = NDSolveValue[{
Derivative[1, 0][t][x, y] == if[(x + y)/2, y], t[1, y] == if[1, y]
}, t, {x, 1, 101}, {y, 1, 101},
MaxStepSize -> {2, 0.5}, InterpolationOrder -> All
]
(* -> InterpolatingFunction[{{1., 101.}, {1., 101.}}, <>] *)

ByteCount[res] (* -> 202664 *)


The result seems (mostly) reasonable:

Plot[res[60, y] - res[40, y], {y, 1, 101}]


• Mathematica does all the work by itself
• Intelligent sampling and construction of InterpolatingFunctions--copes well with nonuniform (e.g. polar) resamplings (at least in principle)
• Can obtain definite integrals directly from the result without further processing

• Can't use it for resampling without integration
• Does not properly capture the structure of the original interpolation without setting MaxStepSize on the order of the feature length--note the shape of the troughs between peaks in the above figure. This is less of a problem for InterpolationOrder of the original ListInterpolation higher than 1, but that introduces artefacts into the resampled image that do not correspond to any feature of the original, and incorrectly rounds out the peaks and troughs of the integral.
• Due to small MaxStepSize, it produces a larger result, and more slowly, than seemingly it could. I admit I'm not very familiar with NDSolve Method tuning, but nothing I tried yielded any progress toward a smaller, quicker, and/or more accurate result.

### Method 2: gridding and reinterpolation

Resampling over a grid and then interpolating with ListInterpolation is much faster than FunctionInterpolation and requires fewer sampling points for a reasonable result. The main problem is fighting against the strange interpretation of Listable that InterpolatingFunction seems to have, which would be more of a nuisance for interpolations in greater than two dimensions.

pts = Compile[{},
Table[{(x + y)/2, y}, {x, 1., 101., 2.}, {y, 1., 101., 0.5}]
][];
resampled = if @@@ Transpose[pts, {1, 3, 2}]; (* unpacks, unfortunately *)

if2 = ListInterpolation[
resampled, {{1, 101}, {1, 101}},
InterpolationOrder -> 1, Method -> "Hermite"
]
(* -> InterpolatingFunction[{{1., 101.}, {1., 101.}}, <>] *)

ByteCount[if2] (* -> 167568 *)


This can be integrated normally and the result is a little better than that produced by NDSolve:

int2 = Integrate[if2[x, y], {x, 40, 60}]
(* -> InterpolatingFunction[{{1., 101.}}, <>][y] *)
Plot[int2, {y, 1, 101}]


• About twice as fast as NDSolve, and much faster than FunctionInterpolation, while being more accurate than either with fewer sampling points (in this case)
• InterpolatingFunction is not conveniently Listable, requiring a Transpose and Apply to produce the 2-d resampling. Ignoring listability to take the samples point-by-point is many times slower, so this is necessary.
The second method works well enough for my real application, but I would like to be able to do a better job with NDSolve. Does anyone have suggestions on how the sampling can be optimized with the latter, to capture the peaks and troughs successfully (without rounding) while spending fewer samples on the smooth sides of the peaks?
• @J.M. yes, me too. I have reported the misbehavior of Integrate as a bug today, and suggested that a desired result would be more like that produced by NDSolve. – Oleksandr R. Dec 15 '15 at 20:12