[R] Double integration - Gauss Quadrature

[Earl F. Glynn]

"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