# SPECIAL DOUBLE INTEGRAL PROCEDURE.  Integrates x^2 * y^3 for x from 1 to 2
# and y from 3 to 4, using 100 points for all midpoint rule applications.

spec.int <- function ()
{ 
  midpoint.rule (spec.int.y, 3, 4, 100)
}

spec.int.y <- function (y)
{
  midpoint.rule (spec.int.xy, 1, 2, 100, y)
}

spec.int.xy <- function (x, y)
{
  x^2 * y^3
}


# GENERAL DOUBLE INTEGRATION USING THE MIDPOINT RULE.  Integrates f(x,y,...) 
# from ax to bx and ay to by using the midpoint rule with nx points along
# the x-axis and ny points along the y-axis.  Additional arguments are
# passed as additional arguments to f.

dbl.midpoint.rule <- function(f,ax,bx,ay,by,nx,ny,...)
{
  midpoint.rule ( function (y,f,ax,bx,nx,...) midpoint.rule(f,ax,bx,nx,y,...),
   ay, by, ny, f, ax, bx, nx, ...)
}


# NUMERICAL INTEGRATION USING THE MIDPOINT RULE.  The arguments are the
# function to integrate, the low and high bounds, and the number of points
# to use.  Additional arguments to pass to the function being integrated
# can follow.

midpoint.rule <- function(f,a,b,n,...)
{ 
  h <- (b-a)/n
  s <- 0

  for (i in 1:n)
  { s <- s+f(a+h*(i-0.5),...)
  }

  h*s
}


# TESTS OF DOUBLE INTEGRATION.

print( spec.int() )

print( dbl.midpoint.rule (function (x,y) x^2 * y^3, 1, 2, 3, 4, 100, 100) )

print( dbl.midpoint.rule (function (x,y) x^2 * y^3, 1, 2, 3, 4, 200, 200) )

print( dbl.midpoint.rule (function (x,y) x^2 + y^3, 1, 2, 3, 4, 100, 100) )

print( dbl.midpoint.rule (function (x,y) x^2 + y^3, 1, 2, 3, 4, 200, 200) )