-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 计算为对数，然后返回它们的指数和。所以如果你这样做：

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