ERROR, The integral is probably divergent-CodePudding

Haii guys, i have problem with code in R i wanna running the integral code, but the code its doesnt work. is there anyone can help me?

retensi = 1136074
b = 1/1.230219e-07  
sx = function(x)
{exp(-x/b)}
integrate(sx, retensi, Inf)

when i run this code, the ruslt is error Error in integrate(sx, retensi, Inf) : the integral is probably divergent i dont get it, where the mistake that i made can anyone help me? Please tell me detail of the correct code guys

the pov is how to solve this integral in R studio

CodePudding user response：

I think the problem is that the expression exp(-x/b) has such a small value at higher values of x that you are running into the limits of floating point arithmetic. For example:

retensi = 1136074
b = 1/1.230219e-07 
sx <- function(x) exp(-x/b)

sx(1e9)
#> [1] 3.734803e-54
sx(1e10)
#> [1] 0

In fact, a quick manual binary search shows that numbers above 6,056,915,224 will return 0

sx(6056915224)
#> [1] 4.940656e-324

sx(6056915225)
#> [1] 0

This means that you will get the best approximation of you integral if the upper limit of your integral is set to 6,056,915,224:

integrate(sx, retensi, 6056915224) 
#> 7068377 with absolute error < 19

We can confirm this is right by simply finding the indefinite integral of your expression, which is: -b e^(-x/b) c and noting that this is 0 when x is infinite, so the manual calculation is:

0 - (-b * exp(-retensi/b))
#> [1] 7068377

And if we are still unsure, we can confirm this in Wolfram Alpha

CodePudding user response：

1) Separate into two parts Sum the integral from retensi to a plus the integral from a to Inf for an a that works. We can find such as a by trying 10^i for i = 7, 8, ... The code stops at the first a that works.

retensi <- 1136074
b <- 1/1.230219e-07  
sx <- function(x) exp(-x/b)
for(a in 10^(7:12)) {
  res <- integrate(sx, a, Inf, stop.on.error = FALSE)
  if (res$message == "OK") break
}

a
## [1] 1e 09

res
## 5.21431e-51 with absolute error < 9.7e-51

so at a value of a = 10^9 the integral from a to Inf is essentially zero so we can just calculate the integral from retensi to a = 10^9

res <- integrate(sx, retensi, a); res
## 7068377 with absolute error < 0.0044

2) Symbolic integration We can calculate the integral symbolically

library(Ryacas0)
x <- Sym("x")
b <- Sym("b")
Integrate(exp(-x/b), x)
## yacas_expression(-(exp(-x/b) * b))

and using that symbolic result

b <- 1/1.230219e-07  
retensi <- 1136074
(-exp(-Inf/b) * b) - (-exp(-retensi/b) * b)
## [1] 7068377

3) Change of variables Another method is to use a change of variables replacing x with a function that goes to infinity when the input to that function goes to some finite value. Trying x = tan(y) we have dx = dy/cos(y)^2:

D(quote(tan(y)), "y")
## 1/cos(y)^2

The answer below checks with the above two.

b <- 1/1.230219e-07  
retensi <- 1136074
sy <- function(y) exp(-tan(y)/b) / cos(y)^2
integrate(sy, atan(retensi), pi/2)
## 7068377 with absolute error < 15