16 Point Gaussian Quadrature for the HP-41C/CV/CX (updated)

[Namir]

16 Point Gaussian Quadrature for the HP-41C/CV/CX (updated)

This program is a slightly enhanced version of the one posted by Tony Udell. The program includes LBL 00 that initializes the values for the quadrature nodes and weights. I have also remapped the memory registers into a contiguous sequence from 00 to 19.

You can download the .raw program file by clicking here.

To use the program:

1) Press [A] to perform a regular integration between A and B. The program prompts you for the values of A and B. Enter these values and press [R/S] to calculate the integral.

2) Press [B] to perform a special integration between A and infinity. The program prompts you for the value of A. Enter the value and press [R/S] to calculate the integral.

Notes:
1) Label "FX" (also label E) is where you enter the statements to evaluate the integrated function.
2) The program automatically calls subroutine at LBL 00 to initialize the nodes and weights for the quadrature.


Memory Map
==========

Code:

R00 = integral
R00 = integral
R01 = a
R02 = b
R03 to R18 x(1),w(1),x(2),w(2),...x(8),w(8)
R19 = index used to access x(i) and w(i)

Flags used 00 and 01

Program Listing
============

Code:

01 LBL "A-B"
02 LBL A
03 XEQ 00
04 FIX 05
05 "a?"
06 PROMPT
07 STO 01
08 "b?"
09 PROMPT
10 STO 02
11 0
12 STO 00
13 CF 00
14 LBL 01
15 18
16 STO 19
17 XEQ 15
18 XEQ 15
19 XEQ 15
20 XEQ 15
21 XEQ 15
22 XEQ 15
23 XEQ 15
24 XEQ 15
25 FS? 00
26 GTO 02
27 SF 00
28 GTO 01
29 LBL 02
30 RCL 02
31 RCL 01
32 -
33 2
34 ÷
35 RCL 00
36 ×
37 RTN
38 LBL 15
39 RCL IND 19
40 FS? 00
41 +/-
42 RCL 02
43 RCL 01
44 -
45 ×
46 RCL 02
47 +
48 RCL 01
49 +
50 2
51 ÷
52 XEQ "FX"
53 DSE 19
54 RCL IND 19
55 DSE 19
56 ×
57 STO+ 00
58 RTN
59 LBL "A-"
60 LBL B
61 XEQ 00
62 FIX 05
63 "a?"
64 PROMPT
65 STO 01
66 0
67 STO 00
68 CF 00
69 LBL 03
70 18
71 STO 19
72 XEQ 16
73 XEQ 16
74 XEQ 16
75 XEQ 16
76 XEQ 16
77 XEQ 16
78 XEQ 16
79 XEQ 16
80 FS? 00
81 GTO 04
82 SF 00
83 GTO 03
84 LBL 04
85 RCL 00
86 2
87 ×
88 RTN
89 LBL 16
90 RCL IND 19
91 FS? 00
92 CHS
93 1
94 +
95 2
96 X<>Y
97 ÷
98 RCL 01
99 +
100 1
101 -
102 XEQ "FX"
103 RCL IND 19
104 FS? 00
105 CHS
106 1
107 +
108 X^2
109 DSE 19
110 RCL IND 19
111 DSE 19
112 X<>Y
113 ÷
114 ×
115 STO+ 00
116 RTN
117 LBL 00
118 FS? 01
119 RTN
120 2.715245941E-2
121 STO 03
122 9.89400935E-1
123 STO 04
124 6.225352394E-2
125 STO 05
126 9.445750231E-1
127 STO 06
128 9.515851168E-2
129 STO 07
130 8.656312024E-1
131 STO 08
132 1.246289713E-1
133 STO 09
134 7.554044084E-1
135 STO 10
136 1.495959888E-1
137 STO 11
138 6.178762444E-1
139 STO 12
140 1.691565194E-1
141 STO 13
142 4.580167777E-1
143 STO 14
144 1.82603415E-1
145 STO 15
146 2.816035508E-1
147 STO 16
148 1.894506105E-1
149 STO 17
150 9.501250984E-2
151 STO 18
152 SF 01
153 RTN
154 LBL "FX"
155 LBL E
156 1/X
157 RTN

Example 1
=========

To calculate the integral of 1/x from 1 to 2 (which is equal to ln(2)):

1. In program mode insert the command 1/x after LBL E and make sure it is followed by RTN.
2. In run mode (user mode on) clear flag 1 using [f][CF][0][1]. Skip this step if you are not running the program the first time.
3. Press the [A] key.
4. The program prompts you to enter the value of "a". Enter 1 and press the [R/S] key.
5. The program prompts you to enter the value of "b". Enter 2 and press the [R/S] key.
6. The program displays 0.69315 (FIX 9 displays 0.693147181) as the value of the integral.

Example 2
========

To calculate the integral of exp(-x)*x^0.8 from 0 to inifinity (which is equal to gamma(1.8)):

1. In program mode insert the commands CHS, EXP, LASTX, CHS, 0.8, Y^X and, * after LBL E and make sure it is followed by RTN.
2. In run mode (user mode on) clear flag 1 using [f][CF][0][1]. Skip this step if you are not running the program the first time.
3. Press the [A] key.
4. The program prompts you to enter the value of "a". Enter 0 and press the [R/S] key.
5. The program displays 0.93138 (FIX 9 displays 0.931378447) as the value of the integral.

[Albert Chan]

(11-09-2023 09:23 PM)Namir Wrote:  2) Press [B] to perform a special integration between A and infinity.

What transform is done to the integrand, to map limits (A .. ∞) to (-1 .. 1) ?

[Valentin Albillo]

.
Hi, Namir,

There's a number of errors in your program's listing, namely wrong numeration of the lines marked in red below:

Namir Wrote:
105 CHS
106 1
101 +
102 X^2
109 DSE 19
110 RCL IND 19
[...]
123 STO 04
124 6.219352394E-2
119 STO 05
126 9.445750231E-1
127 STO 06

Besides, several of your 10-digit numeric coefficients below marked in red are wrong or inaccurate ...

Namir Wrote:120 2.715245941E-2
121 STO 03
122 9.89400935E-1
123 STO 04
124 6.219352394E-2
119 STO 05
126 9.445750231E-1
127 STO 06
128 9.515851168E-2
129 STO 07
130 8.656312024E-1
131 STO 08
132 1.246289713E-1
133 STO 09
134 7.554044024E-1
135 STO 10
136 1.495959888E-1
137 STO 11
138 6.178762444E-1
139 STO 12
140 1.691565194E-1
141 STO 13
142 4.580167777E-1
143 STO 14
144 1.82603415E-1
145 STO 15
146 2.816035502E-1
147 STO 16
148 1.894506105E-1
149 STO 17
150 9.501190984E-2
151 STO 18

... so you'd do well getting the correct/accurate values from a reputable source and/or typing them correctly.

Finally, if you'd take a friendly advice, when posting a program which performs such a computation, it would be proper (nay, mandatory) to include some numeric results generated by the program, if nothing else just that people trying it out can be reassured that it's correctly loaded and properly working by simply running your program and actually getting exactly your results.

Regards.
V.

[Namir]

Thanks Valentin. The errors were generated when I was renaming registers. This process afected some line numbers and some of the numerical constants. Fortunately, I wrote the program first, saved it as a .raw file, read it with a free42 emulator, print it, and then edited it. So the .raw file in the link should be correct.

Will be adding examples for using the program.

Again, Thanks!

Namir

[Namir]

(11-09-2023 10:39 PM)Albert Chan Wrote:  

(11-09-2023 09:23 PM)Namir Wrote:  2) Press [B] to perform a special integration between A and infinity.

What transform is done to the integrand, to map limits (A .. ∞) to (-1 .. 1) ?

Th approximation for the integral uses:

integral = 2 * sum of (w(i)/(1+x(i))^2 * fx(2/(1+x(i)) + A - 1)) for i=1,..., 6

where x(i) are the node values and w(i) are the corresponging weights. The limit 6 can be different, adding more points to get more accurate approximations for the integral.

[Namir]

I added two examples.

[rawi]

Hi Namir,

thank you for sharing. Example 1 works fine but I cannot get to a result of example 2.

When using example 2 with SST at line 101 I get -0.89474 in the x-register. This goes to the function FX resulting in (-0.89474)^0.8 terminating the program with "DATA ERROR". Can you help?

Thank you very much.

[Namir]

(11-10-2023 05:28 PM)rawi Wrote:  Hi Namir,

thank you for sharing. Example 1 works fine but I cannot get to a result of example 2.

When using example 2 with SST at line 101 I get -0.89474 in the x-register. This goes to the function FX resulting in (-0.89474)^0.8 terminating the program with "DATA ERROR". Can you help?

Thank you very much.

Did you deete the command 1/X in the original listing, before you enter the commands for Example 2?

You might also want to try downloading the .raw file and importing it into the HP-41CX emulator from TOS. You can also try to import the code with free42. In the latter case you may want to edit labels A and B to add another letter/digit like AA and BB pr A1 and B1.

Namir

[Valentin Albillo]

(11-10-2023 03:18 AM)Namir Wrote:  Will be adding examples for using the program.

Again, Thanks!

You're welcome. Also, you posted these two results:

Quote:Example 1
=========

To calculate the integral of 1/x from 1 to 2 (which is equal to ln(2)):
[...]
The program displays 0.69315 as the value of the integral,

Example 2
========

To calculate the integral of exp(-x)*x^0.8 from 0 to inifinity (which is equal to gamma(1.8)):
[...]
The program displays 0.93138 as the value of the integral.

Both were posted as 5-digit values. Could you please post these two results to full 10-digit accuracy ? It would be necessary to gauge the correctness of your program/data in order to detect any further errors and, again, for people using your program to check that they do exactly obtain your results to full 10-digit precision.

Regards.
V.

[rawi]

Thanks for your fast reply!

Namir wrote:

Quote:Did you deete the command 1/X in the original listing, before you enter the commands for Example 2?

Yes, I did.

Quote:You might also want to try downloading the .raw file and importing it into the HP-41CX emulator from TOS. You can also try to import the code with free42. In the latter case you may want to edit labels A and B to add another letter/digit like AA and BB pr A1 and B1.

I tried to download but got a message, that it cannot be downloaded safely.

[Namir]

(11-11-2023 08:38 AM)rawi Wrote:  Thanks for your fast reply!

Namir wrote:

Quote:Did you deete the command 1/X in the original listing, before you enter the commands for Example 2?

Yes, I did.

Quote:You might also want to try downloading the .raw file and importing it into the HP-41CX emulator from TOS. You can also try to import the code with free42. In the latter case you may want to edit labels A and B to add another letter/digit like AA and BB pr A1 and B1.

I tried to download but got a message, that it cannot be downloaded safely.

I checked the listed code (for the integration between A and infinity) with that in my 41CX emulator and found no errors or differences.

Namir

[John Keith]

(11-11-2023 08:38 AM)rawi Wrote:  I tried to download but got a message, that it cannot be downloaded safely.

Some browsers will give that message for almost any site that is not a large commercial entity. If you trust the site, you can usually override the message. If not, try another browser.

[rawi]

Thank you Namir and John,

now the download of the .raw file was successful and the program works fine.
Especially the integration to infinity is somethong that is very useful since it is not included in the integration programs used in HP 15C, HP 42S or the Advantage module. And it is very exact: The upper tail of the normal distribution from 1.96 to infinity is computed as 0.02499792 whereas the exact value is 0.02499790.

[Namir]

(11-11-2023 03:41 PM)rawi Wrote:  Thank you Namir and John,

now the download of the .raw file was successful and the program works fine.
Especially the integration to infinity is somethong that is very useful since it is not included in the integration programs used in HP 15C, HP 42S or the Advantage module. And it is very exact: The upper tail of the normal distribution from 1.96 to infinity is computed as 0.02499792 whereas the exact value is 0.02499790.

Good to hear that! Thanks!

Namir

[Valentin Albillo]

.
Hi, Namir,

Would you please post the two values as I asked in Message #9 above ?

Thanks in advance.
V.

[Namir]

(11-12-2023 01:49 AM)Valentin Albillo Wrote:  .
Hi, Namir,

Would you please post the two values as I asked in Message #9 above ?

Thanks in advance.
V.

I updated the original post to indicate what FIX 9 would show for both examples.

[Valentin Albillo]

(11-12-2023 04:22 AM)Namir Wrote:  

(11-12-2023 01:49 AM)Valentin Albillo Wrote:  Would you please post the two values as I asked in Message #9 above ?

I updated the original post to indicate what FIX 9 would show for both examples.

Yes, I noticed, but I need to know all 10-digit mantissas, not just 9-digit. A simple way to obtain them in FIX 9 is once your program outputs rhe result, multiply it by 10 (i.e. 10, *). This will show all 10 digits present. Please do and post a new message here with both 10-digit results, don't update the original post.

Again, thanks in advance.
V.

[Namir]

(11-12-2023 05:32 AM)Valentin Albillo Wrote:  

(11-12-2023 04:22 AM)Namir Wrote:  I updated the original post to indicate what FIX 9 would show for both examples.

Yes, I noticed, but I need to know all 10-digit mantissas, not just 9-digit. A simple way to obtain them in FIX 9 is once your program outputs rhe result, multiply it by 10 (i.e. 10, *). This will show all 10 digits present. Please do and post a new message here with both 10-digit results, don't update the original post.

Again, thanks in advance.
V.

Example 1 0.6931471810
Example 2 0.9313784470

[Valentin Albillo]

(11-12-2023 08:38 AM)Namir Wrote:  

(11-12-2023 05:32 AM)Valentin Albillo Wrote:  Yes, I noticed, but I need to know all 10-digit mantissas, not just 9-digit. A simple way to obtain them in FIX 9 is once your program outputs rhe result, multiply it by 10 (i.e. 10, *). This will show all 10 digits present. Please do and post a new message here with both 10-digit results, don't update the original post.

Again, thanks in advance.
V.

Example 1 0.6931471810
Example 2 0.9313784470

Both values end in 0 ? Really ? It's a 1-in-100 chance. Did you multiply them by 10 as I suggested ? If yes, then they'd be like 6.93... and 9.31... respectively, not the 0.693... and 0.931... values you posted.
Could you please check them out one final time ?

Thanks and regards.
V.

[Namir]

(11-12-2023 03:24 PM)Valentin Albillo Wrote:  

(11-12-2023 08:38 AM)Namir Wrote:  Example 1 0.6931471810
Example 2 0.9313784470

Both values end in 0 ? Really ? It's a 1-in-100 chance. Did you multiply them by 10 as I suggested ? If yes, then they'd be like 6.93... and 9.31... respectively, not the 0.693... and 0.931... values you posted.
Could you please check them out one final time ?

Thanks and regards.
V.

Yes I multiplied both results by 10 ans then moved the decimal to the right, so they correspond to the original answers (now with 10 digits)..

Using MATLAB, I calculated log(2) (the integral of 1/x from 1 to 2) and got 0.693147180559945. Comparing this result with what I got on the HP-41CX, looks like the 41CX rounded the part of 8055 into 810. I assume the calculator might have done similar rounding to the result of Example 2.