# [R] Double integration - Gauss Quadrature

Earl F. Glynn efg at stowers-institute.org
Mon Sep 29 22:44:14 CEST 2008

"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
>

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