# Double integration - Gauss Quadrature

3 messages
Open this post in threaded view
|

## Double integration - Gauss Quadrature

 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
 "Susanne Pfeifer" <[hidden email]> wrote in message news:[hidden email]... > Hi, > > I would like to solve a double integral of the form . . . > 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: Maybe there's some way to use sapply as your code suggests, but I'm not sure where you defined the $value that is being returned in your inner call: gauss_legendre(function(x) x*y, 0, 1, nodes, weights)$value I converted some old IDL code to do this 2D integral but without trying to use your sapply: # For N=5, see "known" values here: # http://mathworld.wolfram.com/Legendre-GaussQuadrature.htmllibrary(statmod) N <- 5 GL <- gauss.quad(N) nodes   <- GL$nodes weights <- GL$weights ############################################## # 1D Gauss-Legendre gauss_legendre <- function(f, a, b, nodes, weights) {   C <- (b - a) / 2   D <- (b + a) / 2   sum <- 0.0   for (i in 1:length(nodes))   {     sum <- sum + weights[i] * f(nodes[i]*C + D)   }   return(C * sum) } ############################################## gauss_legendre2D_helper <- function(f, x, a2,b2, nodes, weights) {   C <- (b2 - a2) / 2   D <- (b2 + a2) / 2   sum <- 0.0   for (i in 1:length(nodes))   {     y <- nodes[i]*C + D     sum <- sum + weights[i] * f(x,y)   }   return(C * sum) } gauss_legendre2D <- function(f, a1,b1, a2,b2, nodes, weights) {   C <- (b1 - a1) / 2   D <- (b1 + a1) / 2   sum <- 0.0   for (i in 1:length(nodes))   {     x <- nodes[i]*C + D     sum <- sum + weights[i] * gauss_legendre2D_helper(f, x, a2, b2, nodes, weights)   }   return(C * sum) } ############################################## # 1D Test: gauss_legendre(function(x) {x}, 0.0, 1.0, nodes, weights) # 2D Test: gauss_legendre2D(function(x,y) {x*y}, 0.0, 1.0, 0.0, 1.0, nodes, weights) # Minimal testing here: > # 1D Test: > gauss_legendre(function(x) {x}, 0.0, 1.0, nodes, weights) [1] 0.5 > > # 2D Test: > gauss_legendre2D(function(x,y) {x*y}, 0.0, 1.0, 0.0, 1.0, nodes, weights) [1] 0.25 BTW:  I don't think you need N as large as you're using. The advantage of Gauss-Legendre quadrature is fairly high precision without that many function evaluations. Formulas for those who may be interested: 1D Gauss-Legendre Quadrature 2D Gauss-Legendre Quadrature This can be extended to a 3D integral evaluations, too. efg Earl F Glynn Scientific Programmer Stowers Institute for Medical Research   begin 666 clip_image002.gif M1TE&.#EA% $K'<,2'^&E-O9G1W87)E.B!-:6-R;W-O9G0@3V9F:6-E"'Y M!$+ (!0/2(@ +_A(^IR^V/@Q2T0JSINL^K!U@PUEE MB(JCX3EK"L?;-I[R;0/E&\%\5_O=2$'69XA,BD*YG)*6<)(8+Z$I:CTRB[.G M-Y4U@J3?Z3B\[*AG9R1:!U>AWN7Y-((Y2?)QUMZ_1^ZOO%K^-) MA\6.R9I9AJ3)^:.FBAXZ^.3"RTO&CGZN^K7&__N_G3[^ ;I%Y]O/_]X_/K[4Q7G[U^ 3S0$H8('$G6=@@O9UH^ 3!0.P  end ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
 In reply to this post by Susanne Pfeifer "Susanne Pfeifer" <[hidden email]> wrote in message news:[hidden email]... > 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.