R语言的集成

Question

我试图计算在下面R-代码中给出的函数的1和一些截止的“切”之间的积分，作为'int‘。这取决于在进行集成之前定义的两个参数dMi和dLambdaj，对于每一对参数，我都将结果保存在向量‘vec’中：

vec = c() #vector for INT values: this is our goal
dM = seq(from = 0, to = 3, by = 0.01) #vector for mass density parameter
dLambda = seq(from = -1.5, to = 3, by = 0.01) #vector for vacuum energy density parameter

for (i in 1:length(dM)) {
  for (j in 1:length(dLambda)) {

    int = function(x) ((dM[i]*x^4*(x - 1) + dLambda[j]*x^2*(1 - x^2) + x^4)^(-1/2))
    cut = 30
    INT_data = integrate(int, 1, cut)
    INT = INT_data$value
    vec = c(vec, INT)
  }
}

但是当我运行这个脚本时，我会得到一个错误：“集成中的错误(int，1，cut)：非有限函数值”。尽管如此，如果我尝试了以下代码

int = function(x) ((0*x^4*(x - 1) -1.5*x^2*(1 - x^2) + x^4)^(-1/2))
cut = 30
INT_data = integrate(int, 1, cut)
INT = INT_data$value
vec = c(vec, INT)

我得到了正确的结果，没有任何错误。所以上面的错误是不正确的，它可以计算积分，但是如果我用2‘for’-循环，似乎R不能算出它。我如何重写代码，以便为我想要的dMi和dLambdaj计算所有不同的值？

Answer 1

您的函数只为dM和dLambda的某些值定义。您可以使用try()函数尝试计算，但不能在发生错误时停止。

预先分配对象来保存结果的效率也要高得多；运行vec = c(vec, INT)会逐渐增长，这是非常慢的，因为R需要继续创建新的向量，只需要比最后一个元素多一个元素。

此代码修复了这两个问题，然后绘制结果：

dM <- seq(from = 0, to = 3, by = 0.01) #vector for mass density parameter
dLambda <- seq(from = -1.5, to = 3, by = 0.01) #vector for vacuum energy density parameter
result <- matrix(NA, length(dM), length(dLambda))

for (i in 1:length(dM)) {
  for (j in 1:length(dLambda)) {

    int <- function(x) ((dM[i]*x^4*(x - 1) + dLambda[j]*x^2*(1 - x^2) + x^4)^(-1/2))
    cut <- 30
    INT_data <- try(integrate(int, 1, cut), silent = TRUE)
    if (!inherits(INT_data, "try-error")) 
      result[i, j] <- INT_data$value
  }
}
image(dM, dLambda, result)

编辑了添加：，这是如何工作的。如果integrate在原始代码中发出错误信号，则循环将停止。try()阻止了这一点。如果没有错误，则返回integrate调用的结果。如果存在错误，则返回包含错误信息的对象。该对象有类"try-error"，所以检查if (!inherits(INT_data, "try-error"))基本上是在问“有错误吗？”如果出现错误，什么都不会发生，result的条目在初始化时被保留为NA。然后，循环继续尝试下一个dM，dLambda对。

​

​

Answer 2

这个问题是数学问题，而不是编码问题。该函数没有为您正在集成的整个域定义。当dM1 =0和dLambda > 1时，您的表达式

(dM[i]*x^4*(x - 1) + dLambda[j]*x^2*(1 - x^2) + x^4)^(-1/2)

简化为

(dLambda[j] * x^2 * (1 - x^2) + x^4)^(-1/2)

让我们把dLambdaj设为1.01，这是您的计算停止的地方：

(1.01 * x^2 * (1 - x^2) + x^4)^(-1/2)

这就是

(1.01 * x^2 - 1.01 * x^4 + x^4)^(-1/2)

或

(1.01 * x^2 - 0.01 x^4)^(-1/2)

现在，您正在计算1到30之间的x。那么当x= 11时会发生什么呢？

(1.01 * 121 - 0.01 * 14641)^(-1/2)

这只剩下你了

(122.21 - 146.41)^(-1/2)

这相当于

1/sqrt(-24.2)

因此，出现此错误的原因是，您正在将一个函数集成到一个未定义的域中。

这个函数对于dM的其他值也表现不好，在范围的中间有无限个峰值，所以即使使用integrate(..., stop.on.error = F)选项也不允许您继续计算，因为您将得到无穷大的和。