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 |
7th Order and Expanding the Range
The 7th 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:
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.
