R语言集成答案

【问题标题】：Integration in R languageR语言集成
【发布时间】：2020-04-23 00:20:33
【问题描述】：

我正在尝试计算 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] 的所有不同值？

【问题讨论】：

您应该在循环中打印出 i 和 j 以找出导致错误的原因。可能只有某些特定值失败，即函数中分母为零的值。
@user2554330 一旦我达到 dM[i] = 0 和 dLambda[j] = 1 我就会失败。我想这是因为 x = 0 点，但我已经删除了用 reg = 0.001 改变 x -> x + reg 的可能性，我仍然收到一个错误，说产生了 NaN。但我不明白，因为使用 reg 我已经消除了 NaN 的可能性
@user2554330 好吧，似乎这些值的积分确实是不同的。但是我怎样才能避免这个错误，所以我可以继续计算呢？我的意思是，程序不会停止
您可以使用try() 函数继续跟踪错误，即INT_data <- try(integrate(int, 1, cut)); if (inherits(INT_data, "try-error")) INT_data <- NA。

标签： r math numeric

【解决方案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 对。

【讨论】：

感谢您的完整回答，但我仍有 3 个疑问：1）try() 到底在做什么？如果对于某些“x”，integrate() 以 NaN 结尾，是否只是忽略了 NaN，但继续使用相同的积分（amin，相同的 dM[i]，dLambda[j]）？ 2）我一直在寻找信息。关于“继承”，我只是发现它与类有关。它是如何工作的？ 3) 'try-error' 是为了……？
添加了一些解释。
我一直在为您在答案中得到的图形寻找某种温度图图例，但我不知道如何制作它。我的目标是获得该图形的图例，告诉我们积分的值
contour() 或 filled.contour() 可能是你想要的。

【解决方案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)

让我们以 1.01 的 dLambda[j] 为例，这就是您的计算停止的地方：

(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) 选项也不允许您继续计算，因为您将得到一个无限的总和。

【讨论】：