1
$\begingroup$

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.)

Improve this question
$\endgroup$
8
  • $\begingroup$ 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. $\endgroup$ Commented Mar 15, 2015 at 21:02
  • $\begingroup$ @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... $\endgroup$ Commented Mar 15, 2015 at 21:20
  • $\begingroup$ You can put the spikes in with code. What about the 2nd suggestion? $\endgroup$ Commented Mar 15, 2015 at 21:20
  • $\begingroup$ 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. $\endgroup$ Commented Mar 15, 2015 at 21:25
  • $\begingroup$ 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. $\endgroup$ Commented Mar 15, 2015 at 21:26

2 Answers 2

Reset to default
2
$\begingroup$

(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]

Mathematica graphics

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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! $\endgroup$ Commented Mar 15, 2015 at 22:58
  • $\begingroup$ @LaurenPearce You're welcome. Yes, something like Evaluate@Flatten@{x, Union[{0., 1000.}, Flatten@ifn1["Grid"], Flatten@ifn2["Grid"]]} might work. $\endgroup$ Commented Mar 15, 2015 at 23:07
0
$\begingroup$

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}}];
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.