Dynamic Gaussian Quadrature code in Excel VBA « Next Oldest | Next Newest »

 ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 07-29-2013, 08:14 AM 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 Explicit Function 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 Integer L0 = 1 If Order = 0 Then LegendrePoly = L0 Exit Function End If L1 = 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 Function Function LegendrePolyDeriv(ByVal X As Double, ByVal Order As Integer) As Double ' First derivative of Legendre polynomials. Dim h As Double h = 0.001 * (1 + Abs(X)) LegendrePolyDeriv = (LegendrePoly(X + h, Order) - LegendrePoly(X - h, Order)) / 2 / h End Function Function 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 Double Const 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 Integer N = 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 If End 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 Integer L0 = 1 If Order = 0 Then LaguerrePoly = L0 Exit Function End If L1 = 1 - X If Order = 1 Then LaguerrePoly = 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 LaguerrePoly = L2 End Function Function 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 = 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, h As Double, Diff As Double Dim N As Integer, Iter As Integer N = 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 If End Function Function Factorial(ByVal N As Integer) As Double Factorial = Application.WorksheetFunction.Fact(N) End Function Function 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 Integer H0 = 1 If Order = 0 Then HermitePoly = H0 Exit Function End If H1 = 2 * X If Order = 1 Then HermitePoly = H1 Exit Function End If K = 1 For I = 2 To Order H2 = 2 * X * H1 - 2 * K * H0 K = K + 1 H0 = H1 H1 = H2 Next I HermitePoly = H2 End Function Function 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 Double Const 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 Integer N = 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 If End Function ``` Enjoy! Edited: 29 July 2013, 3:44 p.m. ▼ Kimberly Thompson Member Posts: 97 Threads: 1 Joined: Jun 2013 07-30-2013, 09:50 AM Namir Very interesting; thanks for the submittal! SlideRule ps: I'm s-l-o-w-l-y working my way thru, thanks again. ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 07-30-2013, 11:20 AM I guess one can first translate the code to HP-71B BASIC. I have no illusions that an HP-41C/42 version will run slow. In the case of calculators I suggest having the part that calculators the roots of the polynomial and their associated weights separate. This scheme allows you to reuse the roots/weights for the same polynomial order for different integral ranges (this applies to Gauss-Legendre quadrature for the ranges of [A, B]). The Gaussian quadrature is an interesting numerical analysis algorithm that requires the calculations of ALL of a polynomial roots. This is where a method like my Scan Range method comes in handy. Namir Edited: 30 July 2013, 3:09 p.m. ▼ Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 07-30-2013, 05:15 PM The 34S should be able to do this relatively well. The orthogonal polynomials are built ins, although actually implemented as key stroke programs. - Pauli ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 07-30-2013, 07:37 PM It will be on my list of machines to use!!! :-) Namir

 Possibly Related Threads... Thread Author Replies Views Last Post HP PRIME: APP program code DISAPPEARS !! Joseph Ec 0 691 11-25-2013, 11:35 AM Last Post: Joseph Ec Where to the 32-bit version of User Code Utiltiy for HP-41 ? Olivier (Wa) 2 1,099 09-26-2013, 01:55 AM Last Post: Olivier (Wa) A HP42S Code Editor Andreas 9 2,058 09-22-2013, 03:17 AM Last Post: Andreas HP-65 Morse Code Dan Lewis 7 1,766 01-29-2013, 05:22 PM Last Post: Mike T. OT: Excel regression modeling using VBA Namir 0 628 12-06-2012, 10:05 AM Last Post: Namir [WP34s] undefined OP-code fhub 21 3,440 04-28-2012, 04:09 AM Last Post: Marcus von Cube, Germany OT: To All Number Crunchers: Don't miss this VBA numeric library Namir 0 568 04-04-2012, 11:25 AM Last Post: Namir 32-bit version of User Code Utiltiy for HP-41 MichaelG 4 1,234 02-08-2012, 07:12 AM Last Post: MichaelG Interest in 32-bit version of User Code Utiltiy for HP-41 MichaelG 16 2,903 01-30-2012, 02:01 PM Last Post: MichaelG HP-41 printer interface question. (To the M-code gurus out there) Diego Diaz 2 1,018 01-29-2012, 12:18 PM Last Post: Diego Diaz

Forum Jump: