HP 71B and HP Prime Python: Gaussian Quadrature

HP 71B and HP Prime Python: Gaussian Quadrature

A Very Brief Introduction

This is a top-level overview of the Gaussian Quadrature method. For more details, please check out the Sources section.

The design of GQ is to calculate:

∫( f(x) dx, -1, 1) ≈ Σ( wi * f(xi), i = 1, n)

Where do xi and wi come from?

* The xi values come from find the roots of the polynomials P’(x)= 0

* Pn’(x) are the derivatives of the Legendre polynomials Pn(x)

* After getting the roots, xi, the weights, w are calculated by the formula:

wi = 2 / ((1 -xi)^2 * Pn’(xi)^2)

Weights and Points of orders 3 and 5:

Among 3 points:

Roots (x) come from the derivative of the Legendre polynomial of the 3^rd order.

Pn’(x) = 3/2 * (5*x^2 – 1)

Point xi	Weight wi
-√(3/5) ≈ -0.7745 9669 2 = -√0.6	5/9
0	8/9
√(3/5) ≈ 0.7745 9669 2 = √0.6	5/9

Among 5 points:

Roots (x) come from the derivative of the Legendre polynomial of the 5^th order.

Pn’(x) = 15/8 * (21*x^4 – 14*x^2 + 1)

Point xi	Weight wi
-0.90617 98459	0.23692 68851
-0.53846 93101	0.47862 86705
0	128/225 ≈ 0.56888 88889
0.53846 93101	0.47862 86705
0.90617 98459	0.23692 68851

We can fit any interval [a, b] by this conversion:

∫( f(x) dx, a, b) = (b – a) / 2 * ∫( f( (b – a) / 2* x + (b + a) / 2 dx, -1, 1)

= (b – a) / 2 * Σ( wi * f( (b – a) / 2 * xi + (b + a) / 2 ), i = 1, n)

However, the programs on this blog entry will focus on the basic version of the Gaussian Quadrature.

The Code

The following code works with the following integrals:

∫( f(x) dx, -1, 1) ≈ Σ( wi * f(xi), i = 1, n)

HP 71B Basic

HP 71B: GQ3 (Gaussian Quadrature – 3 Point: Integrals from -1 to 1)

10 DEF FNF(X) = insert function of X here

15 S = 0

20 RADIANS

25 FOR I = 1 TO 3

30 READ X, W

31 DATA -SQR(0.6),5/9

32 DATA 0,8/9

33 DATA SQR(0.6),5/9

40 S = S + W * FNF(X)

55 NEXT I

50 DISP “ INTEGRAL= ”; S

HP 71B: GQ5 (Gaussian Quadrature – 5 Point: Integrals from -1 to 1)

10 DEF FNF(X) = insert function of X here

15 S = 0

20 RADIANS

25 FOR I = 1 TO 5

30 READ X, W

31 DATA -.9061798459, .2369268851

32 DATA -.5384693101, .4786286705

33 DATA 0, 128/255

34 DATA .5384693101, .4786286705

35 DATA .9061798459, .2369268851

40 S = S + W * FNF(X)

45 NEXT I

50 DISP “ INTEGRAL= ”; S

Integral Test: ∫( f(x) dx, -1, 1) ≈ Σ( wi * f(xi), i = 1, n)

Function f(x)	Actual (approx)	Gaussian 3 point (GQ3)	Gaussian 5 point (GQ5)
(x-3)*(x+4)*(x-5)	352/3 ≈ 117.333333333	117.333333333	117.333333334
-x^2 + 9*x + 10	58/3 ≈ 19.33333333	19.3333333334	19.3333333333
1.5 * cos(x)	2.524412954	2.52450532159	2.52441295561
exp(-x)	2.350402387	2.35033692868	2.35040238646
exp(sin(x))	2.283194521	2.28303911433	2.28319665353
sin(x)	0	0	0
1/(x - 2)	-1.09861228867	-1.09803921569	-1.09860924181

HP Prime Python App

HP Prime Python: gq3.py

# Gaussian Quadrature 3 point

from math import *

def f(x):

  # insert Function here

  return 1/(x-2)

s=0

x=[-sqrt(0.6),0,sqrt(0.6)]

w=[5/9,8/9,5/9]

for i in range(3):

  s=s+w[i]*f(x[i])

print("Integral: ",s)

HP Prime Python: gq5.py

# Gaussian Quadrature 5 point

from math import *

def f(x):

  # insert Function here

  return sin(x)

s=0

x=[0,-.5384693101056831,.5384693101056831,-.9061798459386640, .9061798459386640]

w=[.5688888888889,.4786286704993665,.4786286704993665,.2369268850561891,.2369268850561891]

for i in range(5):

  s=s+w[i]*f(x[i])

print("Integral: ",s)

Function f(x)	Actual (approx)	Gaussian 3 point (GQ3)	Gaussian 5 point (GQ5)
(x-3)*(x+4)*(x-5)	352/3 ≈ 117.333333333	117.333333333	117.333333333
-x^2 + 9*x + 10	58/3 ≈ 19.33333333	19.3333333333	19.3333333334
1.5 * cos(x)	2.524412954	2.52450532159038	2.52441295561081
exp(-x)	2.350402387	2.35033692868001	2.35040238646284
exp(sin(x))	2.283194521	2.2830391143268	2.28319665352905
sin(x)	0	0.0	0.0
1/(x - 2)	-1.09861228867	-1.09803921568627	-1.09860924181248

Sources

Kamermans, Mike “Pomax”. “Gaussian Quadrature Weights and Abscissae” https://pomax.github.io/bezierinfo/legendre-gauss.html 2011. Retrieved January 19, 2025.

Wikipedia. “Gauss-Legendre quadrature” https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature Last edited January 19, 2025. Retrieved January 19, 2025.

Wikipedia. “Legendre Polynomials” https://en.wikipedia.org/wiki/Legendre_polynomials Last edited December 4, 2024. Retrieved January 19, 2025.

When it comes to evaluating integrals numerically, is Gaussian Quadrature better than Simpson’s Method?

Eddie

 All original content copyright, © 2011-2025. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

The content on this blog is 100% generated by humans. The author does not use AI engines and never will.

Casio fx-CG 50 and Python: Gaussian Quadrature

 Casio fx-CG 50 and Python:  Gaussian Quadrature

The Method

Gaussian quadrature is a method of numerical integration of a function f(x) over an interval [a, b] where a < b. Quadrature works similar to the Trapezoid Rule or Simpson's Rule, but instead of equally spaced points of measure, points are determined by the roots of a certain polynomial. Polynomials used include the Legendre, Jacobi, Chebyshev, Hermite, and Laguerre. Each point, xi, will be assigned a corresponding certain weight, wi, designed to make the approximation the most accurate as possible.

In general, ∫( f(x) dx, x = a to x = b) ≈ Σ( wi * f(xi) from i = 1 to n)

For today’s blog, the Legendre polynomial is used and the domain of integration is restricted to the interval [ -1, 1 ]. Hence:

∫( f(x) dx, x = -1 to x = 1) ≈ Σ( wi * f(xi) from i = 1 to n)

Steps:

1. Determine the number of points (order) is needed.

2. Find the roots of the Legendre polynomial, P_n(x) = 0. This determines the points xi.

3. Calculate the weights as follows:

wi = 2 / ((1 – xi^2) * (P’_n(xi))^2)

where P’_n(xi) is the numerical derivative of the Legendre polynomial at point x = xi.

4. Approximate ∫( f(x) dx, x = -1 to x = 1) ≈ Σ( wi * f(xi) from i = 1 to n).

Legendre Polynomials

There are several ways to determine the Legendre polynomials of order n, these are just a few:

P_n(x) = 1 / (2^n * n!) d^n/dx^n (x^2 – 1)^n

P_n(x) = Σ( comb(n, k), * comb(n + k, k) * comb(x -1, 2)^k for k = 0 to n)

where comb is the combination method

Recursive method:

(n + 1) * P_n(x) = (2*n + 1) * x * P_n-1(x) – n * P_n-2(x)

where P_0(x) = 1 and P_1(x) = x.

A table of Legendre polynomials of orders 0 to 10 is found here (scroll down a bit):

https://en.wikipedia.org/wiki/Legendre_polynomials

An Example: Order 3

Let’s determine the points xi and the weights wi for order 3 (n = 3).

The Legendre polynomial of order 3 is:

P_3(x) = 1 / 2 * (5 * x^3 – 3 * x) = 5 / 2 * x^3 – 3 / 2 * x

And it’s derivative:

P’_3(x) = 15 / 2 * x^2 – 3 / 2

Solve P_3(x) = 0 yields:

5 / 2 * x^3 – 3 / 2 * x = 0

x * (5 / 2 * x^2 – 3 / 2) = 0

x= 0 or 5 / 2 * x^2 – 3 / 2 = 0:

5 / 2 * x^2 – 3 / 2 = 0

5 / 2 * x^2 = 3 / 2

x^2 = 3 / 5

x = ±√(3 / 5)

The roots are, from lowest value to highest value are: x1 = - √(3 / 5), x2 = 0, x3 = √(3 / 5)

Determining the associated weights yields: w1 = 5 / 9 , w2 = 8 / 9, w3 = 5 / 9

Our points and weights are:

x1 = -√(3 / 5), w1 = 5 / 9

x2 = 0, w2 = 8 / 9

x3 = √(3 / 5), w3 = 5 / 9

Casio fx-CG 50 Program: QUADRAT0

(The last character in the file name is a zero)

“3PT QUADRATURE”

“-1 TO 1”

“NO QUOTES”

“F(X)”? → fn1

∫(fn1, A, B) → D

{-√(3 _| 5), 0, √(3 _| 5) } → List 1

{ 5 _| 9, 8 _| 9, 5 _| 9 } → List 2

0 → S

For 1 → I To 3

List 1[ I ] → X

List 2[ I ] → W

S + W × fn1 → S

Next

ClrText

Black Locate 1, 3, “∫(fn1 DX)”

Blue Locate 1, 4, “D=”

Red Locate 4, 4, D

Blue Locate 1, 5, “S=”

Red Locate 4, 5, S

_| is the fraction symbol, pressed by the [ [] / [] ] key.

If you have a monochrome calculator, such as the fx-9750G series or the fx-9860G series, leave out the color commands (Black, Blue, Red) in front of the Locate commands.

The function is stored in function memory, slot 1, which is found in the OPTN menu.

Variables:

D = actual integral, using the integral command of the calculators

S = approximation determined the quadrature

Examples

Radians angle assumed	D (actual)	S (approximate)
∫( sin x dx, x = -1 to x = 1)	0	0
∫( e^(.5 * x) dx, x = -1 to x = 1)	2.084381222	2.084380222
∫( ln(x + 4) / 2 dx, x = -1 to x – 1)	1.365058335	1.365060354

7^th Order and Expanding the Range

The 7^th order Legendre polynomial is:

P_7 = 1 / 16 * (429 * x^7 – 693 * x^5 + 315 * x^3 – 35 * x)

The associated roots (xi) and weights (wi) are:

Point, xi =	Weight, wi =
-0.949107912342759	0.129484966168870
-0.741531185599394	0.279705391489277
-0.405845151377397	0.381830050505119
0	0.417959183673469
0.405845151377397	0.381830050505119
0.741531185599394	0.279705391489277
0.949107912342759	0.129484966168870

These values are from the table presented in “Gauss-Kronrod quadrature formula” (see Source). I verified the roots with Wolfram Alpha, which you can see here:

https://www.wolframalpha.com/input?i2d=true&i=Divide%5B1%2C16%5D+*++%5C%2840%29429+*+Power%5Bx%2C7%5D+%E2%80%93+693+*+Power%5Bx%2C5%5D+%2B+315+*+Power%5Bx%2C3%5D+%E2%80%93+35+*+x%5C%2841%29%3D0

We are not limited to the interval [-1, 1]. To use any interval [a, b], we can scale the points as so:

xi ‘ = (b – a) / 2 * xi + (b + a) / 2

wi’ = (b – a) / 2 * wi

Casio fx-CG50 Program: GAUSS7

This the 7-point Gaussian quadrature, the calculator basic version.

“7PT GAUSS RULE”

“NO QUOTES”

“F(X)”? → fn1

“A”? → A

“B”? → B

∫(fn1, A, B) → D

{ -0.949107912342759,-0.741531185599394,-0.405845151377397,0,0.405845151377397,0.741531185599394,0.949107912342759 } → List 1

{ 0.12948496618870,0.279705391489277,0.381830050505119,0.417959183673469,0.381830050505119,0.279705391489277,0.129484966168870 } → List 2

0 → S

For 1 → I To 7

(B – A) ÷ 2 × List 1[ I ] + (B + A) ÷ 2 → X

List 2[ I ] × (B – A) ÷ 2 × fn1 + S → S

Next

ClrText

Black Locate 1, 3, “∫(fn1 DX)”

Blue Locate 1, 4, “D=”

Red Locate 4, 4, D

Blue Locate 1, 5, “S=”

Red Locate 4, 5, S

Python Code: gauss7.py

The only module needed is the math module, so it should work with all calculators. This script was created with the fx-CG 50.

from math import *

# 7 point Gauss Integral Approximation

def f(x):

return (x-1/2)*(x+1/4)

# xi and wi

xi=[-0.949107912342759,-0.741531185599394,-0.405845151377397,0,0.405845151377397,0.741531185599394,0.949107912342759]

wi=[0.12948496618870,0.279705391489277,0.381830050505119,0.417959183673469,0.381830050505119,0.279705391489277,0.129484966168870]

# sum

s=0

# limits

a=eval(input("a? "))

b=eval(input("b? "))

# integral approximation

for i in range(7):

  x=(b-a)/2*xi[i]+(b+a)/2

  w=(b-a)/2*wi[i]

  s+=w*f(x)

print("Approx. Integral:\n",s)

Examples

Radians mode	D (actual)	S (Casio Basic)	S (Python)
∫( 1 / √(x + 3) dx, 1, 5)	1.656854249	1.658042205	1.656854249504868
∫( .03 * e^(-x^2) dx, -1, 1)	0.04480944797	0.04484387349	0.04480944866632411
∫( e^(-x) dx, -5, 0)	147.4131591	147.5268171	147.4131590824608
∫( .3 * ( cos x – 0.05) dx, 0, π / 2)	0.2764380551	0.2767068146	0.2764380551025117
∫( (x – 1 / 2) * (x + 1 / 4) dx, -1, 1)	0.4166666667	0.4168576653	0.41666666666867568

Notice a difference between the Basic and Python?

Observations

* The program calls for approximations of xi and wi that are carried out 15 digits. Since the roots of the Legendre polynomials are most likely to be irrational, the more digits used, the better (generally).

* Most calculators will store up to 10 to 13 digits internally.

* By contrast, Python will store many more decimal places, at least 15 to 17.

* The approximation uses a lot of calculations. Round off errors can affect the final results, even when the round off is small.

* I recommend the Quadrature method if you are using Python or any other platform that allows for a lot of decimal places.

* Thankfully, most calculators have the numeric integration function.

Sources

“Guass-Kronrod quadrature formula” Wikipedia. Last Edited December 27, 2023. Accessed July 3, 2024. https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula

“Gaussian quadrature” Wikipedia. Last Edited July 7, 2024. Accessed July 7, 2024. https://en.wikipedia.org/wiki/Gaussian_quadrature

This is a long one, folks. Thank you and have a great day. To all the students, may your year of studies be full of success, discovery, and joy.

Until next time,

Eddie

All original content copyright, © 2011-2024. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.