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.

Measuring a Distance with an Object in the Way

 Measuring a Distance with an Object in the Way

A Surveying Problem

Refer to the picture above. Here’s the situation: We are given a line and we are to determine a distance, X, from the Observer to the Survey Point. There is one problem though, there is a tree in the way, and the distance we are to measure goes directly through the tree.

Perhaps we can estimate the distance. Pick a point where the line of sight is not blocked. Let this new point be the Auxiliary Point and let C be the distance between the Observer and the Auxiliary Point.

We can also measure the distance between the Auxiliary Point and Survey Point, name this distance H. We can use the “right triangle” created to find the original distance X.

By the Pythagorean Theorem, we have:

X^2 + H^2 = C^2

X^2 = C^2 – H^2

X = √(C^2 – H^2) (only positive roots make geometric sense)

What if the calculator device does not have the square root function? Or if you are working by hand? In the 1938 book, The Principles and Practice of Surveying, Volume I (see the Source section), we can estimate the length of X without the need of taking a square root. Here is the derivation:

X^2 + H^2 = C^2

H^2 = C^2 – X^2

Observe that for any t, s: t^2 – s^2 = (t – s) × (t + s), and:

H^2 = (C – X) × (C + X)

If H is small, then C is close to X. Assume that H is small. As a result, C ≈ X. This derivation applies this approximation only to the term C + X:

H^2 ≈ (C – X) × (C + C)

Solve for X:

H^2 ≈ (C – X) × 2 × C

H^2 ÷ (2 × C) ≈ C – X

X ≈ C - H^2 ÷ (2 × C)

How good is the estimated formula?

Testing the Approximation

I used the HP 15C Collector’s Edition for testing the approximation formula. The code used to compare the approximate to the actual result:

Line #; Key Code; Key

Stack: y: hypotenuse (C), x: side length (H); (y > x)

01; 42, 21, 11; LBL A // approximate calculation

02; __, 44, _1; STO 1

03; __, __, 34; x<>y

04; __, 44, _2; STO 2

05; __, 45, _2; RCL 2

06; __, 45, _1; RCL 1

07; __, 43, 11; x^2

08; __, __, _2; 2

09; __, __, 10; ÷

10; __, 45, _2; RCL 2

11; __, __, 10; ÷

12; __, __, 30; -

13; __, 44, _3; STO 3

14; __, __, 31; R/S

15; __, 45, _2; RCL 2 // actual calculation

16; __, 43, 11; x^2

17; __, 45, _1; RCL 1

18; __, 43, 11; x^2

19; __, __, 30; -

20; __, __, 11; √

21; __, 44, _4; STO 4

22; __, __, 31; R/S

23; __, 45, _3; RCL 3 // absolute error

24; __, 45, _4; RCL 4

25; __, __, 30; -

26; __, 43, 16; ABS

27; __, 43, 32; RTN

(27 steps, 34 bytes)

Sample Data (rounded to 7 digits)

C	H	Approx.	Actual	Abs. Error
5	1	4.9	4.8989795	0.0010205
6	1	5.9166667	5.9160798	0.0005869
7	1	6.2985174	6.9282032	0.0003682
8	1	7.9375	7.9372539	0.0002461
9	1	8.9444444	8.9442719	0.0001725
10	1	9.95	9.9498744	0.0001256
5	2	4.6	4.5825757	0.0174243
6	2	5.6666667	5.6556842	0.0098124
7	2	6.7142857	6.7082039	0.0060818
8	2	7.75	7.7459667	0.0040333
9	2	8.7777778	8.7749644	0.0028134
10	2	9.8	9.7979590	0.0020410
5	3	4.1	4	0.1
6	3	5.25	5.1961524	0.0538476
7	3	6.3571429	6.3245553	0.0325875
8	3	7.4375	7.4161985	0.0213015
9	3	8.5	8.4852814	0.0147186
10	3	9.55	9.5393920	0.0106080
20	19	10.975	6.2449980	4.7300020
20	1	19.975	19.9749844	0.0000156
50	49	25.99	9.9498744	16.0401256
50	1	49.99	49.9899990	0.0000010

From the data above, it seems that the greater the difference between C and H, in general, the better the approximation.

Source

Breed, Charles B. and George L. Hosmer. The Principles and Practice of Surveying: Volume I. Elementary Surveying. 7^th Edition. John Wiley & sons, Inc. London, Chapman & Hall, Limited. 1938. pg. 14

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.

HP Prime: Interval Arithmetic

HP Prime:  Interval Arithmetic

Interval Arithmetic

In computing, we are working with approximate calculation.  Often a number can be stated as a small interval to account for round off error.  

An interval is a range of two numbers from a to b:  [a, b].  The variable a represents the lower bound, with the variable b represents the upper bound.  There is no restriction on how large or small the interval is.

For example, the interval [ 81.8, 82.2 ] represents the range from 81.8 to 82.2.  This could stand for a calculated value of 82.0 with a plus-or-minus error of 0.2.  

Applying a general function of one function, f(x) to interval, results in this:

f([a, b]) = [ min(f(a), f(b)), max(f(a), f(b)) ]

Applying this to the square root function to intervals:

√([a,b]) = [ min(√a, √b), max(√a, √b) ]

The arithmetic functions for intervals are defined as:

[a, b] + [c, d] = [a + c, b + d]

[a, b] - [c, d] = [a - d, b - c]

[a, b] × [c, d] = [min(a × c, a × d, b × c, b × d), max(a × c, a × d, b × c, b × d)]

[a, b] ÷ [c, d] = [min(a/c, a/d, b/c, b/d), max(a/c, a/d, b/c, b/d)]

Scalar multiplication, given the scalar is positive:

s × [a, b] = [ s × a, s × b ]

HP Prime Interval Programs

INTERVALHELP:  help screen for interval arithmetic functions

INTERVALADD:  adding two intervals

INTERVALSUB:  subtracting two intervals

INTERVALMUL:  multiply two intervals

INTERVALDIV:  divide two intervals

In the Home mode, use the curly brackets (used for lists { }), to type the intervals.  For example, the interval [ 81.8, 82.2 ] would be typed {81.8,82.2}.

EXPORT INTERVALHELP()

BEGIN

// 2022-12-22 EWS

PRINT();

PRINT("Interval Arithmetic");

PRINT("Use curly brackets {}");

PRINT("");

PRINT("INTERVALADD({a,b},{c,d})");

PRINT("add intervals");

PRINT("");

PRINT("INTERVALSUB({a,b},{c,d}");

PRINT("subtract intervals");

PRINT("");

PRINT("INTERVALMUL({a,b},{c,d})");

PRINT("multiply intervals");

PRINT("");

PRINT("INTERVALDIV({a,b},{c,d})");

PRINT("divide intervals");

END;

EXPORT INTERVALADD(v1,v2)

BEGIN

RETURN v1+v2;

END;

EXPORT INTERVALSUB(v1,v2)

BEGIN

RETURN v1-REVERSE(v2);

END;

EXPORT INTERVALMUL(v1,v2)

BEGIN

LOCAL v3;

v3:=CONCAT(v1*v2,v1*REVERSE(v2));

RETURN {MIN(v3),MAX(v3)};

END;

EXPORT INTERVALDIV(v1,v2)

BEGIN

INTERVALMUL(v1,1/REVERSE(v2));

END;

Download a zip file here:  https://drive.google.com/file/d/1L8mmnm3xx58PVjavDROCiZAafBOCsH9o/view?usp=share_link

Example

Interval 1 = {13, 29}

Interval 2 = {10, 14}

Add:  {23, 43}

Subtract:  {-1, 19}

Multiply:  {130, 406}

Divide:  {0.928571428572, 2.9}

Source

Moore, Ramon E.  Mathematical Elements of Scientific Computing  Holt, Rinehart, and Winston, Inc: New York.  1975.  ISBN 0-03-088125-0

Until next time and wish you all well,

Eddie  

All original content copyright, © 2011-2023.  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.