[R] Double integration - Gauss Quadrature

[Susanne Pfeifer]

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.

I know that I can use R's integrate function to calculate it:

integrate(function(y) {
sapply(y, function(y) {
integrate(function(x) x*y, 0, 1)$value
})
}, 0, 1)

but I would like to use Gauss Quadrature to do it.
I have written the following code (using R's statmod package) which
works fine for one integral but it doesn't work for a double one:

# Gauss-Legendre abscissas
nodes <- gauss.quad.prob(25,dist="uniform",l=-1,u=1)$nodes
# and weights
weights  <- gauss.quad.prob(25,dist="uniform",l=-1,u=1)$weights
weights <- weights*2

# Approximate integral of f from a to b using Gauss-Legendre
gauss_legendre<-function(f,a,b,nodes,weights) {

	# change of variables from [-1,1] to [a,b]
	ab_nodes <- a + (b-a)*(nodes+1)/2;
	ab_weights <- weights*(b-a)/2;

	# apply Gauss-Legendre rule
	sum <- 0
	for(i in 1:length(nodes)){
        sum <- (sum + ab_weights[i]*f(ab_nodes[i]))}
	return(sum)
}

gauss_legendre(function(y) {
sapply(y, function(y) {
gauss_legendre(function(x) x*y, 0, 1, nodes, weights)$value
})
}, 0, 1, nodes, weights)

Can anybody tell me where the problem lies?
Thank you for your help!

Tiffy