The Museum of HP Calculators
The Museum of HP Calculators
This program is Copyright © 2006 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.
1°) One-dimensional Problem
2°) Two-dimensional Problem
3°) N-dimensional Problem
-The 3 following programs find a minimum of a given function
f.
-Search for a minimum of (-f) if you want to know a maximum of the
function.
1°) One-dimensional Problem
Data Registers:
• R00 = function name ( this register is to be initialized
before executing "EX" )
R01 = x ; R02 = f(x) ; R03 = h ; R04 = f(x+h)
; R05 = f(x-h)
Flags: /
Subroutine: A program which computes
f(x) assuming x is in X-register upon entry.
01 LBL "EX"
02 X<>Y
03 STO 03
04 X<>Y
05 STO 01
06 XEQ IND 00
07 STO 02
08 LBL 01
09 RCL 01
10 RCL 03
11 +
12 XEQ IND 00
13 STO 04
14 LBL 02
15 RCL 01
16 RCL 03
17 -
18 XEQ IND 00
19 STO 05
20 LBL 03
21 VIEW 01
22 RCL 04
23 X>Y?
24 X<>Y
25 RCL 02
the program compares f(x) , f(x+h) , f(x-h)
26 X>Y?
if f(x+h) or f(x-h) is inferior to f(x),
27 X<>Y
x+h or x-h becomes the "new" x
28 RCL 02
29 X=Y?
30 GTO 05
31 CLX
32 RCL 04
33 X=Y?
34 GTO 04
35 RCL 05
36 X<> 02
37 STO 04
38 RCL 03
39 ST- 01
40 GTO 02
41 LBL 04
42 X<> 02
43 STO 05
44 RCL 03
45 ST+ 01
46 RCL 01
47 +
48 XEQ IND 00
49 STO 04
50 RCL 05
51 GTO 03
52 LBL 05
if f(x) is inferior to f(x+h) and f(x-h)
53 PI
then h is divided by PI
54 ST/ 03
( I don't know which value is the best statistically )
55 RCL 03
56 RND
the iteration stops when the rounded h-value equals zero
57 X#0?
so, the accuracy is controlled by the display format.
58 GTO 01
59 RCL 02
60 RCL 01
61 END
( 83 bytes / SIZE 006 )
| STACK | INPUTS | OUTPUTS |
| Y | h | min(f) = f(x) |
| X | x0 | x |
Example: Find the minimum of the function f(x) = | x2 - 2 |
-We load the short routine
LBL "T" ( any global label, maximum 6 characters
)
X^2
2
-
ABS
RTN in program memory, "T" ASTO
00 and if we choose x0 =
0 as initial guess and h = 1
FIX 9
1 ENTER^
0 XEQ "EX" >>>> the successive x-values are
displayed and eventually,
we have: x = 1.141213652 = R01
( execution time = 46 seconds )
RDN fmin
= 10 -9
= R02
-Here, the solution was trivial: x = SQRT(2) , fmin = 0
y
| *
| *
*
*
|
* *
* *
|
*
* *
|
* *
|
*
--|--------x1-------------x2----------------------------
x
-If the function has the graph above and if you choose x = 0
as initial guess,
the program will of course give the solution x = x1
-Choose another guess to find x2
2°) Two-dimensional Problem
Data Registers:
• R00 = function name ( this register is to be initialized
before executing "EXY" )
R01 = x ; R02 = y ; R03 = f(x) ; R04 = h
; R05: temp
Flags: /
Subroutine: A program which computes
f(x,y) assuming x is in X-register and y is in Y-register upon entry.
01 LBL "EXY"
02 X<> Z
03 STO 04
04 RDN
05 STO 02
06 X<>Y
07 STO 01
08 XEQ IND 00
09 STO 03
10 LBL 01
11 VIEW 01
12 3
13 STO 05
14 RCL 02
15 RCL 01
16 RCL 04
17 +
18 XEQ IND 00
19 RCL 03
20 X>Y?
21 GTO 02
22 RCL 02
23 RCL 01
24 RCL 04
25 -
26 XEQ IND 00
27 RCL 03
28 X>Y?
29 GTO 03
30 DSE 05
31 GTO 01
32 LBL 02
33 X<>Y
34 STO 03
35 RCL 04
36 ST+ 01
37 GTO 01
38 LBL 03
39 X<>Y
40 STO 03
41 RCL 04
42 ST- 01
43 LBL 01
44 RCL 02
45 RCL 04
46 +
47 RCL 01
48 XEQ IND 00
49 RCL 03
50 X>Y?
51 GTO 02
52 RCL 02
53 RCL 04
54 -
55 RCL 01
56 XEQ IND 00
57 RCL 03
58 X>Y?
59 GTO 03
60 DSE 05
61 GTO 01
62 LBL 02
63 X<>Y
64 STO 03
65 RCL 04
66 ST+ 02
67 GTO 01
68 LBL 03
69 X<>Y
70 STO 03
71 RCL 04
72 ST- 02
73 LBL 01
74 DSE 05
75 GTO 01
76 PI
77 ST/ 04
78 RCL 04
79 RND
80 X#0?
81 GTO 01
82 RCL 03
83 RCL 02
84 RCL 01
85 END
( 118 bytes / SIZE 006 )
| STACK | INPUTS | OUTPUTS |
| Z | h | min(f) = f(x,y) |
| Y | y0 | y |
| X | x0 | x |
Example: Find the minimum of the
function f(x,y) = | x.y - PI | + | x2 - 2
|
LBL "T" ( any global label, maximum 6 characters
)
ST* Y
X^2
2
-
ABS
X<>Y
PI
-
ABS
+
RTN "T" ASTO 00
and if we choose x0 = y0 = 0 as initial
guesses and h = 1
1 ENTER^
0 ENTER^ XEQ "EXY" >>>> the successive
x-values are displayed and finally,
we have: x = 1.141213652
= R01 ( execution
time = 138 seconds )
RDN y
= 2.221441470 = R02
RDN fmin
= 10 -9
= R03
-The solution was again trivial: x = SQRT(2) , y = PI/SQRT(2) , fmin = 0
Warning: The usual method to find f(x,y) minimum is to solve the system: fx' = 0 ; fy' = 0 ( equating the partial derivatives to zero )
-But the solution of this system is not always a minimum of f
!
-Likewise, "EXY" may converge to any values ( x0 ;
y0 )
such that r(x) = f ( x ; y0
) has a minimum in x0
and
s(y) = f ( x0 ; y ) ----------------- y0
-So don't use this program blindly, especially when "abolute values"
appears in its definition!
3°) N-dimensional Problem
Data Registers:
• R00 = function name
( R00 and R11 to R10+n are to be initialized before executing "EXY"
)
R01 thru R05: unused
R06 & R07: temp
R08 = h ; R09 = n ( number of unknowns )
R10 = f ( x1 ; x2 ; ...... ; xn )
• R11 = x1 • R12 = x2 ; ...........
; • R10+n = xn
Flags: /
Subroutine: A program which computes
f ( x1 ; x2 ; ...... ; xn ) assuming x1
is in R11 , x2 is in R12 , ............ , xn is in
R10+n
01 LBL "EXN"
02 STO 09
03 X<>Y
04 STO 08
05 XEQ IND 00
06 STO 10
07 10.01
08 STO 07
09 LBL 01
10 RCL 09
11 ST+ 07
12 STO 06
13 ISG 06
14 LBL 02
15 VIEW 11
16 RCL 08
17 ST+ IND 07
18 XEQ IND 00
19 RCL 10
20 X>Y?
21 GTO 03
22 RCL 08
23 ST+ X
24 ST- IND 07
25 XEQ IND 00
26 RCL 10
27 X>Y?
28 GTO 03
29 RCL 08
30 ST+ IND 07
31 DSE 06
32 GTO 04
33 LBL 03
34 X<>Y
35 STO 10
36 LBL 04
37 DSE 07
38 GTO 02
39 DSE 06
40 GTO 01
41 PI
42 ST/ 08
43 RCL 08
44 RND
45 X#0?
46 GTO 01
47 RCL 10
48 CLD
49 END
( 81 bytes / SIZE nnn+11 )
| STACK | INPUTS | OUTPUTS |
| Y | h | / |
| X | n | min(f) |
Example: Suppose we want to solve the overdetermined
system: x2 + y2 + z2 = 4
; x2 - y2 + z = 2 ; x + y + z2
= 1 ; x - y - z = 0
by the least squares method ( the system itself has no solution )
-Let f(x,y,z) = ( x2 + y2 + z2 - 4 )2 + ( x2 - y2 + z - 2 )2 + ( x + y + z2 - 1 )2 + ( x - y - z )2
LBL "T"
+
-
X^2
-
+
-
-
RTN
RCL 11
RCL 13 X^2
-
X^2
RCL 13 X^2
RCL 13
X^2
X^2 RCL 11
RCL 13 +
X^2
+
-
RCL 12
+
X^2
+
RCL 11 +
RCL 11 X^2
X^2
4
RCL 12 2
RCL 12 1
RCL 12 +
"T" ASTO 00 CLX STO 11 STO 12 STO 13 if we choose x = y = z = 0 as initial guesses
-Here, n = 3 and if we take again h = 1 and for instance FIX 4
FIX 4
1 ENTER^
3 XEQ "EXN" >>>> the successive
x-values are displayed ( note that 1 is displayed several times ) and
4mn35s later:
min(f) = 1.4344 = X-register = R10 and R11 = x = 1.1055 , R12 = y = -0.9168 , R13 = z = 1.2630
-If you want an extra decimal, simply press:
FIX 5 XEQ 01 >>>> min(f) = 1.43438 , R11 = x = 1.10544 , R12 = y = -0.91683 , R13 = z = 1.26297
-Obviously, the method is quite elementary and ( very ) long execution
time is to be expected
if the function is complicated and/or for large n-values...
Warning: The same limitations described for
"EXY" apply a fortiori to "EXN"
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall