HP Forums

Full Version: Numerical Integration using chained Gauss-Legendre
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
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

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
Is this related to hp 41?
(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
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :