The Museum of
HP Calculators
This program is by Jean-Marc Baillard and is used here by permission.
This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.
Overview
1°) Romberg-Programs
a) Program#1 (
improved )
b) Program#2 (
new )
2°) Arc Length of a Curve
3°) Area of a Surface of Revolution
4°) Area of a 3-dimensional Surface
(
new )
5°) Integrals
a) Program#1
b) Program#2 (
new )
c) Program#3 (
new )
6°) Double Integrals
7°) Triple Integrals
1°) Romberg-Programs
a) Program#1
-The Romberg's method is an "extrapolation to the limit" algorithm which can be applied to many problems:
-Suppose that a sequence In tends to I
as n tends to infinity and that the "errors"
I - In are nearly proportional to 1/nk
If errors were exactly proportional to 1/nk :
I = In + A/nk = I2n + A/(2n)k
, whence I = ( 2k I2n
- In ) / ( 2k - 1 )
When errors are only nearly proportional to 1/nk ,
the errors in this last formula are proportional to 1/nk+2
- "ROM" calculates I1 , I2 , I4
, I8 , ..... ( contimually doubling n )
The formula above is applied to each pair of consecutive In
the results are
J1 , J2 , J4 , ....
and "ROM" applies recursively the previous formula
producing
K1,K2, ....
...., .... until two rounded consecutive approximations
are equal. Thus, accuracy is determined by the display setting.
-This program calls a subroutine which calculates In
for a given n-value. n is stored in register R03 and R03 must be
unchanged by the subroutine.
-The name of this subroutine ( global label, 6 characters
or less ) is to be stored into register R20
-The integer k ( the "order" of the method ) is to be stored
into register R21
Data Registers:
R03: n ( change line 03 and line
10 if you store n in another register )
• R20: the name of the subroutine ( for instance
an integrator ... )
• R21: k = the order of the method used by the subroutine
( k = 2 for the trapezoidal rule,
k = 4 for Simpson's rule,
k = 6 for the 3-point Gauss-Legendre formula ...
>>>> However, even Gaussian integration can behave like a first order
method with singular integrals! )
R22: 2k
R24: counter
R23: 2k , 22k
, ... R25 , R26: .... successive approximations.
Flags: /
Subroutines: It depends on the problem you
have to solved
01 LBL "ROM"
02 1
03 STO 03
04 25
05 STO 24
06 XEQ IND 20
07 STO 25
08 LBL 01
09 2
10 ST* 03
11 RCL 21
12 Y^X
13 STO 22
14 STO 23
15 RCL 24
16 INT
17 E3
18 /
19 25
20 +
21 STO 24
22 XEQ IND 20
23 LBL 02
24 ENTER^
25 ENTER^
26 X<> IND 24
27 -
28 RCL 23
29 RCL 22
30 ST* 23
31 SIGN
32 -
33 /
34 +
35 ISG 24
36 GTO 02
37 STO IND 24
38 VIEW X
39 RND
40 X<>Y
41 RND
42 X#Y?
43 GTO 01
44 RCL IND 24
44 END
( 76 bytes SIZE 026+ ? )
| STACK | INPUTS | OUTPUTS |
| X | / | Limit |
N.B:
-The subroutines called by "ROM" are supposed to keep registers R20
R21 ... unchanged.
-Otherwise, replace for instance these registers by R30 R31 ... in
the "ROM" listing, and replace line 04 and line 19 by 35.
-The Lagrange interpolation formula can also be used to obtain the
same results ( cf the last example at the bottom of this page )
-In fact, "ROM" stops when the rounded approximations obtained with
n = 2 , 4 , 8 , ....... , 2m and with n = 1 , 2 , 4 ,
8 , ....... , 2m are equal.
-This trick generally avoids an unuseful extra-computation with
n = 2m+1
but unfortunately, it may sometimes cause a premature
stop.
-The variant hereafter really compares 2 successive approximations obtained
with n = 1 , 2 , 4 , 8 , ....... , 2m and n = 1
, 2 , 4 , 8 , ....... , 2m+1
b) Program#2
Data Registers:
R03: n ( change line 03 and line
10 if you store n in another register )
• R20: the name of the subroutine ( for instance
an integrator ... )
• R21: k = the order of the method used by the subroutine
R22: 1 , 1/2 , 1/4 , 1/8 , ......
R24 , R25: .... successive approximations.
R23: counter
Flags: /
Subroutines: It depends on the problem you
have to solved
01 LBL "ROM2"
02 1
03 STO 03
04 24
05 STO 23
06 XEQ IND 20
07 STO 24
08 LBL 01
09 2
10 ST* 03
11 SIGN
12 STO 22
13 RCL 23
14 INT
15 E3
16 /
17 24
18 +
19 STO 23
20 XEQ IND 20
21 LBL 02
22 ENTER^
23 X<> IND 23
24 STO Z
25 -
26 2
27 ST/ 22
28 CLX
29 RCL 22
30 RCL 21
31 Y^X
32 1
33 -
34 /
35 -
36 ISG 23
37 GTO 02
38 STO IND 23
39 VIEW X
40 RND
41 X<>Y
42 RND
43 X#Y?
44 GTO 01
45 RCL IND 23
46 END
( 77 bytes / SIZE 025+? )
| STACK | INPUTS | OUTPUTS |
| X | / | Limit |
N.B:
-The termination criterion is of course more rigorous, but the execution
time is longer.
-Moreover, there is now a risk of infinite loop, especially if the
display setting is FIX 9 or SCI 9.
-The following modification is a partial remedy:
Add LASTX 9 / +
after line 39
and ENTER^ after line 37
-Now let's apply the Romberg method to solve a few problems:
( all the following formulas are second-order methods: R21 =
2 )
2°) Arc length of a curve
-The arc length of the curve y = f(x) ( a < x <
b ) is given by L = §ab
( 1 + (y')2 )1/2 dx
"LNG" doesn't use this formula and avoids the calculation of
dy/dx . It simply applies Pythagoras' theorem.
-This is a second order method: k = 2
STO 21
Data Registers:
• R00: f name R04:
L
Registers R00 thru R03 are to be initialized before executing "LNG"
• R01: a
R05: counter ( "ROM" initializes
R03 automatically )
• R02: b
R06: (b-a)/n
• R03: n
R07- R08: temp.
Flags: /
Subroutine: a program which takes
x in X-register and calculates f(x) in X-register.
01 LBL "LNG"
02 RCL 02
03 RCL 01
04 STO 07
05 -
06 RCL 03
07 STO 05
08 /
09 STO 06
10 ST+ 07
11 RCL 01
12 XEQ IND 00
13 STO 08
14 CLX
15 STO 04
16 LBL 01
17 RCL 07
18 XEQ IND 00
19 ENTER^
20 X<> 08
21 -
22 X^2
23 RCL 06
24 ST+ 07
25 X^2
26 +
27 SQRT
28 ST+ 04
29 DSE 05
30 GTO 01
31 RCL 04
32 END
( 48 bytes / SIZE 009 )
| STACK | INPUTS | OUTPUTS |
| X | / | Arc Length |
Example:
-Calculate the arc length of the curve y = ln x 1 < x < 3
LBL "T"
LN
RTN
"LNG" ASTO 20 2 STO
21
"T" ASTO 00
1 STO 01 3 STO 02
FIX 9 XEQ "ROM"
the following approximations are displayed:
2.300459681
2.301931038
2.301986835
2.301987537
2.301987533 ( 76 seconds
) The exact result is 2.301987535
( rounded to 9 decimals )
Notes:
1- If you execute "ROM" in FIX 5 the program gives 2.30199
At this step, if you want a better accuracy, you
don't have to start all over again: simply modify the display format (
for instance FIX 9 ) and press XEQ 01
2-If the equation is given in polar coordinates, replace lines 21-24
by ST+ Z - X^2 X<>Y 2
/ RCL 06 ST+ 07 * and add
ENTER^ after line 19
( this is a second order method too: k = 2
STO 21 )
3-This method can be generalized to parametric equations ...
3°) Area of a surface of revolution
-The rotation of the curve y = f(x) ( a < x <
b ) around x-axis generates a surface of revolution given by
S = 2.pi §ab y ( 1 + (y')2
)1/2 dx
"SRV" doesn't use this formula and avoids the calculation of
dy/dx : the area of a truncated cone is used instead.
-It is also a second order method : k = 2 STO 21
Data Registers:
• R00: f name R04:
S
Registers R00 thru R03 are to be initialized before executing "SRV"
• R01: a
R05: counter ( "ROM"
initializes R03 automatically )
• R02: b
R06: (b-a)/n
• R03: n
R07 R08 temp.
Flags: /
Subroutine: a program which takes
x in X-register and calculates f(x) in X-register.
01 LBL "SRV"
02 RCL 02
03 RCL 01
04 STO 07
05 -
06 RCL 03
07 STO 05
08 /
09 STO 06
10 ST+ 07
11 RCL 01
12 XEQ IND 00
13 STO 08
14 CLX
15 STO 04
16 LBL 01
17 RCL 07
18 XEQ IND 00
19 ENTER^
20 ENTER^
21 X<> 08
22 ST+ Z
23 -
24 X^2
25 RCL 06
26 ST+ 07
27 X^2
28 +
29 SQRT
30 *
31 ST+ 04
32 DSE 05
33 GTO 01
34 PI
35 ST* 04
36 RCL 04
37 END
( 55 bytes / SIZE 009 )
| STACK | INPUTS | OUTPUTS |
| X | / | Area |
Example:
-Evaluate the area of the surface of revolution generated by the rotation of the curve y = sin x ( 0 < x < pi ) around the x-axis.
LBL "T"
SIN
RTN
"SRV" ASTO 20 2 STO
21
"T" ASTO 00
0 STO 01 PI STO 02
set the HP-41 in RAD mode ( XEQ RAD )
FIX 9 XEQ "ROM" the
following approximations are displayed:
15.59985804
14.26493243
14.42990383
14.42356623
14.42359884
14.42359950
( 177 seconds )
4°) Area of a 3-Dimensional Surface
-This routine computes the area of a 3D-surface defined by: z = f(x,y) a < x < b , c < y < d
-The result could be obtained by the double integral §ab §cd ( 1 + fx2 + fy2 )1/2 dx dy
where fx = ¶f/¶x and fy = ¶f/¶y are the partial derivatives with respect to x and y respectively.
-But "SKS" avoids the calculus of the partial derivatives:
-The intervals [a,b] and [c,d] are divided into n parts, and the approximate
area is the sum of the areas of triangles obtained by Heron's formula.
-This is a second order method.
Data Registers:
• R00: f name •
R04: c
Registers R00 thru R05 are to be initialized before executing "SKS"
• R01: a
• R05: d
( "ROM" initializes R03 automatically )
• R02: b
R06: Area
• R03: n
R07 thru R18 temp.
Flags: /
Subroutine: A program
which takes x in X-register and y in Y-register and calculates f(x,y) in
X-register.
01 LBL "SKS"
02 RCL 02
03 RCL 01
04 STO 09
05 -
06 STO 07
07 RCL 05
08 RCL 04
09 -
10 RCL 03
11 ST/ 07
12 STO 17
13 /
14 STO 08
15 CLX
16 STO 06
17 LBL 01
18 RCL 03
19 STO 18
20 RCL 04
21 STO 10
22 RCL 09
23 XEQ IND 00
24 STO 11
25 RCL 10
26 RCL 09
27 RCL 07
28 +
29 XEQ IND 00
30 STO 12
31 LBL 02
32 RCL 08
33 ST+ 10
34 RCL 10
35 RCL 09
36 RCL 07
37 +
38 XEQ IND 00
39 STO 14
40 RCL 10
41 RCL 09
42 XEQ IND 00
43 STO 13
44 RCL 12
45 -
46 X^2
47 RCL 07
48 X^2
49 +
50 RCL 08
51 X^2
52 +
53 SQRT
54 STO 15
55 RCL 12
56 RCL 11
57 -
58 X^2
59 RCL 07
60 X^2
61 +
62 SQRT
63 STO 16
64 +
65 RCL 13
66 RCL 11
67 -
68 X^2
69 RCL 08
70 X^2
71 +
72 SQRT
73 STO Z
74 +
75 2
76 /
77 ENTER^
78 ENTER^
79 R^
80 -
81 *
82 X<>Y
83 RCL 15
84 -
85 *
86 X<>Y
87 RCL 16
88 -
89 *
90 SQRT
91 ST+ 06
92 RCL 14
93 RCL 12
94 -
95 X^2
96 RCL 08
97 X^2
98 +
99 SQRT
100 STO 16
101 RCL 14
102 STO 12
103 RCL 13
104 STO 11
105 -
106 X^2
107 RCL 07
108 X^2
109 +
110 SQRT
111 STO Z
112 +
113 RCL 15
114 +
115 2
116 /
117 ENTER^
118 ENTER^
119 ST- 15
120 R^
121 -
122 *
123 RCL 15
124 *
125 RCL 16
126 R^
127 -
128 *
129 SQRT
130 ST+ 06
131 DSE 18
132 GTO 02
a three-byte GTO
133 RCL 07
134 ST+ 09
135 DSE 17
136 GTO 01
a three-byte GTO
137 RCL 06
138 END
( 170 bytes / SIZE 019 )
| STACK | INPUTS | OUTPUTS |
| X | / | Area |
Example: Calculate the area of the surface
defined by z = exp ( -x2-y )
0 < x < 2 , 0 < y < 3
LBL "Z"
X^2
+
CHS
E^X
RTN
"Z" ASTO 00
0 STO 01 STO 04
2 STO 02
3 STO 05
"SKS" ASTO 20
-This is a second order method: 2 STO
21
-For instance FIX 3 XEQ "ROM" the following approximations are displayed:
6.262
6.284
6.284 ( in
6mn46s )
-The last approximation is actually 6.283584
in X-register
whereas the exact result is 6.2833787 rounded to
7 decimals.
5°) Integrals §ab
f(x) dx
a) Program#1
-The following program uses the very simple "midpoint" formula: §ab f(x) dx = (b-a).f((a+b)/2) It is a second order method ( k = 2 STO 21 )
-It does not have the refinements of the PPC ROM program "IG" ( cf PPC-ROM
user's manual page 220-225 )
therefore, the convergence may be disappointing for singular
integrals,
but it's a good example of how the Romberg's method can produce
superb acceleration.
-Furthermore, it's probably one of the shortest integration program
for the HP-41.
Data Registers:
• R00: f name R04:
§ab f(x) dx
Registers R00 thru R03 are to be initialized before executing "IG"
• R01: a
R05: counter
( "ROM" initializes R03 automatically )
• R02: b
R06: (b-a)/n
• R03: n
R07 temp.
Flags: /
Subroutine: a program which takes
x in X-register and calculates f(x) in X-register.
01 LBL "IG"
02 RCL 02
03 RCL 01
04 STO 07
05 -
06 RCL 03
07 STO 05
08 /
09 STO 06
10 2
11 /
12 ST+ 07
13 CLX
14 STO 04
15 LBL 01
16 RCL 07
17 XEQ IND 00
18 ST+ 04
19 RCL 06
20 ST+ 07
21 DSE 05
22 GTO 01
23 ST* 04
24 RCL 04
25 END
( 40 bytes / SIZE 008 )
| STACK | INPUTS | OUTPUTS |
| X | / | Integral |
Example:
-Evaluate I = §12 xx dx to 7 decimals
LBL "T"
ENTER^
Y^X
RTN
"IG" ASTO 20 2 STO 21
"T" ASTO 00 1
STO 01 2 STO 02
FIX 7 XEQ "ROM" the following values
are displayed:
2.0438805
2.0503885
2.0504461
2.0504462 ( in 40 seconds
)
-In fact, the result is correct to almost 9 decimals: I = 2.050446234
( the last digit should be a 5 )
b) Program#2
-This routine uses the Romberg method with the trapezoidal rule.
-So, f(a) and f(b) must be finite.
-The previous sums are used the next iterations.
-The Romberg extrapolation is performed inside the program itself:
it doesn't use "ROM"
Data Registers: R00 thru R09 are unused ( Register R10 is to be initialized before executing "IG2" )
• R10: f name R13-R14:
temp
R11: a
R15: counter
R17: 22 , 42 , 82 , ....
R12: (b-a)/2n
R16: 2n
R18: counter
R19 , R20 , .... = successive approximations
Flags: /
Subroutine: a program which calculates
f(x) assuming x is in X-register upon entry.
01 LBL "IG2"
02 STO 12
03 X<>Y
04 STO 11
05 X<>Y
06 XEQ IND 10
07 STO 14
08 RCL 11
09 ST- 12
10 XEQ IND 10
11 RCL 14
12 +
13 2
14 /
15 STO 14
16 RCL 12
17 *
18 STO 19
19 19
20 STO 18
21 SIGN
22 STO 16
23 LBL 01
24 RCL 16
25 STO 15
26 RCL 11
27 RCL 12
28 2
29 /
30 +
31 LBL 02
32 STO 13
33 XEQ IND 10
34 ST+ 14
35 RCL 12
36 RCL 13
37 +
38 DSE 15
39 GTO 02
40 2
41 ST/ 12
42 ST* 16
43 X^2
44 STO 17
45 RCL 18
46 INT
47 E3
48 /
49 19
50 +
51 STO 18
52 RCL 12
53 RCL 14
54 *
55 LBL 03
56 ENTER^
57 ENTER^
58 X<> IND 18
59 STO 13
60 -
61 RCL 17
62 4
63 ST* 17
64 SIGN
65 -
66 /
67 +
68 ISG 18
69 GTO 03
70 STO IND 18
71 VIEW X
72 RND
73 RCL 13
74 RND
the accuracy is controlled by the display setting
75 X#Y?
the iterations stop when 2 rounded successive approximations are equal.
76 GTO 01
77 RCL IND 18
78 END
( 114 bytes / SIZE 021+? )
| STACK | INPUTS | OUTPUTS |
| Y | a | / |
| X | b | §ab f(x) dx |
Example: Evaluate §13 x ( 1+x3 )1/2 dx
LBL "T" X^2
SIGN *
"T" ASTO 10
ENTER^ *
+ RTN
ENTER^ ENTER^ SQRT
FIX 7 1 ENTER^
3 XEQ "IG2" >>>> 13.7629072
13.7691196
13.7693295
13.7693320
13.7693320 ( in 36 seconds )
-The last value is in fact 13.76933202 all
the decimals are correct.
c) Program#3
-Like "IG" , "IG3" uses the midpoint formula,
but in order to use the previous evaluations of the function,
the number of steps is trippled with each iteration.
Data Registers: R00 thru R09 are unused ( Register R10 is to be initialized before executing "IG3" )
• R10: f name
R13: 3n
R16: counter and 32 , 92 , 272 ,
....
R11: a
R14: sums
R17: counter
R12: (b-a)/3n
R15: x
R18, R19 , .... = successive approximations
Flags: /
Subroutine: a program which calculates
f(x) assuming x is in X-register upon entry.
01 LBL "IG3"
02 STO 12
03 X<>Y
04 STO 11
05 ST- 12
06 +
07 2
08 /
09 XEQ IND 10
10 STO 14
11 RCL 12
12 *
13 STO 18
14 18
15 STO 17
16 SIGN
17 STO 13
18 LBL 01
19 RCL 13
20 STO 16
21 3
22 ST/ 12
23 ST* 13
24 RCL 11
25 RCL 12
26 2
27 /
28 +
29 LBL 02
30 STO 15
31 XEQ IND 10
32 ST+ 14
33 RCL 12
34 ST+ X
35 ST+ 15
36 RCL 15
37 XEQ IND 10
38 ST+ 14
39 RCL 12
40 RCL 15
41 +
42 DSE 16
43 GTO 02
44 9
45 STO 16
46 RCL 17
47 INT
48 E3
49 /
50 18
51 +
52 STO 17
53 RCL 12
54 RCL 14
55 *
56 LBL 03
57 ENTER^
58 ENTER^
59 X<> IND 17
60 STO 15
61 -
62 RCL 16
63 9
64 ST* 16
65 SIGN
66 -
67 /
68 +
69 ISG 17
70 GTO 03
71 STO IND 17
72 VIEW X
73 RND
74 RCL 15
75 RND
the accuracy is controlled by the display setting
76 X#Y?
the iterations stop when 2 rounded successive approximations are equal.
77 GTO 01
78 RCL IND 17
79 END
( 117 bytes / SIZE 020+? )
| STACK | INPUTS | OUTPUTS |
| Y | a | / |
| X | b | §ab f(x) dx |
Example: Evaluate §13 x ( 1+x3 )1/2 dx
LBL "T" X^2
SIGN *
"T" ASTO 10
ENTER^ *
+ RTN
ENTER^ ENTER^ SQRT
FIX 7 1 ENTER^
3 XEQ "IG3" >>>> 13.7718432
13.7693525
13.7693320
13.7693320 ( in 64 seconds )
-The last value is 13.76933201 ( the last decimal
should be a "2" )
-Actually, the next to last result was exact: 13.76933202
( in R15 and in L-register )
-Unlike "IG" , "IG3" doesn't loose the previous function evaluations,
but step trippling often leads to long execution times.
6°) Double Integrals §ab
§u(x)v(x)
f(x,y) dx dy
-The mid-point formula is used in both x and y-axis.
( a second order method : k = 2 STO 21 )
Data Registers:
• R00: f name •
R04: u name
Registers R00 thru R05 are to be initialized before executing "DIN"
• R01: a
• R05: v name
( "ROM" initializes R03 automatically )
• R02: b
R06 = §§
• R03: n
R07 thru R13: temp.
Flags: /
Subroutines: A program
which takes x in X-register and y in Y-register and calculates f(x,y) in
X-register.
A program which takes x in X-register and calculates u(x) in X-register.
A program which takes x in X-register and calculates v(x) in X-register.
01 LBL "DIN"
02 RCL 02
03 RCL 01
04 STO 09
05 -
06 RCL 03
07 STO 07
08 /
09 STO 08
10 2
11 /
12 ST+ 09
13 CLX
14 STO 06
15 LBL 01
16 RCL 09
17 XEQ IND 04
18 STO 10
19 RCL 09
20 XEQ IND 05
21 RCL 10
22 -
23 RCL 03
24 STO 11
25 /
26 STO 12
27 2
28 /
29 ST+ 10
30 CLX
31 STO 13
32 LBL 02
33 RCL 10
34 RCL 09
35 XEQ IND 00
36 ST+ 13
37 RCL 12
38 ST+ 10
39 DSE 11
40 GTO 02
41 RCL 13
42 *
43 ST+ 06
44 RCL 08
45 ST+ 09
46 DSE 07
47 GTO 01
48 ST* 06
49 RCL 06
50 END
( 72 bytes / SIZE 014 )
| STACK | INPUTS | OUTPUTS |
| X | / | Integral |
Example:
-Calculate §23 §xx^2 ( 1 + xy )1/2 dx dy to 7 decimals
LBL "T"
*
1
+
SQRT
RTN
LBL "U"
RTN
LBL "V"
X^2
RTN
"DIN" ASTO 20 2 STO
21
"T" ASTO 00
2 STO 01 3 STO 02
"U" ASTO 04 "V" ASTO
05
FIX 7 XEQ "ROM" the HP-41 displays:
13.7800635
13.7747491
13.7746576
13.7746565 ( 263 seconds
)
7°) Triple Integrals §ab
§u(x)v(x)
§w(x,y)t(x,y)
f(x,y,z) dx dy dz
-The mid-point formula is again used ( in x, y and z-axis ).
( a second order method : k = 2 STO 21 )
Data Registers:
• R00: f name •
R04: u name
Registers R00 thru R05 are to be initialized before executing "TIN"
• R01: a
• R05: v name
( "ROM" initializes R03 automatically )
• R02: b
• R06: w name
• R03: n
• R07 t name
R08 = §§§ , R09 thru R19 temp.
Flags: /
Subroutines: A program which
takes x in X-register, y in Y-register and z in Z-register and calculates
f(x,y,z) in X-register.
A program which takes x in X-register and calculates u(x) in X-register.
A program which takes x in X-register and calculates v(x) in X-register.
A program which takes x in X-register and y in Y-register and calculates
w(x,y) in X-register.
A program which takes x in X-register and y in Y-register and calculates
t(x,y) in X-register.
01 LBL "TIN"
02 RCL 02
03 RCL 01
04 STO 11
05 -
06 RCL 03
07 STO 09
08 /
09 STO 10
10 2
11 /
12 ST+ 11
13 CLX
14 STO 08
15 LBL 01
16 RCL 11
17 XEQ IND 04
18 STO 12
19 RCL 11
20 XEQ IND 05
21 RCL 12
22 -
23 RCL 03
24 STO 13
25 /
26 STO 14
27 2
28 /
29 ST+ 12
30 CLX
31 STO 15
32 LBL 02
33 RCL 12
34 RCL 11
35 XEQ IND 06
36 STO 16
37 RCL 12
38 RCL 11
39 XEQ IND 07
40 RCL 16
41 -
42 RCL 03
43 STO 17
44 /
45 STO 18
46 2
47 /
48 ST+ 16
49 CLX
50 STO 19
51 LBL 03
52 RCL 16
53 RCL 12
54 RCL 11
55 XEQ IND 00
56 ST+ 19
57 RCL 18
58 ST+ 16
59 DSE 17
60 GTO 03
61 RCL 19
62 *
63 ST+ 15
64 RCL 14
65 ST+ 12
66 DSE 13
67 GTO 02
68 RCL 15
69 *
70 ST+ 08
71 RCL 10
72 ST+ 11
73 DSE 09
74 GTO 01
75 ST* 08
76 RCL 08
77 END
( 114 bytes / SIZE 020 )
| STACK | INPUTS | OUTPUTS |
| X | / | Integral |
Example: Evaluate §01 §0x §-x-y-y 1 / ( 1 + x + y + z ) dx dy dz to 5 places
LBL "P"
+
+
1
+
1/X
RTN
LBL "U"
CLX
RTN
LBL "V"
RTN
LBL "W"
+
CHS
RTN
LBL "T"
X<>Y
CHS
RTN
"TIN" ASTO 20 2
STO 21
"P" ASTO 00
0 STO 01 1 STO 02
"U" ASTO 04
V ASTO 05
"W" ASTO 06
T ASTO 07
FIX 5 XEQ "ROM" the calculator
displays:
0.24838
0.24998
0.25000 ( 9 minutes )
( exact value = 1/4 )
The last approximation is actually 0.249999724 and was obtained from the following results:
with n = 1 "TIN" yields 0.200000000
---- n = 2 ------------ 0.236284830
---- n = 4 ------------ 0.246481100
---- n = 8 ------------ 0.249114231
-If you use these 4 numbers as explained in "Lagrange Interpolation
Formula for the HP-41" (example2) you'll get 0.249999724 too
-With n = 16 , the execution time increases very much , but you could
perhaps XEQ "TIN" with n = 10 ( STO 03 ):
the result is 0.249432635
and then, the Lagrange polynomial gives a slightly better accuracy:
with these 5 numbers, you'll find: L(0) = 0.249999997
-"ROM" has the advantage to be an automatic process, but execution time
may become a problem because n is continually doubled.
- In the Lagrange approach, you have the opportunity to choose the
n-values by yourself in a possibly more economical way.
References:
[1] "Numerical Analysis" - F. Scheid - Mc Graw Hill
ISBN 07-055197-9
[2] "PPC-ROM user's manual" ( page 220 ..... )
[3] "Extended Functions made easy" - Keith Jarett ( synthetix
)
[4] "Numerical Recipes in C++" - Press , Vetterling , Teukolsky
, Flannery - Cambridge University Press - ISBN 0-521-75033-4
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall