The Museum of HP Calculators

# Extrema for the HP-41

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.

## Overview

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 |

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"