为什么在 R 中使用 Composite Simpson 规则（数值积分）我的近似值太大？

Question

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

,

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

为此，我在 R 中将 Composite Simpson 规则实现为一个函数，它以函数 (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

Answer 1

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

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)))
}

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

为什么在 R 中使用 Composite Simpson 规则（数值积分）我的近似值太大？

Why is my approximation too large using Composite Simpson's rule in R (numerical integration)?

r

numerical-methods

numerical-integration

simpsons-rule