3

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

First double integral

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

Second double integral

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
r2evans
  • 141,215
  • 6
  • 77
  • 149
Dessel
  • 33
  • 4
  • 2
    Welcome to SO, Dessel! This question is lacking a lot of information for us to provide any real assistance. Please make this question *reproducible*. This includes sample code (including listing non-base R packages), sample *unambiguous* data (e.g., `dput(head(x))` or `data.frame(x=...,y=...)`), and expected output. Refs: https://stackoverflow.com/questions/5963269, https://stackoverflow.com/help/mcve, and https://stackoverflow.com/tags/r/info. – r2evans Dec 05 '19 at 00:13
  • 3
    Now that you've added the picture ... what have you tried so far? (Saying *"with no success"* doesn't help us much. If you get errors, include code and the errors and we can help. If you get warnings and/or the wrong results, we need the code, the actual output, and why it is incorrect.) – r2evans Dec 05 '19 at 00:15
  • 2
    I updated my question with the correct integral and code. I apologize for the inconvenience. – Dessel Dec 05 '19 at 00:31
  • 1
    Much better, thank you Dessel! – r2evans Dec 05 '19 at 00:32

1 Answers1

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
IRTFM
  • 258,963
  • 21
  • 364
  • 487