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
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
"Susanne Pfeifer" <tiffy at tiffy.it> wrote in message news:48DE3BE0.5020200 at tiffy.it...> 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.html library(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 at 3V9F:6-E`"'Y M! $`````+ (`!0`/`2(`@ ````````+_A(^IR^V/`@Q2T0JSINL^K!U at PUEE MB(JCX3EK"L?;-I[R;0/E&\%\5_O=2$'69XA,BD*YG)*6<)(8+Z$I:CTRB[.G M-Y4U at J3?Z3B\[*AG9R1:!U>AWN7Y-((Y2?)QUMZ_1^<G)[846+%3L]:#10C' MUX;((<6E`VEH*3G(6/<DM**'&<=GB8=#^FC4=ZEHX7JEFCJ:9D/6YW at Y&_O7 MVCFHUP3\4[6K^E?,B7?1\DJ*N&IZN[J\'!N%_(S:JT(="CV+6BHM&#,I/JK] M*FU,JSMV/?3,+C]M;7_/0Y_-GM:_R.A,5[UUY,"\:X?PSK6#WO*\*RC*6#ID MR?S=*T0KX":,__\BHNLWL6*=@QX%;EQW;J"S8_3*D=R!*9=%?%1>ZOO%K^-) MA\6.R9I9AJ3)</9 *6/I(:FFE&".HIND;YM(E&?,93HD+A&5?]4D0=TE=22\ MK(1:2H7H"^T$.V-]K4WK%E^PJE-YQ;6C5 at M;GW?5Y'71-VC-OUYL)8$8*# 1 MMX05_W0,.;+DR;@:'[9,.;/FS=%07JW&.;3HT:-;DCZ-.G5HTP:[N7X-.[;L MV;1KV[Z-.[?NW;Q[/W3A6K7PX<1ED,%</+GRU<N;.W^NMQ'TZ=1)\U5773FS M[)D->^:.FBAXZ^.3"RTO&CGZN^K7&__N_G3[^ ;I%Y]O/_]X_/K[4Q7G[U^ 3S0$H8('$G6=@@O9UH^ 3!0``.P`` ` end
"Susanne Pfeifer" <tiffy at tiffy.it> wrote in message news:48DE3BE0.5020200 at tiffy.it...> 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