[R] Double integration - Gauss Quadrature

[Earl F. Glynn]

"Susanne Pfeifer" <tiffy at tiffy.it> wrote in message 
news:48DE3BE0.5020200 at tiffy.it...
> Hi,
>
> I would like to solve a double integral of the form
>
> \int_0^1 \int_0^1 x*y dx dy
>
> using Gauss Quadrature.

OK, I felt guilty for using a for loop, so here's something that should be 
close to what you want using sapply:

library(statmod)
N <- 5
GL <- gauss.quad(N)
nodes   <- GL$nodes
weights <- GL$weights

gauss_legendre2D <- function(f, a1,b1, a2,b2, nodes, weights)
{
  C2 <- (b2 - a2) / 2
  D2 <- (b2 + a2) / 2
  y <- nodes*C2 + D2

  C1 <- (b1 - a1) / 2
  D1 <- (b1 + a1) / 2
  x <- nodes*C1 + D1

  C2*sum(weights *
      sapply( y,
          function(y) { C1 * sum( weights * f(x, y)) } ) )
}

# your problem:  area = 0.25
gauss_legendre2D(function(x,y) {x*y}, 0.0, 1.0, 0.0, 1.0, nodes, weights)

# test:  area of unit circle = integral 0..2pi integral 0..1 r*dr *dTheta
gauss_legendre2D(function(x,y) {x}, 0.0, 1.0, 0.0, 2*pi, nodes, weights)

efg
Earl F Glynn