Numerical Integration using chained Gauss-Legendre

[Namir]

Back last year, a member of this website pointed to the simplicity of using the Gauss-Legendre quadrature (with a 3rd order Legendre polynomial) with vintage and new HP calculators. This prospect made me think of using a "chained" version of that type of quadrature to yield relatively good results. Here are my preliminary results usin Excel VBA.

The following Excel VB Code compares the chained Simpson's rule with a "chained" Gauss-Legendre quadrature using a 3rd order Legendre polynomial. The following is the configuration contents of the Excel sheet:

Code:

Cell        Contents
----------------------------------
A1        "A"
A2        "B"
A3        "N"
A4        "FX"
A5        "Simpson"
A6        ""
A7        Exact

B1        value for A, e.g. 1
B2        value for B, e.g. 2
B3        value for N, e.g. 19
B4        string for FX, e.g. "1/X"

C5        Formula for % error for Simpson's rule
C6        Formula for % error for Gauss-Legendre quadrature

Here is the VBA code:

Code:

Function Fx(ByVal sFx As String, ByVal X As Double) As Double
  sFx = Replace(sFx, "EXP(", "!")
  sFx = Replace(sFx, "X", "(" & X & ")")
  sFx = Replace(sFx, "!", "EXP(")
  Fx = Evaluate(sFx)
End Function

Sub go3()
  Const MAX = 3
  Dim A As Double, B As Double, Xa As Double, Xb As Double
  Dim N As Integer, I As Integer, J As Integer, C As Integer, h As Double
  Dim Sum As Double, Sum2 As Double, X As Double
  Dim Xar(MAX) As Double, Wt(MAX) As Double, T1 As Double, T2 As Double
  Dim sFx As String
  
  A = [B1].Value
  B = [B2].Value
  N = [B3].Value
  sFx = [B4].Value
  sFx = UCase(Replace(sFx, " ", ""))
  
  ' Simpson's rule
  If N Mod 2 = 0 Then N = N + 1
  h = (B - A) / (N + 1)
  X = A + h
  Sum = Fx(sFx, A) + Fx(sFx, B) + 4 * Fx(sFx, X)
  C = 2
  For I = 2 To N
    X = X + h
    Sum = Sum + C * Fx(sFx, X)
    C = 6 - C
  Next I
  [B5].Value = h / 3 * Sum
  
  ' Gauss Quadrature
   N = [B8].Value * [B3].Value
   Xar(1) = 0
   Wt(1) = 8 / 9
   Xar(2) = Sqr(3 / 5)
   Wt(2) = 5 / 9
   Xar(3) = -Xar(2)
   Wt(3) = Wt(2)
   
    If N Mod 2 = 0 Then N = N + 1
    h = (B - A) / (N + 1)
    Xa = A
    Xb = Xa + h
    Sum = 0
    Do
      T1 = (Xb - Xa) / 2
      T2 = (Xb + Xa) / 2
      Sum2 = 0
      For J = 1 To MAX
        Sum2 = Sum2 + Wt(J) * Fx(sFx, T1 * Xar(J) + T2)
      Next J
      Sum = Sum + T1 * Sum2
      Xa = Xb
      Xb = Xb + h
    Loop Until Xa >= B
    [A6].Value = "GS" & MAX
    [B6].Value = Sum
End Sub

As you experiment with different functions and integration ranges, you should see that the chained Gauss-Legendre quadrature is significantly more accurate than Simpson's rule. Both methods use three points per divided interval.

Enjoy!

Namir

[Tugdual]

Is this related to hp 41?

[lcwright1964]

(03-22-2016 06:40 AM)Tugdual Wrote:  Is this related to hp 41?

I think that Namir introduced his VBA code with a comment that this approach to quadrature might be well suited to keystroke programming. I discerned at least an implied challenge that someone adapt this to HP41/42 code--unless Namir is working on that himself already.

I have to admit I am more of a Romberg man myself...

Les