Double Integral implementation in R

Question

I've working on double integrals using mosaicCalc package in R. I'm having trouble getting the correct result of the second double integral.

This is the code of the first double integral which yields to the correct result (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

This is the second double integral which, according to Wolfram, its correct result should be 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

Answer 1

First I needed to convince myself that Mathematica was correct. Yeah, I suppose that stems from my distrust of authority, but it was not difficult. It required realizing that the integral of unity from upper to lower was simply the difference of upper minus lower, so it could be reduced to a single variable problem and solved with R's integrate function:

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

So that led me to try the same strategy in the 'mosaicCalc' framework:

> 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

That got the right answer but it didn't seem to properly represent and preserve the "functional cascade" (to invent a term I've never heard used before). I wanted the limits to appear in a manner that reflected a more general functional set of calls, so eventually came up with this which seems satisfactory:

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

This does support more complex functional integrands. I don't have Mathematica set up to do that integration with anything but with the mosaicCalc setup above I get 1.284286 for x^2 as an integrand function. You might want to check.

In working this problem it did seem to me that the order of integration was reversed. In my hazy memory of these problems from 40, no 50 years ago, it seemed I had always used the dx calculation as the inner one, but I do realize that is arbitrary. At any rate the roles you assigned to the x and y values in the second anti-derivative didn't seem appropriate. You were getting the result of a 2D integration of unity from lower limits of x=pi/6 and y=0.5 to upper limits of x=5*pi/6, y=1)

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