CODEWITHSUNDEEP
r
RobH
RobH

Reputation: 13

Problems calculating a double integral of a product of functions in R

I'm trying to calculate a double integral with the function "integral2" from the R package "pracma". I'm having issues calculating 

integral2(function(x,y){ X(x)*R(x,y)*X(y) }, 0, 10, 0, 10)

where 

X <- function(t) {
  -0.4*sqrt(2)*sin(pi*1*t)+0.016*sqrt(2)*sin(pi*2*t)-0.01*sqrt(2)*sin(pi*3*t)
}

and

R <- function(x,y){(1/2*(x^2-x+1/6))*(1/2*(y^2-y+1/6))-
           (1/24*((abs(x-y)^4)-2*(abs(x-y)^3)+(abs(x-y)^2)-1/30))}.

My result for the double integral in r is 

integral2(function(x,y){ X(x)*R(x,y)*X(y) }, 0, 10, 0, 10)$Q = 80.77929,

but if I calculate the same integral in Maple, the result is 87.911.

Upvotes: 1

Views: 128

Answers (1)

Rui Barradas
Rui Barradas

Reputation: 76460

The following will work.

X <- function(t){-0.4*sqrt(2)*sin(pi*1*t)+0.016*sqrt(2)*sin(pi*2*t)-
    0.01*sqrt(2)*sin(pi*3*t)}

R <- function(x,y){(1/2*(x^2-x+1/6))*(1/2*(y^2-y+1/6))-
    (1/24*((abs(x-y)^4)-2*(abs(x-y)^3)+(abs(x-y)^2)-1/30))}

f <- function(x, y){X(x)*R(x, y)*X(y)}

integral2b <- function(f, lower, upper){
  integrate( function(y) {
    sapply(y, function(y) {
      integrate(function(x) f(x,y), lower[1], upper[1])$value
    })
  }, lower[2], upper[2])
}

integral2b(f, c(0, 0), c(10, 10))
#84.94517 with absolute error < 0.0081

See R-Help, this answer was adapted from that thread.

Upvotes: 1

Related Questions