area - Hard Coding a Double Reimann Sum in R

Question

I am trying to code up a double Reimann sum in R for the function f(x,y) = x^5 * y^2 where 0 < x < 2 and 4 < y < 6. This is what I have so far, but I am having a hard time finding the correct answer since this code outputs an incorrect value.

a <- 0 # lower bound for x
b <- 2 # upper bound for x
c <- 4 # lower bound for y
d <- 6 # upper bound for y
xsize <- (b-a)/1000 # step size in x-direction
ysize <- (d-c)/1000 # step size in y-direction
x <- seq(a + (xsize/2), b - (xsize/2), xsize)
y <- seq(c + (ysize/2), d - (ysize/2), ysize)

f <- x^5 * y^2

ans <- sum(f*xsize*ysize)
ans

question from:https://stackoverflow.com/questions/65830100/hard-coding-a-double-reimann-sum-in-r

Answer 1

f <- outer(x, y, function(x, y) x^5 * y^2)
ans <- sum(f)*xsize*ysize
ans
#> [1] 540.4438

Although obviously in this case you could do the x and y sums independently and not have to do the outer product.