# Initial Step Size for NIntegrate

I am using Mathematica to numerically integrate a numerically-defined function (from an interpolation and data). This function has a sequence of somewhat unevenly spaced sharp peaks which dominate the integral. When I use NIntegrate, the initial sampling doesn't pick up all of the the sharp peaks, and hence, the total value of the integral is wrong.

I would like to set the initial step size used by NIntegrate to be something small enough that it will find all of the relevant peaks. However, I would prefer not to use a numerical integration method with a fixed sampling size, as the adaptive features are beneficial in actually determining the contribution from each peak.

One option, of course, is to break the integral up into a series of sequentially evaluated integrals; however, as the peaks are unevenly spaced, this would be time-consuming. I'd much prefer setting the initial spacing if that's possible.

Any ideas? Thank you!

(The specific function to be integrated is not very important, and since it is interpolated, I did not attach it here but if there's a need for it I can provide a sample.)

• If the peaks are at x1, x2, x3,..., try NIntegrate[f[x], {x, a, x1,x2, x3, ..., b}]. Also try Integrate[f[x], {x, a, b}], if it's an InterpolatingFunction. – Michael E2 Mar 15 '15 at 21:02
• @MichaelE2 There's about 100 spikes, so putting the list in by hand is not very feasible. (Besides, I would like something that I can easily scale as I use different numerically defined functions.) They're also not truly singularities, just sharp peaks, so I'm not sure why using NIntegrate's singularity testing will help... – Lauren Pearce Mar 15 '15 at 21:20
• You can put the spikes in with code. What about the 2nd suggestion? – Michael E2 Mar 15 '15 at 21:20
• If f[x] is not an InterpolatingFunction, then perhaps MinRecursion -> 2 or whatever lowest number works. There are other parameters to play with, described in "tutorial/NIntegrateIntegrationRules" and elsewhere. – Michael E2 Mar 15 '15 at 21:25
• Using N[Integrate[ ---- ]] runs for a long time without returning a result. Changing minrecursion to a small number gives ncvb errors, because of the spikes. – Lauren Pearce Mar 15 '15 at 21:26

(If this doesn't help, then I think you will have to post your own example.)

Random data:

SeedRandom[1];
data = Sort@
Transpose[{{0, 1000}~Join~RandomReal[{0, 1000}, 10000 - 2],
RandomReal[{-5, -4}, 10000 - 100]~Join~RandomReal[{4, 5}, 100]}];

ifn = Interpolation[data]


For the integrand I'll use ifn[x]^2. It does have about 100 spikes randomly placed:

Plot[ifn[x]^2, {x, 0, 100}, PlotRange -> {0, 100000}, PlotPoints -> 200]


A naive application of NIntegrate complains:

NIntegrate[ifn[x]^2, {x, 0, 1000}]


NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {460.922}. NIntegrate obtained 44561.490415711516and 6980.389951991658 for the integral and error estimates. >>

(*  44561.5  *)


Increasing MinRecursion changes the value -- there must be a problem:

NIntegrate[ifn[x]^2, {x, 0, 1000}, MinRecursion -> 1]


NIntegrate::ncvb:...

(*  54356.8  *)


At MinRecursion -> 3 we get a NIntegrate::slwcon warning, which I take as a good sign and NIntegrate is getting some traction, and the value has crept up to 61409.. At MinRecursion -> 4, the value jumps up to 3.05019*10^8 - wow. There is another way to go...

Using the interpolating function

An InterpolationFunction is basically a piecewise polynomial, a cubic one by default. The method ifn["Grid"] tells you points where the polynomial changes. So if f[ifn[x]] is well-behave (in our case f[t] = t^2), then the following should give good results, perhaps slowly:

NIntegrate[ifn[x]^2, Evaluate@Flatten@{x, ifn["Grid"]}] // AbsoluteTiming

(*  {8.024718, 2.07559*10^9}  *)


If the spikes are infrequent, we might be able to get away skipping some of the grid points, since NIntegrate is probably using the Gauss-Kronrod rule. How much we can skip will depend on where the spikes are and the integrand. It's hard to give a good rule of thumb. The more we skip, the more likely that NIntegrate will have to recursively subdivide some parts of the domain of integration. In this particular case, skipping every 10 works and speeds the calculation up by a factor of 8.

skip = 10;
NIntegrate[ifn[x]^2,
Evaluate@Flatten@{x, 0, ifn["Grid"][[2 ;; -2 ;; skip]],
1000}] // AbsoluteTiming

(*  {0.999823, 2.07559*10^9}  *)


Back to MinRecursion

We have to bump MinRecursion up to 10 to get a good result, that's still a wee bit off.

NIntegrate[ifn[x]^2, {x, 0, 1000}, MinRecursion -> 10, MaxRecursion -> 20] // AbsoluteTiming

(*  {3.825013, 2.07558*10^9}  *)


At MinRecursion -> 7 it converges to NIntegrate's satisfaction, but the value 2.07401*10^9 is not too good.

Hope that helps.

• This is helpful, although since I have a product of interpolating functions not immediately...although I can probable merge the grids of the two...thanks! – Lauren Pearce Mar 15 '15 at 22:58
• @LaurenPearce You're welcome. Yes, something like Evaluate@Flatten@{x, Union[{0., 1000.}, Flatten@ifn1["Grid"], Flatten@ifn2["Grid"]]} might work. – Michael E2 Mar 15 '15 at 23:07

As the following code suggests, a lot of spikes in the integrand means that several different numerical answers should be obtained when representing the data by means of more or less smooth functions (in the code the function is obtained by Interpolation and its smoothness depends on the parameter named order):

a = -50;
b = 50;
Table[SeedRandom[5]; datax = RandomSample[Range[a, b, 0.1], 75];
SeedRandom[5]; datay = RandomSample[Range[a, b, .5], 75];
data = Transpose[{datax, datay}];
p = Interpolation[data, InterpolationOrder -> order];
Print[Plot[Evaluate[p[x]], {x, a, b}, Frame -> True]];
Print["int=", NIntegrate[p[x], {x, First[datax], Last[datax]}]],
{order, { 2, 3, 4 , 5, 6}}];