为什么使用 R 中的复合辛普森规则数值积分我的近似值太大？

如何解决为什么使用 R 中的复合辛普森规则数值积分我的近似值太大？

我正在尝试使用 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,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,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,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

解决方法

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

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,t2,100000)))
}

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