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

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

CLICK HERE to find out more related problems solutions.

Leave a Comment

Your email address will not be published.

Scroll to Top