1

我正在尝试使用 R 中的数值积分来近似以下积分：

待计算的积分,

其中函数mu由以下公式定义：

在此处输入图像描述

为此，我将复合辛普森规则实现为 R 中的一个函数，它以函数 ( integrand )、积分区间 ( [a,b] ) 和所需子区间数 ( n ) 作为参数。

我已经在各种不同的数学函数上测试了我的代码，它似乎工作得很好。但是，当我尝试逼近图中所示的积分时，逼近变得很大。

我的方法是首先根据其复合辛普森近似值将内部积分定义为 R 中t的函数。然后，再次使用复合辛普森规则，以便通过将内部近似视为函数来计算外部积分被整合。

这样做时，内部近似值在自行计算时是正确的，正如预期的那样，但是整个表达式的近似值变得太大，我似乎无法弄清楚为什么。

我正在将近似值与 Maple 给出的近似值进行比较；使用t=20自己计算的内部表达式应该给出0.8157191，整个表达式应该是12.837。R 正确计算0.8157191，但给出整个表达式的32.9285 。

我尝试使用许多不同的数学函数进行简化，并使函数独立于 R 中的t，但似乎都导致相同的错误。所以，总而言之，我的问题是，为什么只有外部积分被错误地近似？

我将非常感谢任何提示或指针 - 我在此处包含了说明问题的代码：

compositesimpson <- function(integrand, a, b, n) {
  
  h<- (b-a)/n #THE DEFINITE INTERVAL IS SCALED BY
  #THE DESIRED NUMBER OF SUBINTERVALS
  
  xi<- seq.int(a, b, length.out = n+1) #DIVIDES THE DEFINITE INTERVAL INTO THE
  xi<- xi[-1]                          #DESIRED NUMBER OF SUBINTERVALS, 
  xi<- xi[-length(xi)]                 #EXCLUDING a AND b
  
  #THE APPROXIMATION ITSELF
  approks<- (h/3)*(integrand(a) + 2*sum(integrand(xi[seq.int(2, length(xi), 2)])) + 
                     4*sum(integrand(xi[seq.int(1, length(xi), 2)])) + integrand(b))
  return(approks)
  
}

# SHOULD YIELD -826.5755 BY Maple, SO THE FUNCTION IS WORKING HERE
ftest<- function(x) {
  return(exp(2*x)*sin(3*x))
}
compositesimpson(ftest, -4, 4, 100000)


# MU FUNCTION FOR TESTING
mu.01.kvinde<- function(x){ 0.000500 + 10^(5.728 + 0.038*(x+48) -10)}


#INNER INTEGRAL AS A FUNCTION OF ITS APPROXIMATION
indreintegrale.person1<- function(t){
  indre<- exp(-compositesimpson(mu.01.kvinde, 0, t, 100000))
  return(indre)
}

indreintegrale.person1(20) #YIELDS 0.8157191, WHICH IS CORRECT

compositesimpson(indreintegrale.person1, 20, 72, 100000) #YIELDS 32.9285,
#BUT SHOULD BE 12.837 ACCORDING TO MAPLE
4

1 回答 1

0

这与尝试在两个递归级别上使用向量化有关，并且它没有按照您的意愿进行操作。例如比较

indreintegrale.person1(20) 
#> [1] 0.8157191
indreintegrale.person1(c(20, 72))
#> [1] 0.8157191 0.4801160
indreintegrale.person1(72) 
#> [1] 2.336346e-10

我认为中间的答案是错误的，但其他两个是正确的。

最快的修复，做这个替换：

indreintegrale.person1 <- function(t){
  sapply(t, function(t2) exp(-compositesimpson(mu.01.kvinde, 0, t2, 100000)))
}

现在它给出了您期望的答案（但需要更长的时间来计算！）。

于 2021-02-21T20:13:37.743 回答