Re: Points sampled by NIntegrate
- To: mathgroup at smc.vnet.net
- Subject: [mg67464] Re: Points sampled by NIntegrate
- From: "antononcube" <antononcube at gmail.com>
- Date: Mon, 26 Jun 2006 00:13:18 -0400 (EDT)
- References: <e7g9mj$2na$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
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.