# Numerical integration with Gaussian Quadrature

#pragma rtGlobals=3     // Use modern global access method and strict wave access.

//Perform Gaussian Quadrature Integration of a given function.
//This is slightly different to the inbuilt Integrate1D in that one can pass in a wave containing wave references as extra input
//to the function to be integrated!
//The function to be integrated can also be threadsafe to speed up calculation.

//There is also a threadsafe version of the integration code (integrate1D is not threadsafe). In this case the function to be integrated also needs to be threadsafe.

//Lookup tables are used to store commonly used values of the Gauss-Legendre nodes.
//The GaussLegendre node calculation is based on code from the NIST SANS macros, which was originally based on code from
//Numerical recipes.
//Type test() for a quick test of the code.
//Written by Andrew Nelson July 2013

Function GQIntegrate1D(fcn, min_x, max_x, w, N, [tol, rtol, adaptive, startiter])
//a function for performing Gaussian Quadrature integration with a set order of quadrature points
//INPUT
//fcn           - the function to be integrated.  Has the signature:
//                      Function fcn(w, y, x)
//                          wave/wave w
//                          wave y, x
//                          y = f(x)
//                      End
//                  If the function is declared threadsafe, then the integration is done in a parallel manner.   The function is written
//                  in much the same manner as an all-at-once fit function.  i.e. All the abscissa (xpoints) are passed in at once, and the function
//                  is supposed to populate all the y values at once.
//min_x     - the lower limit for the integration
//max_x     - the upper limit for the integration
//w         - wave containing wave references to pass into the function fcn.  Use this to pass data into your function
//N             - number of quadrature points to be used. If the adaptive variable is set then this is the maximum order examined.
//
//tol           - the Iteration stops when the error between successive iterates is less than tol OR the relative change is less than rtol. default = 1e-13
//rtol          - the Iteration stops when the error between successive iterates is less than tol OR the relative change is less than rtol.

string fcn
Variable min_x, max_x
Wave w
Variable N, tol, rtol, adaptive, startiter

// local variables
variable ii, jj, va, vb, isthreadsafe, NORD, err, val, maxiter
string info = "", funcname = "", fri = ""

if(paramisdefault(tol))
tol = 1.e-13
endif
if(paramisdefault(rtol))
rtol = 1.49e-08
endif
if(paramisdefault(startiter))
startiter = 5
endif

maxiter = N
startiter = N
make/free/n=(maxiter - startiter + 1)/d W_gaussQuad
else
maxiter = N
make/o/n=(maxiter - startiter + 1)/d W_gaussQuad
endif

funcname = fcn
info = functioninfo(funcname)

//limits of integration are input to function
va = min_x
vb = max_x

for(NORD = startiter, jj = 0 ; NORD <= maxiter ; NORD += 1,  jj += 1)
//calculate the Gauss Legendre abscissae and weights
Wave GLpoints = gauss_legendre_tbl(NORD)

//place to put evaluated data  and x point to be evaluated
make/n=(NORD)/d/free yy, zi
yy = 0

// calculate Gauss points on integration interval (q-value for evaluation)
Multithread zi[] = ( GLpoints[p][0] * (vb - va) + vb + va ) / 2.0

//evaluate function at all the points
fri = funcrefinfo(cfcn)
if(numberbykey("ISPROTO", fri) == 0)
cfcn(w, yy, zi)
else
print "Function does not satisfy the signature:"
return NaN
endif
else
fri = funcrefinfo(cfcn1)
if(numberbykey("ISPROTO", fri) == 0)
cfcn1(w, yy, zi)
else
print "Function does not satisfy the signature:"
return NaN
endif
endif

//multiply by weights to calculate partial sum
W_gaussQuad[jj] = (vb - va) / 2.0 * sum(yy)
if(jj)
err = abs(val - W_gaussQuad[jj - 1])
if(err < tol || err < rtol * abs(val))
break
endif
endif
endfor

//calculate value to return
return  val
End

//INPUT
//fcn           - the function to be integrated.  Has the signature:
//                      Threadsafe Function fcn(w, y, x)
//                          wave/wave w
//                          wave y, x
//                          y = f(x)
//                      End
//min_x         - the lower limit for the integration
//max_x     - the upper limit for the integration
//w         - wave containing wave references to pass into the function fcn.  Use this to pass data into your function
//N             - number of quadrature points to be used. If the adaptive variable is set then this is the maximum order examined.
//
//tol           - the Iteration stops when the error between successive iterates is less than tol OR the relative change is less than rtol. default = 1e-13
//rtol          - the Iteration stops when the error between successive iterates is less than tol OR the relative change is less than rtol.

string fcn
Variable min_x, max_x
Wave w
Variable N, tol, rtol, adaptive, startiter

// local variables
variable ii, jj, va, vb, isthreadsafe, NORD, err, val, maxiter
string info = "", funcname = "", fri = ""

if(paramisdefault(tol))
tol = 1.e-13
endif
if(paramisdefault(rtol))
rtol = 1.49e-08
endif
if(paramisdefault(startiter))
startiter = 5
endif

maxiter = N
startiter = N
make/free/n=(maxiter - startiter + 1)/d W_gaussQuad
else
maxiter = N
make/o/n=(maxiter - startiter + 1)/d W_gaussQuad
endif

//limits of integration are input to function
va = min_x
vb = max_x

for(NORD = startiter, jj = 0 ; NORD <= maxiter ; NORD += 1,  jj += 1)
//calculate the Gauss Legendre abscissae and weights
Wave GLpoints = gauss_legendre(NORD)

//place to put evaluated data  and x point to be evaluated
make/n=(NORD)/d/free yy, zi

// calculate Gauss points on integration interval (q-value for evaluation)
Multithread zi[] = ( GLpoints[p][0] * (vb - va) + vb + va ) / 2.0

//evaluate function at each point
fri = funcrefinfo(cfcn)
if(numberbykey("ISPROTO", fri) == 0)
cfcn(w, yy, zi)
else
print "Function does not satisfy the signature:"
return NaN
endif

//multiply by weights to calculate partial sum
W_gaussQuad[jj] = (vb - va) / 2.0 * sum(yy)

if(jj)
err = abs(val - W_gaussQuad[jj - 1])
if(err < tol || err < rtol * abs(val))
break
endif
endif
endfor

//calculate value to return
return  val
End

Wave/wave w
Wave y, x
return Inf
end

Wave/wave w
wave y, x
return Inf
end

Function/wave gauss_legendre_tbl(N)
variable N
///routine from Numerical Recipes to calculate weights and abscissae for
// Gauss-Legendre quadrature, over the range (-1, 1)
//
//Use this function in preference to gauss_legendre. This function creates a lookup table and is faster
//to use in the longrun, compared to gauss_legendre which calculates new values each time
//INPUT - order for Gauss Legendre (>> 1)
//
//RETURNS - free wave containing Gauss abscissa (first column) and weights (second column).
//          The wave has dimensions [N][2]

//  variable timer = startmstimer
variable newsize, original_idx_size

DFREF saveDFR = GetDataFolderDFR()
newdatafolder/o root:packages
make/n=0/i/u/o idx
make/n=(0, 0)/d/o wt = 0, zi = 0
SetDataFolder saveDFR
endif

findvalue/I=(N)/z idx
if(V_Value == -1)
//haven't created those values yet
original_idx_size = numpnts(idx)
if(N > dimsize(wt, 0))
newsize = N
else
newsize = dimsize(wt, 0)
endif
redimension/n=(original_idx_size + 1) idx
redimension/n=(newsize, original_idx_size + 1) zi, wt
Wave vals = gauss_legendre(N)
idx[numpnts(idx) - 1] = N
zi[0, N - 1][numpnts(idx) - 1] = vals[p][0]
wt[0, N - 1][numpnts(idx) - 1] = vals[p][1]
//      print stopmstimer(timer)/1e6
return vals
else
make/n=(N, 2)/d/free GLpoints
GLpoints[][0] = zi[p][V_Value]
GLpoints[][1] = wt[p][V_Value]
//      print stopmstimer(timer)/1e6
return GLpoints
endif
End

///routine from Numerical Recipes to calculate weights and abscissae for
//Code modified from GaussUtils_v40.ipf, written by Steve Kline @NCNR NIST
//
//INPUT - order for Gauss Legendre (>> 1)
//
//RETURNS - free wave containing Gauss abscissa (first column) and weights (second column).
//          The wave has dimensions [N][2]

variable N

Variable x1, x2
variable m, jj, ii
variable z1, z, xm, xl, pp, p3, p2, p1
Variable eps = 3e-11

make/free/n=(N, 2)/d GLpoints

x1 = -1
x2 = 1

m = (N+1)/2
xm = 0.5 * (x2 + x1)
xl = 0.5 * (x2 - x1)

for (ii = 1; ii <= m; ii += 1)
z = cos(pi * (ii - 0.25)/(N + 0.5))
do
p1 = 1.0
p2 = 0.0
for (jj = 1;jj <= N;jj += 1)
p3 = p2
p2 = p1
p1 = ((2.0 * jj - 1.0) * z * p2 - (jj - 1.0) * p3) / jj
endfor
pp = N * (z * p1 - p2) / (z * z -1.0)
z1 = z
z = z1 - p1 / pp
while (abs(z - z1) > EPS)
GLpoints[ii - 1][0] = xm - xl * z
GLpoints[N  - ii][0] = xm + xl * z
GLpoints[ii - 1][1] = 2.0 * xl / ((1.0 - z * z) * pp * pp)
GLpoints[N  - ii][1] = GLpoints[ii - 1][1]
Endfor

return GLpoints
End

//Lets integrate x^2
Wave/wave w
wave y, x
y = x^2
End

Function test_sine(w, y, x)
//Lets integrate x^2
Wave/wave w
wave y, x
y = sin(1/x)
End

Function test_sine_IGOR(x)
variable x
return sin(1/x)
End

Function test()
variable x1, x2, true, err
Wave/wave w
x1 = 0
x2 = 1
print "Integrating x^2 between", x1, "and", x2
print/d GQIntegrate1D("test_quadratic", x1, x2, w, 100)
print/d GQIntegrate1D("test_quadratic", x1, x2, w, 100, adaptive = 1, startiter = 3)

//sin(x)
x1 = 0.001
x2 = pi / 2
print "Integrating sin(1/x) between", x1, "and", x2
print "IGOR's version"
true = integrate1D(test_sine_igor, 0.001, pi/2, 2, 0)
print/d true
err = GQIntegrate1D("test_sine", x1, x2, w, 100)
print/d true - err
err =  GQIntegrate1D("test_sine", x1, x2, w, 1000, adaptive = 1, startiter = 2)
print true - err
End

Forum

Support

Gallery