R语言集成

Question

我正在尝试计算以下R代码中以int形式给出的函数的1和某些截断值“ cut”之间的整数。它取决于我进行积分之前定义的2个参数dM [i]和dLambda [j]，对于每一对，我将结果保存在向量'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)
  }
}

但是当我运行脚本时，出现错误：“ integrate（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似乎无法解决。如何重新编写代码，以便可以计算出我想要的dM [i]和dLambda [j]的所有不同值？

Answer 1

之前定义的2个参数dM [i]和dLambda [j]。

您的功能仅针对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

该问题是数学问题，而不是与编码有关。该功能未针对您要集成的整个域进行定义。在dM [1] = 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)
所以让dLambda [j]处于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)
现在，您正在评估x在1到30之间。那么当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)选项也将不允许您继续计算，因为您将获得无限的总和。