# [R] Double integration - Gauss Quadrature

Susanne Pfeifer tiffy at tiffy.it
Sat Sep 27 16:03:36 CEST 2008

Hi,

I would like to solve a double integral of the form

\int_0^1 \int_0^1 x*y dx dy

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?