Re: Points sampled by NIntegrate
- To: mathgroup at smc.vnet.net
- Subject: [mg67475] Re: Points sampled by NIntegrate
- From: "Andrew Moylan" <andrew.j.moylan at gmail.com>
- Date: Tue, 27 Jun 2006 03:14:51 -0400 (EDT)
- References: <e7g9mj$2na$1@smc.vnet.net><e7no5g$97s$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Ah OK, I understand what is happening now. Thanks for your clear explanation Anton. Cheers, Andrew antononcube wrote: > NIntegrate transforms (on the fly) the integrand > > (1/(r^2 + Sin[r]) > > into > > 1/((1 - r)^2*(1/(1 - r)^2 + Sin[1/(1 - r)])) > > using this transformation > > {rules, jac} = {{r -> 1/(1 - r)}, 1/(1 - r)^2}; > > I.e. let > > g[r_] := Evaluate[(1/(r^2 + Sin[r]) /. rules) jac] > > then > > In[22]:= NIntegrate[g[r], {r, 0, 1}] > Out[22]= 0.816047 > > If we plot g[r] we will see that it oscillates around 1. This can be > seen from these plots: > > Plot[g[r], {r, 0, 1}, PlotPoints -> 1000] > Plot[g[r], {r, 0.9, 1}, PlotPoints -> 1000] > Plot[g[r], {r, 0.99, 1}, PlotPoints -> 1000] > > We can also plot the sampling points used by these commands: > > samplingPoints = Reap[NIntegrate[g[r], {r, 0, 1}, EvaluationMonitor :> > Sow[r]]][[2, 1]]; > Length[samplingPoints] > ListPlot[Transpose[{samplingPoints, Range[Length[samplingPoints]]}]] > > (Here I used EvaluationMonitor, and Sow and Reap, instead of AppendTo > inside the integrand.) > > >From the last plot, one can see how NIntegrate`s global adaptive > algorithms samples the integrand. > > When the precision goal is increased from the 4 to 5, NIntegrate's > global adaptive algorithm needs deeper recursion level of interval > bisection. When this recursion level reaches the value given to the > option SingularityDepth (which is by default 4), the singularity > handling mechanism kicks in. The singularity handling transformation > has the property to cluster the sampling points around the singular > point (see http://library.wolfram.com/infocenter/Conferences/5365/). > This is the reason we get large values for precision goal 5. > for this integrand though, the reason for the slower convergence is not > a singularity at 1, it is the oscillatory nature of the integrand > around 1. So the singularity handling is not necessary. If we prevent > the singularity handling by giving a large value to SingularityDepth, > we get result with smaller number of sampling points, that have smaller > magnitude: > > In[66]:= sampledPoints = {}; > f[r_?NumericQ] := (AppendTo[sampledPoints, r]; 1/(r^2 + Sin[r])) > NIntegrate[f[r], {r, 1, Infinity}, SingularityDepth -> 1000, > PrecisionGoal -> 5] > Take[Sort[sampledPoints], -10] > Length[sampledPoints] > > Out[68]= 0.816048 > > Out[69]= {138.669, 170.539, 251.341, 260.339, 341.078, 502.682, > 682.156, 1005.36, 2010.73, 4021.45} > > Out[70]= 143 > > > Anton Antonov, > Wolfram Research, Inc. > > andrew.j.moylan at gmail.com wrote: > > Hi, > > > > [I apologise for the length of this post; I have tried to state my > > problem as clearly as possible, and the result is long.] > > > > I am trying to understand how NIntegrate decides at which points to > > sample a function when estimating its improper integral over e.g. {1, > > Infinity}. Please consider the following code: > > > > sampledPoints = {}; > > f[r_?NumericQ] := > > ( > > AppendTo[sampledPoints, r]; 1/(r^2 + Sin[r]) > > ) > > NIntegrate[f[r], {r, 1, Infinity}, PrecisionGoal -> 4] > > Take[Sort[sampledPoints], -10] > > Length[sampledPoints] > > > > The function f is just 1/(r^2 + Sin[r]), except it keeps track of the > > points at which it is evaluated, in a list called sampledPoints. After > > NIntegrate is called on f over {1, Infinity}, the last two lines > > returns the 10 largest values of r for which f[r] was evaluated, and > > the number of values of r for which f[r] is evaluated. In the case > > above, the following are returned: > > > > {69.3345, 85.2695, 125.67, 130.17, 170.539, 251.341, 341.078, 502.682, > > 1005.36, 2010.73} > > > > and > > > > 99. > > > > Now consider the following code, which is identical except that the > > PrecisionGoal option to NIntegrate has been increased from 4 to 5: > > > > sampledPoints = {}; > > f[r_?NumericQ] := > > ( > > AppendTo[sampledPoints, r]; 1/(r^2 + Sin[r]) > > ) > > NIntegrate[f[r], {r, 1, Infinity}, PrecisionGoal -> 5] > > Take[Sort[sampledPoints], -10] > > Length[sampledPoints] > > > > The following are now returned by the last two lines: > > > > {1005.363620876827`, 1518.2201755137608`, 2010.727241753654`, \ > > 20094.148354998237`, 34180.18350952692`, 6.859835603121333`*^7, \ > > 1.066219395284406`*^10, 1.931379671660648`*^19, 2.22745842812734`*^55, > > \ > > 8.429342764501086`*^109} > > > > and > > > > 176. > > > > To me, sampling 176 points instead of 99 seems perfectly reasonable; > > but I can't understand the presence of such large numbers as ~10^110 in > > the list of sampled points of f. Can anyone explain why they suddenly > > emerge as the PrecisionGoal increases? > > > > Now, since 1/(r^2 + Sin[r]) is not (I assume) any more expensive to > > compute at ~10^110 than at around ~1, these enormous values at which f > > is sampled do not prevent NIntegrate from converging rapidly. However, > > I wish to integrate elements of a class of functions that _are_ more > > expensive to evaluate at large values of their argument. > > > > The functions I am interested in integrating over {1, Infinity} go to > > zero as their argument goes to Infinity (as they must, if they are to > > converge) and do so monotonically. Their values at ~10^100 might as > > well be zero. By a hack, I could define them to be zero for values of > > their argument larger than some cutoff; or I could use only low options > > for PrecisionGoal in NIntegrate; but I would prefer to understand what > > it is that is making NIntegrate suddenly choose to sample the integrand > > at very large values of the argument. > > > > Thanks for any help. > > > > Cheers, > > > > Andrew.