3

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?

Zebra Propulsion Lab
  • 1,736
  • 1
  • 17
  • 28
  • it is, easy fix is to change the tolerance: `integrate(tail, 5000, Inf,rel.tol =1*10^-5 )` since your default is too large. – J.R. Oct 15 '14 at 00:52
  • Thanks, that works! Just curious, how to decide the value of rel.tol? Can I always assign a small number, such as .Machine$double.eps, to it? – Zebra Propulsion Lab Oct 15 '14 at 04:28

1 Answers1

2

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.

J.R.
  • 3,838
  • 1
  • 21
  • 25