-Inf 和 Inf 之间积分的数值解：错误非有限函数值

Question

我想在 -Inf 和 Inf 之间为 x 集成一个函数。我在 R 中使用 integrate 函数。但是，我收到一条错误消息 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

看来这很可能是因为 x 值朝向 -Inf，a1 变为 Inf 而 b1 是 0。因此，当 a1 和 b1 相乘时，结果是 NaN.

对于如何解决此类数值问题有什么建议吗？

Answer 1

这里有两点需要指出。首先，您的函数需要将要集成的变量作为其第一个参数，因此您需要将函数重写为：

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
}

其次，请记住这是数字积分而不是符号积分，因此您需要为传递给函数的所有其他参数设置值。我不知道你想要这些是什么，所以让我们将它们全部设置为 1:

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

第三，如果您遇到无穷大错误，最好查看您正在集成的内容。如果评估点之间有无穷大的值，那么你会得到一个错误而不是一个数字结果：

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

我们可以看到，如果我们选择 -10 和 10 的极限，我们将得到一个很好的积分近似值：

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

但是，最终导致错误的原因是 a1 离中心峰越远，变得非常快，b1 变得无穷小。尽管它们的乘积几乎为零，但中间计算超出了 R 的数值容差，这就是破坏计算的原因。一旦 a1 超过 10^308，R 将称其为 Inf 并且 a1 * b1 因此也是 Inf.

解决方法是计算 a1 和 b1 作为对数，然后 return 它们的幂和。所以如果你这样做：

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)
}

然后你得到：

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

-Inf 和 Inf 之间积分的数值解：错误非有限函数值

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

r

integral

infinity

numerical-integration