R中的双积分实现

Double Integral implementation in R

我一直在使用 R 中的 mosaicCalc 包处理二重积分。我无法获得第二个二重积分的正确结果。

这是产生正确结果 (pi) 的第一个二重积分的代码。

one = makeFun(1 ~ y + x)
bx.y = antiD(one(y = y, x = x) ~ y)
bx.yx = antiD(bx.y(y = 1+cos(x), x = 2) ~ x)
bx.yx(x = pi) - bx.yx(x = 0)
# [1] 3.141593

这是第二个二重积分,根据 Wolfram,其正确结果应该是 0.684853

one = makeFun(1 ~ y + x)
bx.y = antiD(one(y = y, x = x) ~ y)
bx.yx = antiD(bx.y(y = 1/2, x = sin(x)) ~ x)
bx.yx(x = 5*pi/6) - bx.yx(x = pi/6)
# [2] 1.047198

首先我需要说服自己 Mathematica 是正确的。是的,我想这源于我对权威的不信任,但这并不困难。需要认识到单位从上到下的积分只是上减下的差,所以可以将其简化为单变量问题,用R的integrate函数求解:

integrate( function(x){sin(x)-0.5},lower=pi/6,upper=5*pi/6)
0.6848533 with absolute error < 7.6e-15

所以这促使我在 'mosaicCalc' 框架中尝试相同的策略:

> one = makeFun(1 ~ y + x)
> bx.y = antiD(one(y = y, x = x) ~ y)
> bx.yx = antiD(bx.y(y = sin(x)-1/2, x = x) ~ x)
> bx.yx(x = 5*pi/6) - bx.yx(x = pi/6)
[1] 0.6848533

答案正确,但它似乎没有正确表示和保留 "functional cascade"(发明一个我以前从未听说过的术语)。我希望限制以一种反映更通用的函数调用集的方式出现,所以最终想出了一个看起来令人满意的方法:

one = makeFun(1 ~ y + x)
bx.y = antiD(one(y = y, x = x) ~ y)
bx.yx = antiD(bx.y(x=x, y = sin(x)) - bx.y(x=x,y=1/2) ~ x)
bx.yx(x = 5*pi/6) - bx.yx(x = pi/6)
[1] 0.6848533

这确实支持更复杂的函数被积函数。我没有设置 Mathematica 来与任何东西进行集成,但使用上面的 mosaicCalc 设置我得到 1.284286 for x^2 作为被积函数。你可能想检查一下。

在处理这个问题时,在我看来,积分的顺序确实颠倒了。在我对 40 年而不是 50 年前的这些问题的模糊记忆中,似乎我一直使用 dx 计算作为内部计算,但我确实意识到这是任意的。无论如何,您在第二个反导数中分配给 x 和 y 值的角色似乎并不合适。您得到的是从 x=pi/6 和 y=0.5 的下限到 x=5*pi/6, y=1)

的上限的 2D 统一积分结果
library(cubature) # package capable of 2D-integration with fixed limits
adaptIntegrate(function(x){1}, lower=c(a=pi/6, b=0.5), upper=c(a=5*pi/6,b=1 ))
$integral
[1] 1.047198

$error
[1] 0

$functionEvaluations
[1] 17

$returnCode
[1] 0