I have a function defined as
tail <- function(x) {585.1961*x^-2.592484}
When I tried to get the integral of it from 5000 to Inf, it returns an error:
> integrate(tail, 5000, Inf)
Error in integrate(tail, 5000, Inf) : the integral is probably divergent
However, integrating from 1000 to Inf and from 1000 to 5000 both worked fine:
> integrate(tail, 1000, Inf)
0.006134318 with absolute error < 2.5e-05
> integrate(tail, 1000, 5000)
0.005661634 with absolute error < 4.9e-09
Isn't integrate(tail, 5000, Inf) simply equal to integrate(tail, 1000, Inf) - integrate(tail, 1000, 5000)? Why that resulted in a divergent integral?
Your default tolerance (.Machine$double.eps^0.25) is too large, so we change it:
> tail <- function(x) {585.1961*x^-2.592484}
> integrate(tail, 5000, Inf,rel.tol =.Machine$double.eps^0.5 )
0.0004727982 with absolute error < 1.5e-09
and this is approximately the same as:
> integrate(tail, 1000, Inf)$val-integrate(tail, 1000, 5000)$val
[1] 0.0004726847
Of course you can just set your tolerance to .Machine$double.eps, but this comes at a cost of time:
> library(microbenchmark)
> a<-function(x) for(i in 1:50) integrate(tail, 5000, Inf,rel.tol =.Machine$double.eps^x )
> microbenchmark(a(0.5),a(0.7),a(1))
Unit: milliseconds
expr min lq median uq max neval
a(0.5) 10.44027 10.97920 11.12981 11.40529 19.70019 100
a(0.7) 13.02904 13.69813 13.95942 14.89460 23.02422 100
a(1) 15.14433 15.96499 16.12595 16.38194 26.27847 100
eg. a increase in time around 50 pct.