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:
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.
