"Susanne Pfeifer" <[EMAIL PROTECTED]> wrote in message
news:[EMAIL PROTECTED]
> 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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
