-2

I was doing some integration into a loop using integrate and I come up with an error I can't understand neither get rid of. Here is a MWE I could extract:

b=1/1.230219e-07
f=function(x)
{exp(-x/b)}
integrate(f,0, Inf)

this returns an error "Error in integrate(f, 0, Inf) : the integral is probably divergent" which is obviously false Error in integrate(f, 0, Inf) : the integral is probably divergent

i tried to change the upper and lower and it work

b=1/1.230219e-07
f=function(x)
{exp(-x/b)}
integrate(f,0, 2)

and the reslut 2 with absolute error < 2.2e-14

i think the problem comes from the infinity, but i dont know how to solve it?

Severin Pappadeux
  • 18,636
  • 3
  • 38
  • 64
EDITHA
  • 11
  • 2
  • 1
    It's not C++. Don't spam tags – 273K Dec 27 '22 at 03:38
  • 1
    You're doing numerical integration, not mathematical integration, so theoretical niceties (the function is mathematically convergent) don't apply. Numerical integration involves some strategy of calling the function to sample values (by calling `f()`), and adding up or averaging. The algorithm being used in `integrate()` [whatever it is] is presumably sensitive to the errors produced by each call of `f()` in your case. Bear in mind there is no universal strategy that *always* bounds errors - your function is presumably one for which the strategy used in `integrate()` fails to bound errors. – Peter Dec 27 '22 at 03:55
  • so, whats should i do? – EDITHA Dec 27 '22 at 04:05

1 Answers1

1

Use package cubature, function hcubature.
By hand the integral is equal to b and that's the result hcubature finds.
Also, I have redefined the function to take b as a formal argument.

library(cubature)

f <- function(x, b){exp(-x/b)}

b <- 1/1.230219e-07
b
#> [1] 8128634

eps <- .Machine$double.eps^0.5
hcubature(f, 0, Inf, b = b, tol = eps)
#> $integral
#> [1] 8128634
#> 
#> $error
#> [1] 0.05960908
#> 
#> $functionEvaluations
#> [1] 825
#> 
#> $returnCode
#> [1] 0

Created on 2022-12-27 with reprex v2.0.2

Rui Barradas
  • 70,273
  • 8
  • 34
  • 66