这个数值积分背后发生了什么？

Question

我正在计算积分值

∫₀1 ∫₁² (r - u)r drdu

我想了解 R 中如何处理以下错误。

浏览论坛，我意识到需要对两个积分中的第一个积分（在 r 中）进行向量化来计算第二个积分（在 u 中）。因此我假设以下代码是解决计算的正确方法。

fI <- function(u) {
  integrand <- function(r, u) {
    (r - u)*r
  }
  sapply(X=u, function(z) {
    integrate(f=integrand, lower=1, upper=2, u=z)$value
  })
}

integrate(f=fI, lower=0, upper=1)

返回

1.583333 with absolute error < 1.8e-14.

根据 Mathematica，前者也是正确的值。

但是，我可以注意到，即使在没有矢量化的情况下错误地定义

fI

，一个小技巧也可以计算积分（我认为是错误的）

fI <- function(u) {
  integrand <- function(r, u) {
    (r - u) * r
  }
  integrate(f=integrand, lower=1, upper=2, u=u)$value
}

integrand2 <- function(u) {
  rep(1, length(u))*fI(u=u)
}

integrate(f=integrand2, lower=0, upper=1)

这会带来不同的结果

1.5 with absolute error < 1.7e-14.

任何人都可以解释一下第二种情况发生了什么以及在

u

函数内的

fI(u=u)

中作为

integrand

传递的值是什么？我知道这样

u

不能取一定范围的值。
我需要理解这个机制，因为我的实际计算有点复杂，我需要了解这个简化的情况下发生了什么。

Answer 1

如果将

sapply(u, message)

作为调试语句放入

integrand2()

中，您可以看到

u

值是

lower = 0

和

upper = 1

之间的样本，从中心开始，然后从两侧向中心交替。

fI <- function(u) {
  integrand <- function(r, u) {
    (r - u) * r
  }
  integrate(f=integrand, lower=1, upper=2, u=u)$value
}

integrand2 <- function(u) {
  sapply(u, message)
  rep(1, length(u)) * fI(u = u)
}

integrate(f=integrand2, lower=0, upper=1)
#> 0.5
#> 0.0130467357414141
#> 0.986953264258586
#> 0.0674683166555077
#> 0.932531683344492
#> 0.160295215850488
#> 0.839704784149512
#> 0.283302302935376
#> 0.716697697064624
#> 0.425562830509184
#> 0.574437169490816
#> 0.00217141848709596
#> 0.997828581512904
#> 0.0349212543221459
#> 0.965078745677854
#> 0.109591136706792
#> 0.890408863293208
#> 0.218621432665698
#> 0.781378567334302
#> 0.35280356864927
#> 0.64719643135073
#> 1.5 with absolute error < 1.7e-14

^{创建于 2024-06-28，使用 reprex v2.1.0.9000}