"Susanne Pfeifer" <[EMAIL PROTECTED]> wrote in message news:[EMAIL PROTECTED] > 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;[EMAIL PROTECTED]:6-E`"'Y M! $`````+ (`!0`/`2(`@ ````````+_A(^IR^V/[EMAIL PROTECTED]@PUEE MB(JCX3EK"L?;-I[R;0/E&\%\5_O=2$'69XA,BD*YG)*6<)(8+Z$I:CTRB[.G [EMAIL PROTECTED]:!U>AWN7Y-((Y2?)QUMZ_1^<G)[846+%3L]:#10C' MUX;((<6E`VEH*3G(6/<DM**'&<=GB8=#^FC4=ZEHX7JEFCJ:9D/[EMAIL PROTECTED]&_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":,[EMAIL PROTECTED];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;[EMAIL PROTECTED];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 ______________________________________________ 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.
