# How to speed up integration of interpolation function?

I have a list of data, and I interpolate it to a function. Then I need to do an integration with the interpolating function. But I found that the speed is unacceptably slow.

My data is here (a little long, but don't go away :) ):

data = {-0.00799443, -0.00581522, -0.00557107, -0.00543862, -0.00528042, \
-0.00508618, -0.00486091, -0.00461279, -0.00435009, -0.00408028, \
-0.0038098, -0.00354416, -0.00328805, -0.00304535, -0.00281906, \
-0.0026108, -0.00242045, -0.00224588, -0.00208332, -0.00192845, \
-0.00177801, -0.00163133, -0.00149092, -0.00136145, -0.0012474, \
-0.00115024, -0.00106716, -0.000992246, -0.000919878, -0.000848073, \
-0.000779184, -0.000717175, -0.000663667, -0.000616407, -0.00057173, \
-0.000528424, -0.000488614, -0.000454547, -0.000425288, -0.000397686, \
-0.000370268, -0.00034446, -0.000321488, -0.00030009, -0.000278782, \
-0.000258483, -0.000240725, -0.000224931, -0.000209972, -0.000196452, \
-0.000184918, -0.000174195, -0.000163592, -0.000153752, -0.000144418, \
-0.000134884, -0.000125771, -0.000117444, -0.000109436, -0.000102175, \
-0.0000959463, -0.0000902133, -0.000085125, -0.0000806452, \
-0.0000762082, -0.0000719591, -0.0000677566, -0.000063368, \
-0.0000591507, -0.0000549953, -0.0000510613, -0.0000475951, \
-0.0000444333, -0.0000417897, -0.0000394736, -0.0000373836, \
-0.0000354749, -0.0000334705, -0.0000314543, -0.0000292503, \
-0.0000269879, -0.0000247026, -0.0000224853, -0.0000204942, \
-0.0000187118, -0.0000172668, -0.0000160166, -0.0000149913, \
-0.0000139883, -0.0000129844, -0.000011827, -0.0000105289, \
-9.06132*10^-6, -7.50783*10^-6, -5.94092*10^-6, -4.46213*10^-6, \
-3.15097*10^-6, -2.0399*10^-6, -1.13236*10^-6, -3.57489*10^-7,
3.57489*10^-7, 1.13236*10^-6, 2.0399*10^-6, 3.15097*10^-6,
4.46213*10^-6, 5.94092*10^-6, 7.50783*10^-6,
9.06132*10^-6, 0.0000105289, 0.000011827, 0.0000129844, \
0.0000139883, 0.0000149913, 0.0000160166, 0.0000172668, 0.0000187118, \
0.0000204942, 0.0000224853, 0.0000247026, 0.0000269879, 0.0000292503, \
0.0000314543, 0.0000334705, 0.0000354749, 0.0000373836, 0.0000394736, \
0.0000417897, 0.0000444333, 0.0000475951, 0.0000510613, 0.0000549953, \
0.0000591507, 0.000063368, 0.0000677566, 0.0000719591, 0.0000762082, \
0.0000806452, 0.000085125, 0.0000902133, 0.0000959463, 0.000102175, \
0.000109436, 0.000117444, 0.000125771, 0.000134884, 0.000144418, \
0.000153752, 0.000163592, 0.000174195, 0.000184918, 0.000196452, \
0.000209972, 0.000224931, 0.000240725, 0.000258483, 0.000278782, \
0.00030009, 0.000321488, 0.00034446, 0.000370268, 0.000397686, \
0.000425288, 0.000454547, 0.000488614, 0.000528424, 0.00057173, \
0.000616407, 0.000663667, 0.000717175, 0.000779184, 0.000848073, \
0.000919878, 0.000992246, 0.00106716, 0.00115024, 0.0012474, \
0.00136145, 0.00149092, 0.00163133, 0.00177801, 0.00192845, \
0.00208332, 0.00224588, 0.00242045, 0.0026108, 0.00281906, \
0.00304535, 0.00328805, 0.00354416, 0.0038098, 0.00408028, \
0.00435009, 0.00461279, 0.00486091, 0.00508618, 0.00528042, \
0.00543862, 0.00557107, 0.00581522, 0.00799443}


Interpolate it:

f = Interpolation[data];
Plot[f[x], {x, 1, 200}, PlotRange -> All]


It looks like:

Now I define a function:

Clear[b];
b[x_, y_] := NIntegrate[
Cross[{0, f[s], 0}, {x - s, 0, y}][[{1, 3}]]/((x - s)^2 + y^2),
{s, 1, 200}
];


The integration is so slow!!

b[1., 2.] // AbsoluteTiming

(*{1.17772, {-0.00827965, -0.0104805}}*)


What I want to do is a vector plot:

VectorPlot[b[x, y], {x, -10, 210}, {y, -3, 3}]


But with this kind of slow integration, this is painful. Are there better ways to speed up the integration?

• – Alexey Popkov May 25 '16 at 7:32

The way to deal with this is to use the special setting Method -> "InterpolationPointsSubdivision" of NIntegrate[], which will automagically split the integrand so that an integration rule (by default, "GlobalAdaptive") is only applied within each piecewise polynomial interval of the InterpolatingFunction[] involved. This is akin to the functionality of the old package NumericalMathNIntegrateInterpolatingFunct​.

As a demonstration:

With[{x = 5, y = 5},
NIntegrate[Cross[{0, f[s], 0}, {x - s, 0, y}][[{1, 3}]]/((x - s)^2 + y^2),
{s, 1, 200}]] // AbsoluteTiming
{0.619479, {-0.00929476, -0.00291246}}

With[{x = 5, y = 5},
NIntegrate[Cross[{0, f[s], 0}, {x - s, 0, y}][[{1, 3}]]/((x - s)^2 + y^2),
{s, 1, 200}, Method -> "InterpolationPointsSubdivision"]] // AbsoluteTiming
{0.0798281, {-0.00929476, -0.00291246}}


Options[NIntegrateInterpolationPointsSubdivision] displays the suboptions that can be fed to this method.

• Hi, J.M. This option works like magic. But I still don't understand how does it work. For example Options[NIntegrateInterpolationPointsSubdivision] has option MaxSubregion->1000, what does it mean? To split the region into 1000 piece? And how does this makes interpolation function be integrated faster? Is the calculation times being reduced? – matheorem May 25 '16 at 5:57
• MaxSubregion gives the upper limit of how many subintervals should your integral be split into. As I already said, the splitting is done so that the integration methods are only applied into each of the piecewise polynomial regions; otherwise, it will waste time recursing at those breakpoints. – J. M. is away May 25 '16 at 6:02
• But why there is breakpoints? The data already been smoothed – matheorem May 25 '16 at 6:14
• It is only $C^1$ continuous, at the least. – J. M. is away May 25 '16 at 6:26
• I think I understand. So this option will automatically and exactly cut interpolation function at original data points, right? – matheorem May 25 '16 at 6:46

An alternative to using NIntegrate is to convert the integral into an ODE, and then use NDSolveValue or ParametricNDSolveValue. For example:

pf = ParametricNDSolveValue[
{
{u'[s], v'[s]} == Cross[{0,f[s],0}, {x-s,0,y}][[{1,3}]]/((x-s)^2+y^2),
u[1]==0, v[1]==0
},
{u[200], v[200]},
{s, 1, 200},
{x, y}
];


Compare:

b[1, 2] //AbsoluteTiming
pf[1, 2] //AbsoluteTiming


{1.11616, {-0.00827965, -0.0104805}}

{0.004751, {-0.00827962, -0.0104802}}

Over 200 times faster with basically the same result. VectorPlot visualization:

VectorPlot[
pf[x, y],
{x, -10, 210},
{y, -3, 3},
VectorPoints -> {14, 14}
]


If you are willing to accept some error you can get a faster result by fitting the data to function.

Other than at each extreme the data looks like a shifted and scaled Sinh function.

Using your data (i.e., not re-copied)

sol = FindFit[data,
a Sinh[(x - x0)/b], {{a, 0.000015}, {b, 15}, {x0, 100}}, x]

(* {a -> 0.0000140493, b -> 14.3721, x0 -> 100.5} *)


Using the results define

f[x_] := 1.405*10^-5 Sinh[(x - 100.5)/14.37]


Here is the comparison

Show[
ListPlot[data, PlotStyle -> Red, PlotRange -> All],
Plot[f[x], {x, 0, 200}, PlotStyle -> Black, PlotRange -> All]

]


I think I had to make sure that the inputs to b were numerical.

b[x_?NumericQ, y_?NumericQ] := NIntegrate[
Cross[{0, f[s], 0}, {x - s, 0, y}][[{1, 3}]]/((x - s)^2 + y^2),
{s, 1, 200}]


Using a fitted function enables you to pick up some speed (a factor ~ 80)

b[1., 2.] // AbsoluteTiming

(* {0.0143029, {-0.00872379, -0.0107461}} *)


and plot the vector plot (still slow, but not too painful).

If you find the error is too great you could try a Piecewise function and use the Sinh in the middle and fit a separate function at the end points. Looks like two lines on the right from 185 to 199 and 199 to 200 and on the left another two lines from 1 to 2 and 2 to 15 blended in with the Sinh.

I didn't try it so I don't know how it would affect performance.

Good luck!

• Wow, you really captured the essential feature of my data. Your approach is really fast and the result is almost the same as accurate calculation. Great insight ! What is more, I tried Piecewise but found it generally harm the performance. – matheorem May 25 '16 at 5:51
• If it's really Sinh around 100.5, the integral is essentially odd about the middle of the domain. The integral will be very very small as a result of large cancellations (and that may signal that your data is insufficient to get a good answer). – evanb May 27 '16 at 19:34

If results with lower precision goals are satisfactory (which is an assumption in couple of the answers) then using NIntegrate with a combination of option settings for Method and PrecisionGoal can produce results much quicker.

In the table below the first two options settings are the same as in J.M.'s answer. We can see that the last one is 100 times and 10 times quicker respectively.

opts = {{},
Method ->
"InterpolationPointsSubdivision", {Method -> {Automatic,
"SymbolicProcessing" -> 0}, PrecisionGoal -> 3}};

res =
Table[
With[{x = 5, y = 5, o1 = o},
Prepend[
AbsoluteTiming[
NIntegrate[
Cross[{0, f[s], 0}, {x - s, 0, y}][[{1, 3}]]/((x - s)^2 +
y^2), {s, 1, 200}, o1]], o1]
],
{o, opts}];
TableForm[res,
TableHeadings -> {None, {"Options", "AbsoluteTiming", "Result"}}]