Search code examples
rintegralnumerical-integrationinfinity

Numerical Solution for Integral between -Inf and Inf: Error non-finite function value

I have a function which I would like to integrate for x between -Inf and Inf. I'm using the function integrate in R. However, I do get an error saying Non-finite function value.

random_walk_func<-function(t,A,sigma,y,x){

a1 = (2*A/(sigma))*exp((4*A*(y-x+(4*A*t)))/(sigma))

b1 = erfc((y-x+(8*A*t))/(2*sqrt(sigma*t)))

  return(a1 * b1)
}

integrate(random_walk_func, lower = -Inf , upper = Inf, t,A,sigma,y)$value


Error in integrate(random_walk_func, lower = -Inf, upper = Inf,  : 
  non-finite function value

It appears that this is most likely due to the fact that for x values towards -Inf, a1 becomes Inf whereas b1 is 0. Thus, when a1 and b1 are multiplied, the result is NaN.

Any suggestion on how to solve this kind of numerical issues?

Solution

  • There are a couple of things to point out here. Firstly, your function needs to have as its first argument the variable over which you want to integrate, so you need to rewrite your function as:

    random_walk_func<-function(x, t, A, sigma, y)
{
  a1 <- (2*A/(sigma))*exp((4*A*(y-x+(4*A*t)))/(sigma))
  b1 <- erfc((y-x+(8*A*t))/(2*sqrt(sigma*t)))
  a1 * b1
}

    Secondly, remember that this is numeric rather than symbolic integration, so you need to have values for all the other parameters you are passing to your function. I have no idea what you want these to be, so let's set them all to 1:

    t <- A <- sigma <- y <- 1

    Thirdly, it's a good idea to look at what you're integrating if your are getting infinity errors. If there are infinite values among the evaluated points, then you will get an error rather than a numeric result:

    x <- seq(-10, 10, 0.01)
 
plot(x, random_walk_func(x, t, A, sigma, y), type = "l")

    enter image description here

    We can see that we will get an excellent approximation of the integral if we choose limits of -10 and 10:

    integrate(random_walk_func, lower = -10 , upper = 10, 
          t = t, A = A, sigma = sigma, y = y)$value
#> [1] 1

    However, ultimately the reason why you are getting the error is that a1 gets monstrously large very quickly the further from the central peak that we go, and b1 becomes infintesimal. Even though their product is nearly zero, the intermediate calculations are beyond R's numerical tolerance, which is what breaks the calculation. Once a1 exceed about 10^308, R will call it Inf and a1 * b1 is therefore also Inf.

    The way round this is to calculate a1 and b1 as logs, then return their exponentiated sum. So if you do:

    random_walk_func <- function(x, t, A, sigma, y)
{
  a1 = log(2 * A / sigma) + 4 * A * (y - x + (4 * A * t)) / sigma
  b1 = log(erfc((y - x + 8 * A * t) / (2 * sqrt(sigma * t))))
  exp(a1 + b1)
}

    Then you get:

    integrate(random_walk_func, lower = -Inf, upper = Inf, 
          t = t, A = A, sigma = sigma, y = y)$value
#> [1] 1