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

问题描述

我有一个功能想集成到x和-Inf之间的Inf。我在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,y)$value


Error in integrate(random_walk_func,: 
  non-finite function value

这似乎很可能是由于以下事实：对于x的{​​{1}}值，-Inf变成a1而Inf是{{1 }}。因此，当b1和0相乘时，结果为a1。

关于如何解决此类数字问题的任何建议？

解决方法

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

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,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,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,y = y)$value
#> [1] 1