R语言集成

Question

我正在尝试计算 1 和下面 R 代码中给出的函数的某个截止点 'cut' 之间的积分 'int'。它取决于在我进行积分之前定义的 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)
  }
}

但是当我运行脚本时我得到错误：“集成（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

问题是数学问题，与编码无关。该函数不是为您正在集成的整个域定义的。当 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)

现在，您正在计算 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) 选项也不允许您继续计算，因为您将得到无穷和。

Answer 2

您的函数仅为 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() 阻止了这种情况。如果没有错误，它 returns 是 integrate 调用的结果。如果有错误，它 returns 一个包含错误信息的对象。该对象具有 class "try-error"，因此检查 if (!inherits(INT_data, "try-error")) 基本上是在询问 "Was there an error?" 如果有错误，则什么也不会发生，并且 result 的条目是保留为 NA，因为它已初始化。然后循环继续尝试下一个 dM、dLambda 对。

R语言集成

Integration in R language

math

r

numeric