Skip to content

Double integration - Gauss Quadrature

3 messages · Susanne Pfeifer, Earl F. Glynn

Susanne Pfeifer
# 
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
2 days later
Earl F. Glynn
# 
"Susanne Pfeifer" <tiffy at tiffy.it> wrote in message 
news:48DE3BE0.5020200 at tiffy.it...
. . .
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:
[1] 0.5
[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
Earl F. Glynn
# 
"Susanne Pfeifer" <tiffy at tiffy.it> wrote in message 
news:48DE3BE0.5020200 at tiffy.it...
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