Hi All,

If you love numerical analysis like I do, then you have most likely come across the various Gaussian quadrature algorithms. For many years I implemented these algorithms by using the roots and weights of a related quadrature polynomial (Legendre, Laguerre, Hermite, etc.) listed in tables for a specific polynomial order. Higher orders gave better quadrature results at the price of using more roots and weights. About two years ago I was able to implement dynamic quadrature function in Matlab. Matlab allowed me to get the roots of the quadrature polynomials very quickly and easily. These routines allowed me to specify the order of the polynomial used and have the routine calculate the roots and weights used to calculate the integral numerically. No more hard-coding tables of roots and weights!

A few days ago I was able to implement similar functions in Excel VBA, even though Excel lacked the ability to find the roots of the quadrature polynomials using a simple function call like in Matlab. The approach I used employed my Scan Range Method that I presented in HHC2012. This method finds the roots of any function within a range of values. I was able to adapt this algorithm to the Gauss-Legendre, Gauss-Laguerre, and Gauss-Hermite quadratures to integrate between [A, B], [0, infinity], and [-infinity, infinity] respectively. The functions I present below use the following parameters to find the polynomial roots:

1. sExpress passes a string that specifies the function to be integrated. For example you can pass the argument "EXP($X)-3*$X^2".

2. VarName is the string that specifies the name of the variable in the parameter sExpress. For example you can pass "$X" or just "X". Avoid "X" if your function string includes "EXP".

3. DeltaX to specify the increments of X used in scanning for teh roots.

4. Toler to specify the tolerance used in calculating the roots.

5. MaxX used to specify when to stop scanning if not all the anticipated roots are found.

In implementing the code in other programming languages, you can access the integrated function using user-defined functions that you declare elsewhere in your code.

I recommend that you use Wikipedia to look up the various Gaussian quadrature methods and their related polynomials.

Here is the Excel VBA code that should be easily readable and adaptable to any other programming language:

Option ExplicitFunction LegendrePoly(ByVal X As Double, ByVal Order As Integer) As Double

' Implementation of recursive relation for Legendre polynomials.

Dim L0 As Double, L1 As Double, L2 As Double

Dim I As Integer, K As IntegerL0 = 1

If Order = 0 Then

LegendrePoly = L0

Exit Function

End IfL1 = X

If Order = 1 Then

LegendrePoly = L1

Exit Function

End If

K = 1

For I = 2 To Order

L2 = ((2 * K + 1) * X * L1 - K * L0) / (K + 1)

K = K + 1

L0 = L1

L1 = L2

Next I

LegendrePoly = L2

End FunctionFunction LegendrePolyDeriv(ByVal X As Double, ByVal Order As Integer) As Double

' First derivative of Legendre polynomials.

Dim h As Doubleh = 0.001 * (1 + Abs(X))

LegendrePolyDeriv = (LegendrePoly(X + h, Order) - LegendrePoly(X - h, Order)) / 2 / h

End FunctionFunction GaussLegendreQuad(ByVal sExpress As String, ByVal sVarName As String, _

ByVal A As Double, ByVal B As Double, _

ByVal Order As Integer, _

Optional ByVal DeltaX As Double = 0.1, _

Optional ByVal Toler As Double = 0.00000001, _

Optional ByVal MaxX As Double = 1000) As DoubleConst MAX_ITER = 1000

Dim Sum As Double

Dim Xa As Double, Fa As Double, Xb As Double, Fb As Double, Fx As Double

Dim Xrt As Double, Wt As Double, Diff As Double, Xval As Double

Dim N As Integer, Iter As IntegerN = Order

Sum = 0

GaussLegendreQuad = 0

Xb = 0

Fb = LegendrePoly(Xb, Order)' found root at x=0 (for odd polynomial orders)?

If Abs(Fb) < 0.000001 Then

N = N - 1

Xrt = Xb

Wt = 2 / (1 - Xrt * Xrt) / (LegendrePolyDeriv(Xrt, Order)) ^ 2

Xval = (B - A) * Xrt / 2 + (A + B) / 2

Fx = Evaluate(Replace(sExpress, sVarName, "(" & CStr(Xval) & ")"))

Sum = Sum + Wt * Fx

End If' find roots of the Legendre polynomial using the scan range method

Do

DoEvents

Xa = Xb

Fa = Fb

Xb = Xa + DeltaX

Fb = LegendrePoly(Xb, Order)

' root in interval [Xa, Xb]?

If Fa * Fb < 0 Then

' start Newton;s method

Xrt = (Xa + Xb) / 2

Iter = 0

Do

DoEvents

Iter = Iter + 1

Diff = LegendrePoly(Xrt, Order) / LegendrePolyDeriv(Xrt, Order)

Xrt = Xrt - Diff

Loop Until Abs(Diff) <= Toler Or Iter > MAX_ITER

N = N - 2 ' decrement 2 ... one for positive-value root and one for it's negative twin' calculate weight for postive and negative roots

Wt = 2 / (1 - Xrt * Xrt) / (LegendrePolyDeriv(Xrt, Order)) ^ 2

' calculate transformed x for positive root

Xval = (B - A) * Xrt / 2 + (A + B) / 2

Fx = Evaluate(Replace(sExpress, sVarName, "(" & CStr(Xval) & ")"))

' update the value of the integral

Sum = Sum + Wt * Fx

' calculate transformed x for negative root

Xval = (A - B) * Xrt / 2 + (A + B) / 2

Fx = Evaluate(Replace(sExpress, sVarName, "(" & CStr(Xval) & ")"))

' update the value of the integral

Sum = Sum + Wt * Fx

End If

Loop Until N = 0 Or Xa > MaxX' found all roots?

If N = 0 Then

GaussLegendreQuad = (B - A) / 2 * Sum

End IfEnd Function

Function LaguerrePoly(ByVal X As Double, ByVal Order As Integer) As Double

' Implementation of recursive relation for Laguerre polynomials.

Dim L0 As Double, L1 As Double, L2 As Double

Dim I As Integer, K As IntegerL0 = 1

If Order = 0 Then

LaguerrePoly = L0

Exit Function

End IfL1 = 1 - X

If Order = 1 Then

LaguerrePoly = L1

Exit Function

End IfK = 1

For I = 2 To Order

L2 = ((2 * K + 1 - X) * L1 - K * L0) / (K + 1)

K = K + 1

L0 = L1

L1 = L2

Next I

LaguerrePoly = L2

End FunctionFunction GaussLaguerreQuad(ByVal sExpress As String, ByVal sVarName As String, _

ByVal Order As Integer, _

Optional ByVal DeltaX As Double = 0.1, _

Optional ByVal Toler As Double = 0.00000001, _

Optional ByVal MaxX As Double = 1000) As Double

' Gauss-Laguerre Quadrature

Const MAX_ITER = 1000Dim Sum As Double

Dim Xa As Double, Fa As Double, Xb As Double, Fb As Double, Fx As Double

Dim Xrt As Double, Wt As Double, h As Double, Diff As Double

Dim N As Integer, Iter As IntegerN = Order

Sum = 0

GaussLaguerreQuad = 0

Xb = 0

Fb = LaguerrePoly(Xb, Order)' find the roots of the Laguerre polynomials using the scan range method

Do

DoEvents

Xa = Xb

Fa = Fb

Xb = Xa + DeltaX

Fb = LaguerrePoly(Xb, Order)

' root in interval [Xa, Xb]?

If Fa * Fb < 0 Then

' start Newton's method

Xrt = (Xa + Xb) / 2

Iter = 0

Do

DoEvents

Iter = Iter + 1

h = 0.001 * (1 + Abs(Xrt))

Diff = 2 * h * LaguerrePoly(Xrt, Order) / (LaguerrePoly(Xrt + h, Order) - LaguerrePoly(Xrt - h, Order))

Xrt = Xrt - Diff

Loop Until Abs(Diff) <= Toler Or Iter > MAX_ITER

N = N - 1

' calculaet weight at x root

Wt = Xrt / (Order + 1) ^ 2 / LaguerrePoly(Xrt, Order + 1) ^ 2

' calculate function value

Fx = Evaluate(Replace(sExpress, sVarName, "(" & CStr(Xrt) & ")"))

' update area value

Sum = Sum + Wt * Fx

End If

Loop Until N = 0 Or Xa > MaxX' foudn all roots?

If N = 0 Then

GaussLaguerreQuad = Sum

End IfEnd Function

Function Factorial(ByVal N As Integer) As Double

Factorial = Application.WorksheetFunction.Fact(N)

End FunctionFunction HermitePoly(ByVal X As Double, ByVal Order As Integer) As Double

' Implementation of recursive relation for Hermite polynomials.

Dim H0 As Double, H1 As Double, H2 As Double

Dim I As Integer, K As IntegerH0 = 1

If Order = 0 Then

HermitePoly = H0

Exit Function

End IfH1 = 2 * X

If Order = 1 Then

HermitePoly = H1

Exit Function

End IfK = 1

For I = 2 To Order

H2 = 2 * X * H1 - 2 * K * H0

K = K + 1

H0 = H1

H1 = H2

Next IHermitePoly = H2

End FunctionFunction GaussHermiteQuad(ByVal sExpress As String, ByVal sVarName As String, _

ByVal Order As Integer, _

Optional ByVal DeltaX As Double = 0.1, _

Optional ByVal Toler As Double = 0.00000001, _

Optional ByVal MaxX As Double = 1000) As DoubleConst MAX_ITER = 1000

Dim Sum As Double

Dim Xa As Double, Fa As Double, Xb As Double, Fb As Double, Fx As Double

Dim Xrt As Double, Wt As Double, Diff As Double, h As Double

Dim N As Integer, Iter As IntegerN = Order

Sum = 0

GaussHermiteQuad = 0

Xb = 0

Fb = HermitePoly(Xb, Order)' found root at x=0 (for odd polynomial orders)?

If Abs(Fb) < 0.000001 Then

N = N - 1

Xrt = Xb

Wt = 2 ^ (Order - 1) * Factorial(Order) * Sqr(4 * Atn(1)) / Order ^ 2 / (HermitePoly(Xrt, Order - 1)) ^ 2

Fx = Evaluate(Replace(sExpress, sVarName, "0"))

Sum = Sum + Wt * Fx

End If' find roots of the Hermite polynomial using the scan range method

Do

DoEvents

Xa = Xb

Fa = Fb

Xb = Xa + DeltaX

Fb = HermitePoly(Xb, Order)

' root in interval [Xa, Xb]?

If Fa * Fb < 0 Then

' start Newton;s method

Xrt = (Xa + Xb) / 2

Iter = 0

Do

DoEvents

Iter = Iter + 1

h = 0.001 * (1 + Abs(Xrt))

Diff = 2 * h * HermitePoly(Xrt, Order) / (HermitePoly(Xrt + h, Order) - HermitePoly(Xrt - h, Order))

Xrt = Xrt - Diff

Loop Until Abs(Diff) <= Toler Or Iter > MAX_ITER

N = N - 2 ' decrement 2 ... one for positive-value root and one for it's negative twin' calculate weight for postive and negative roots

Wt = 2 ^ (Order - 1) * Factorial(Order) * Sqr(4 * Atn(1)) / Order ^ 2 / (HermitePoly(Xrt, Order - 1)) ^ 2

Fx = Evaluate(Replace(sExpress, sVarName, CStr(Xrt)))

' update the value of the integral

Sum = Sum + Wt * Fx

Xrt = -Xrt

Fx = Evaluate(Replace(sExpress, sVarName, "(" & CStr(Xrt) & ")"))

' update the value of the integral

Sum = Sum + Wt * Fx

End If

Loop Until N = 0 Or Xa > MaxX' found all roots?

If N = 0 Then

GaussHermiteQuad = Sum

End IfEnd Function

Enjoy!

*Edited: 29 July 2013, 3:44 p.m. *