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