时间戳

首页 > 解决方案 > R语言集成

r - R语言集成

问题描述

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

标签： rmathnumeric

解决方案

Your function is only defined for some values of dM and dLambda. You can use the try() function to attempt evaluation, but not stop in case an error occurs.

It's also a lot more efficient to pre-allocate the object to hold the results; running vec = c(vec, INT) gradually grows it, and that's very slow, because R needs to keep creating new vectors just one element longer than the last one.

This code fixes both issues, and then plots the result:

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)

Edited to add: Here's how this works. If integrate signals an error in your original code, the loop will stop. try() prevents that. If there's no error, it returns the result of the integrate call. If there is an error, it returns an object with information about the error. That object has class "try-error", so the check if (!inherits(INT_data, "try-error")) is basically asking "Was there an error?" If there was an error, nothing happens, and that entry of the result is left as NA, as it was initialized. The loop then goes on to try the next dM, dLambda pair.

推荐阅读