hp41programs

Extrema

Extrema for the HP-41


Overview
 

 1°) One-dimensional Problems

  1-a) Program#1
  1-b) Program#2 - Parabolic Interpolation
  1-c) Discrete Data: 3 tabulated values
  1-d) Discrete Data: 5 tabulated values

 2°) Two-dimensional Problems

  2-a) Program#1
  2-b) Program#2 - Parabolic Interpolation
  2-c) Discrete Data: grid of 9 points
  2-d) Discrete Data: grid of 25 points ( 2 programs )

 3°) Three-dimensional Problems

  3-a) Program#1
  3-b) Discrete Data: grid of 125 points

 4°) N-dimensional Problems

  4-a) Program#1
  4-b) Program#2 - Parabolic Interpolation
 

-These programs seek a local extremum of a function f when the derivative is not available - or is difficult to calculate.
-The parabolic interpolations only work for smooth functions.

>>>  The latest additions are the second program in paragraph 2-d) and paragraph 3°)
 

1°) One-dimensional Problems
 

     1-a) Program#1
 

-This program finds a minimum of a given function f.
-Seek a minimum of (-f)  if you want to know a maximum of the function.
 

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
26  X>Y?
27  X<>Y
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
53  PI
54  ST/ 03
55  RCL 03
56  RND
57  X#0?
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

Note:

-The iteration stops when the rounded h-value equals zero - lines 56-57.
-So, the accuracy is controlled by the display format.
 

     1-b) Program#2 - Parabolic Interpolation
 

-This program uses the parabola  y = a x2 + b x + c   passing through 3 points ( x1 , f(x1) ) , ( x2 , f(x2) ) , ( x3 , f(x3) )
-It calculates the abscissa of the extremum  x4 = -b/(2a)  which is the next approximation

-The algorithm is repeated until 2 consecutive rounded approximations are equal.
 

Data Registers:           •  R00 = function name          ( Register R00 is to be initialized before executing "PEX" )

                                          R01 = x3       R04 = f(x3)        R07 = -a
                                          R02 = x2       R05 = f(x2)
                                          R03 = x1       R06 = f(x1)
Flags: /
Subroutine:   A program that computes f(x) assuming x is in X-register upon entry.
 
 

01  LBL "PEX"
02  STO 01
03  RDN
04  STO 02
05  X<>Y
06  STO 03
07  XEQ IND 00
08  STO 06
09  RCL 02
10  XEQ IND 00
11  STO 05
12  CLX
13  STO 07
14  LBL 01
15  VIEW 01
16  RCL 01
17  XEQ IND 00
18  STO 04
19  RCL 05
20  -
21  RCL 01
22  RCL 02
23  -
24  /
25  RCL 05
26  RCL 06
27  -
28  RCL 02
29  RCL 03
30  -
31  /
32  -
33  RCL 03
34  RCL 01        
35  -
36  /
37  X#0?
38  STO 07
39  RCL 02
40  RCL 01
41  -
42  *
43  +
44  RCL 07
45  ST+ X
46  X#0?
47  /
48  RCL 01        
49  +
50  X<> 01
51  X<> 02
52  STO 03
53  RCL 04
54  X<> 05
55  STO 06
56  RCL 01
57  RND
58  RCL 02 
59  RND
60  X#Y?
61  GTO 01 
62  RCL 07
63  X#0?
64  SIGN
65  RCL 04        
66  RCL 01
67  END

 
   ( 84 bytes / SIZE 008 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             x1          +/-1
           Y             x2           f(x)
           X             x3            x

  if  Z-register = +1  we have a local maximum                           However, roundoff errors can play a trick
  if  Z-register = -1  we have a local minimum                             and Z-output is not completely
  if  Z-register = 0   the conclusion is doubtful                             guaranteed...

Example:            f(x) = ( 1/x5 )/[ exp(1/x) - 1 ]        ( Planck's function )

-We load a routine to compute f(x)

   LBL "T"
   1/X
   E^X-1
   LASTX
   5
   Y^X
   X<>Y
   /
   RTN                 "T"   ASTO 00

-If we choose    0.1  0.2   0.3   as initial guesses and  FIX 9

   0.1   ENTER^
   0.2   ENTER^
   0.3   XEQ "PEX"  the successive x-values are displayed and eventually it yields:

                 x  =  0.201406520       ( execution time = 38 seconds )
    RDN   f(x) = 21.20143565
    RDN    +1                                 so f has a maximum value (  21.20143565 )  for  x  =  0.201406520
 

-In this example, we can of course calculate the first derivative and the true x-value is  0.2014052353  ( f(x) = 21.20143566 )
-As usual in these methods, x cannot be computed very precisely whereas  f(x)  is very accurate! ( often but not always... )
 

     1-c) Discrete Data:  3 Tabulated Values
 

-The following program uses the Lagrange 3-point interpolation formula.
-We assume that the xi's are equally spaced. ( xi+1 - xi = h  ,   i  =  -1 , 0 )
-So,   f(x0+p.h) - f(x0) = A.p + B.p2
 

Data Registers:                                                            ( Registers R01 thru R03 are to be initialized before executing "DEX3" )

                                      •  R01 = f(x-1)    •  R02 = f(x0)    •  R03 = f(x1)                   R00 = h  ,  R04-R05: temp
Flags: /
Subroutines: /
 
 

01  LBL "DEX3"
02  STO 04
03  X<>Y
04  STO 00       
05  RCL 02
06  RCL 01
07  -
08  STO 05       
09  RCL 03 
10  RCL 02
11  -
12  ST- 05
13  +
14  ENTER^
15  X^2
16  RCL 05       
17  8
18  *
19  /
20  RCL 02
21  +
22  X<>Y
23  RCL 05       
24  ST+ X
25  /
26  RCL 00
27  *
28  RCL 04       
29  +
30  END

 
   ( 41 bytes / SIZE 006 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h    extr(f) = f(x)
           X             x0            x

 
Example:     Estimate the minimum of f(x) from these 3 values
 
 

      x       3       5       7
      y       5       3       4

 
  5  STO 01   3  STO 02   4  STO 03

  h = 5-3 = 7-5 = 2  and  x0 = 5

  2  ENTER^
  5  XEQ "DEX3"   >>>>     x   =  5.3333         ( in 1 second )
                              RDN    fmin =  2.9583
 

     1-d) Discrete Data:  5 Tabulated Values
 

-The following program uses the Lagrange 5-point interpolation formula.
-We assume that the xi's are equally spaced. ( xi+1 - xi = h  ,   i  = -2 , -1 , 0 , 1 )
-So,   f(x0+p.h) - f(x0) = A.p + B.p2 + C.p3 + D.p4

Data Registers:                                                                      ( Registers R01 thru R05 are to be initialized before executing "DEX5" )

                                      •  R01 = f(x-2)    •  R02 = f(x-1)    •  R03 = f(x0)    •  R04 = f(x1)    •  R05 = f(x2)               R00 & R06 to R11; temp

Flag:   F01  is cleared
Subroutine:  "P3"  ( Cubic Equations - see for instance "Polynomials for the HP-41" )
 
 

01  LBL "DEX5"
02  STO 10
03  X<>Y
04  STO 00
05  RCL 01       
06  RCL 05 
07  +
08  RCL 02
09  RCL 04
10  +
11  STO 07
12  4
13  *
14  -
15  6
16  /
17  RCL 03
18  +
19  STO 11 
20  4
21  /
22  STO 06       
23  RCL 03
24  +
25  RCL 07
26  2
27  /
28  -
29  STO 08
30  RCL 01
31  RCL 05
32  -
33  RCL 02
34  RCL 04
35  -
36  STO 07       
37  8
38  *
39  -
40  12
41  /
42  STO 09
43  RCL 07
44  2
45  /
46  +
47  CHS
48  STO 07
49  RCL 11 
50  X<>Y
51  3
52  *
53  RCL 08       
54  ST+ X
55  CHS
56  RCL 09
57  XEQ "P3"
58  FC?C 01
59  CLX
60  RCL 07
61  RCL 06
62  RCL Z
63  *
64  +
65  *
66  RCL 08 
67  -
68  *
69  RCL 09       
70  +
71  *
72  RCL 03
73  +
74  X<>Y
75  RCL 00
76  *
77  RCL 10
78  +
79  END

 
  ( 95 bytes / SIZE 012 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h    extr(f) = f(x)
           X             x0            x

 
Example:     Estimate the minimum of f(x) from these 5 values
 
 

      x       1       3       5       7      9
      y       9       5       3       4      7

 
  9  STO 01   5  STO 02   3  STO 03   4  STO 04   7  STO 05

  h = 3-1 = 5-3 = 7-5 = 9-7 = 2  and  x0 = 5

  2  ENTER^
  5  XEQ "DEX5"  >>>>   x   =  5.316624790       ( in 6 seconds )
                             RDN  fmin =  2.960474247

Remark:

-If you don't want to use a subroutine like "P3" , the problem may be solved by the iterative method of the following variant.
 ( but the method fails if C = 0 )
-The accuracy is controlled by the display settings.
 
 

01  LBL "DEX5"
02  STO 10
03  X<>Y
04  STO 00
05  RCL 01       
06  RCL 05
07  +
08  RCL 02
09  RCL 04
10  +
11  STO 07 
12  4
13  *
14  -
15  6
16  /
17  RCL 03
18  +
19  STO 11
20  4
21  /
22  STO 06 
23  RCL 03
24  +
25  RCL 07       
26  2
27  /
28  -
29  STO 08 
30  RCL 01
31  RCL 05
32  -
33  RCL 02
34  RCL 04
35  -
36  STO 07
37  8
38  *
39  -
40  12
41  /
42  STO 09 
43  RCL 07       
44  2
45  /
46  +
47  CHS
48  STO 07 
49  CLST
50  LBL 01
51  CLX
52  RCL 11
53  *
54  RCL 07
55  3
56  *
57  +
58  *
59  *
60  RCL 09 
61  +
62  RCL 08       
63  ST+ X
64  /
65  ENTER^
66  ENTER^
67  R^
68  -
69  RND
70  X#0?
71  GTO 01
72  CLX
73  RCL 07 
74  RCL 06
75  RCL Z
76  *
77  +
78  *
79  RCL 08       
80  -
81  *
82  RCL 09
83  +
84  *
85  RCL 03
86  +
87  X<>Y
88  RCL 00
89  *
90  RCL 10 
91  +
92  END

 
  ( 105 bytes / SIZE 012 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h    extr(f) = f(x)
           X             x0            x

 
Example:     Estimate the minimum of f(x) from these 5 values
 
 

      x       1       3       5       7      9
      y       9       5       3       4      7

 
  9  STO 01   5  STO 02   3  STO 03   4  STO 04   7  STO 05

  h = 3-1 = 5-3 = 7-5 = 9-7 = 2  and  x0 = 5

      FIX 9
  2  ENTER^
  5  XEQ "DEX5"  >>>>   x   =  5.316624790       ( in 6 seconds )
                             RDN  fmin =  2.960474247

-This method is used in "TIDE1"  ( cf "Tides for the HP-41" )
 

2°) Two-dimensional Problems
 

     2-a) Program#1
 

-This program finds a minimum of a given function f.
-Seek a minimum of (-f)  if you want to know a maximum of the function.
 

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:     To find f(x,y) minimum, we usually 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 for x = x0
         and    s(y) =  f ( x0 , y )  has a minimum for y = y0
-So don't use this program blindly, especially if "absolute values" appear in its definition!
 

     2-b) Program#2 - Parabolic Interpolation
 

-The method used in paragraph 1-b) is now generalized to a function of 2 variables.
 

Data Registers:           •  R00 = function name          ( Register R00 is to be initialized before executing "PEXY" )

                                          R01 = x3        R03 = x2          R05 =  x1        R07 = f(x,y)            R10 =  -ax
                                          R02 = y3        R04 = y2          R06 =  y1       R08-R09: temp        R11 =  -ay
Flags: /
Subroutine:   A program that computes  f(x,y) assuming x is in X-register and y is in Y-register upon entry.

-The accuracy is controlled by the display setting.
-Lines 119 & 125 are three-byte GTO 01
 
 

  01  LBL "PEXY"
  02  STO 01
  03  STO 03
  04  STO 05
  05  RDN
  06  STO 02
  07  STO 04
  08  STO 06
  09  X<>Y
  10  ST+ 01
  11  ST+ 02
  12  ST- 05
  13  ST- 06
  14  LBL 01
  15  VIEW 01
  16  RCL 02
  17  RCL 05
  18  XEQ IND 00
  19  STO 09
  20  RCL 02
  21  RCL 03
  22  XEQ IND 00
  23  STO 08
  24  RCL 02
  25  RCL 01
  26  XEQ IND 00
  27  STO 07
  28  RCL 08
  29  -
  30  RCL 01
  31  RCL 03
  32  -
  33  X#0?
  34  /
  35  RCL 08
  36  RCL 09
  37  -
  38  RCL 03
  39  RCL 05
  40  -
  41  X#0?
  42  /
  43  -
  44  RCL 05 
  45  RCL 01 
  46  -
  47  X#0?
  48  /
  49  X#0?
  50  STO 10
  51  RCL 03        
  52  RCL 01
  53  -
  54  *
  55  +
  56  RCL 10
  57  ST+ X
  58  X#0?
  59  /
  60  RCL 01
  61  +
  62  X<> 01
  63  X<> 03
  64  STO 05
  65  RCL 06
  66  RCL 01
  67  XEQ IND 00
  68  STO 09
  69  RCL 04
  70  RCL 01
  71  XEQ IND 00
  72  STO 08
  73  RCL 02
  74  RCL 01
  75  XEQ IND 00
  76  STO 07
  77  RCL 08
  78  -
  79  RCL 02
  80  RCL 04
  81  -
  82  X#0?
  83  /
  84  RCL 08
  85  RCL 09
  86  -
  87  RCL 04
  88  RCL 06
  89  -
  90  X#0?
  91  /
  92  -
  93  RCL 06
  94  RCL 02
  95  -
  96  X#0?
  97  /
  98  X#0?
  99  STO 11 
100  RCL 04        
101  RCL 02
102  -
103  *
104  +
105  RCL 11
106  ST+ X
107  X#0?
108  /
109  RCL 02
110  +
111  X<> 02
112  X<> 04
113  STO 06
114  RCL 01
115  RND
116  RCL 03
117  RND
118  X#Y?
119  GTO 01
120  RCL 02
121  RND
122  RCL 04
123  RND
124  X#Y?
125  GTO 01
126  RCL 10 
127  RCL 11        
128  *
129  X<0?
130  CLX
131  X#0?
132  LASTX
133  X#0?
134  SIGN
135  RCL 07
136  RCL 02
137  RCL 01
138  END

 
   ( 168 bytes / SIZE 012 )
 
 

      STACK        INPUTS      OUTPUTS
           T             /         +/-1
           Z             h         f(x,y)
           Y             y0            y
           X             x0            x

  if  T-register = +1  we have a local maximum                           However - due to roundoff errors -
  if  T-register = -1  we have a local minimum                             these conclusions are not completely
  if  T-register = 0  the conclusion is doubtful                              guaranteed...

Example:        f(x,y) = exp [ 2 - ( x - sqrt(2) )2 - ( y - sqrt(3) )2 ]

  LBL "T"         X<>Y         +                RTN
  2                    3                 2
  SQRT            SQRT         X<>Y
  -                    -                  -
  X^2               X^2            E^X

  "T"  ASTO 00    if we choose  h = 1   x0 = y0 = 2    and  FIX 5

   1  ENTER^
   2  ENTER^  XEQ "PEXY"  the successive x-values are displayed and finally ( in 75 seconds ),

   it gives     x    = 1.414213686             the solutions are of course  x   =  sqrt(2)
    RDN      y    = 1.732055058                                                      y   =  sqrt(3)
    RDN   f(x,y) = 7.389056099             all the digits are right for f(x,y) =  exp(2)
    RDN     T    = +1  ( maximum )
 

     2-c) Discrete Data:  Grid of 9 Points
 
 

-The following program uses the Lagrange 3-point interpolation formula in the x-direction and in the y-direction
-We assume that the xi's are equally spaced and that the yj's are equally spaced. ( xi+1 - xi = h  & yj+1 - yj = k    i , j =  -1 , 0  )

-So,  f(x0+p.h,y0+q.k) - f(x0,y0) = A.p + B.q + C.p2 + D.p.q + E.q2 + F.p2q + G.p.q2 + H.p2q2

-The 8 coefficients A , B , C , D , E , F , G , H  are obtained by the formulae:

       2.A = f10 - f-10                      4.D =  f-1-1 - f1-1 - f-11 + f11           F = D - B - ( f-1-1 - f-11 )/2
       2.B = f01 - f0-1                      2.E =   f0-1 - 2 f00 + f01                  G = D - A - ( f-1-1 + f-11 )/2                             with   fi j = f(x0+i.h,y0+j.k)
       2.C = f-10 - 2 f00 + f10                                                                H = - C - E - f00 + D + ( f-11 + f1-1 )/2

-Equating to zero the partial derivatives with respect to p and q leads to a non-linear system that is re-written:

       p = ( - A - D.q - G.q2 )/( 2.C + 2.F.q + 2.H.q2 )                               and this system is solved
       q = ( - B - D.p - F.p2)/( 2.E + 2.G.p + 2.H.p2 )                                by a successive approximation method.
 

N.B:   The program fails if the denominator(s) equal(s) zero !
 

Data Registers:                                                  ( Registers R01 thru R09 are to be initialized before executing "DEXY9" )

                                      •  R01 = f(x-1,y-1)    •  R04 = f(x-1,y0)    •  R07 = f(x-1,y1)
                                      •  R02 = f(x0,y-1)     •  R05 = f(x0,y0)     •  R08 = f(x0,y1)           R00 & R11 to R23: temp   ( R11 = p , R12 = q )
                                      •  R03 = f(x1,y-1)     •  R06 = f(x1,y0)     •  R09 = f(x1,y1)
Flags: /
Subroutines: /

-The accuracy is controlled by the display format.
 
 

 01  LBL "DEXY9"
  02  STO 19
  03  RDN
  04  STO 20        
  05  RDN
  06  STO 21
  07  X<>Y
  08  STO 22
  09  RCL 04
  10  RCL 06
  11  +
  12  RCL 05
  13  ST+ X
  14  STO 10
  15  -
  16  STO 13
  17  RCL 02
  18  LASTX
  19  -
  20  RCL 08
  21  +
  22  STO 14
  23  RCL 01
  24  RCL 03
  25  -
  26  RCL 07
  27  -
  28  RCL 09
  29  +
  30  2
  31  /
  32  STO 15
  33  ST- 10
  34  RCL 08
  35  RCL 02
  36  -
  37  STO 18        
  38  -
  39  RCL 01 
  40  -
  41  RCL 07 
  42  +
  43  STO 16
  44  RCL 15
  45  RCL 06
  46  RCL 04
  47  -
  48  STO 00
  49  -
  50  RCL 01
  51  -
  52  RCL 03
  53  +
  54  STO 17
  55  RCL 13
  56  RCL 14
  57  +
  58  RCL 03
  59  -
  60  RCL 07
  61  -
  62  ST+ 10
  63  CLX
  64  STO 11
  65  STO 12
  66  LBL 01
  67  VIEW 11
  68  RCL 12
  69  RCL 17        
  70  *
  71  RCL 15 
  72  +
  73  RCL 12 
  74  *
  75  RCL 00
  76  +
  77  RCL 10
  78  RCL 12
  79  *
  80  RCL 16
  81  -
  82  RCL 12
  83  *
  84  RCL 13
  85  -
  86  ST+ X
  87  /
  88  ENTER^
  89  X<> 11
  90  -
  91  ABS
  92  RCL 11
  93  RCL 16
  94  *
  95  RCL 15
  96  +
  97  RCL 11
  98  *
  99  RCL 18
100  +
101  STO 23        
102  RCL 10
103  RCL 11 
104  *
105  RCL 17 
106  -
107  RCL 11
108  *
109  RCL 14
110  -
111  ST+ X
112  /
113  ENTER^
114  X<> 12
115  -
116  ABS
117  +
118  RND
119  X#0?
120  GTO 01
121  RCL 17
122  RCL 10
123  RCL 11
124  *
125  -
126  RCL 11
127  *
128  RCL 14
129  +
130  RCL 12
131  *
132  RCL 23        
133  +
134  RCL 12 
135  *
136  RCL 11 
137  RCL 13
138  *
139  RCL 00
140  +
141  RCL 11
142  *
143  +
144  2
145  /
146  RCL 05
147  +
148  RCL 12
149  RCL 22
150  *
151  RCL 20
152  +
153  RCL 11
154  RCL 21
155  *
156  RCL 19
157  +
158  CLD
159  END

 
   ( 197 bytes / SIZE 024 )
 
 

      STACK        INPUTS      OUTPUTS
           T             k         f(x,y)
           Z             h         f(x,y)
           Y             y0            y
           X             x0            x

        f(x,y) = f-extremum

Example:     Estimate the minimum of f(x,y) from these 9 values
 
 

    x \ y       1       4       7
      1       7       4       6
      3       6       3       4
      5       8       6       9

 
                     7     4    6                        R01   R04   R07
    Store         6     3    4         into         R02   R05   R08         respectively
                     8     6    9                        R03   R06   R09

 We have    x0 = 3    h = 3-1 = 5-3 = 2
                  y0 = 4    k = 4-1 = 7-4 = 3   and if we choose  FIX 9

       FIX 9
   3  ENTER^
   2  ENTER^
   4  ENTER^
   3  XEQ "DEXY9"  >>>>  the successive p-values are displayed and finally

                             x    = 2.507454957                       ( execution time = 12 seconds )
             RDN        y    = 4.784962554
             RDN     f(x,y) = 2.736025820  = f-minimum
 

     2-d) Discrete Data:  Grid of 25 Points
 

-The following program uses the Lagrange 5-point interpolation formula in the x-direction and in the y-direction
-We assume that the xi's are equally spaced and that the yj's are equally spaced. ( xi+1 - xi = h  & yj+1 - yj = k    i , j = -2 , -1 , 0 , 1 )

-We have then

  576 f(x0+p.h,y0+q.k)   =   (q2-1)q(q-2) [ (p2-1)p(p-2) f-2-2 - 4.(p-1)p(p2-4) f-1-2 + 6.(p2-1)(p2-4) f0-2 - 4.(p+1)p(p2-4) f1-2 + (p2-1)p(p+2) f-2-2 ]
                                     - 4.(q-1)q(q2-4) [ (p2-1)p(p-2) f-2-1 - 4.(p-1)p(p2-4) f-1-1 + 6.(p2-1)(p2-4) f0-1 - 4.(p+1)p(p2-4) f1-1 + (p2-1)p(p+2) f-2-1 ]
                                     + 6.(q2-1)(q2-4) [ (p2-1)p(p-2) f-20 - 4.(p-1)p(p2-4) f-10 + 6.(p2-1)(p2-4) f00 - 4.(p+1)p(p2-4) f10 + (p2-1)p(p+2) f-20 ]
                                    - 4.(q+1)q(q2-4) [ (p2-1)p(p-2) f-21 - 4.(p-1)p(p2-4) f-11 + 6.(p2-1)(p2-4) f01 - 4.(p+1)p(p2-4) f11 + (p2-1)p(p+2) f-21 ]
                                      +  (q2-1)q(q+2)  [ (p2-1)p(p-2) f-22 - 4.(p-1)p(p2-4) f-12 + 6.(p2-1)(p2-4) f02 - 4.(p+1)p(p2-4) f12 + (p2-1)p(p+2) f-22 ]

            with   fi j = f(x0+i.h,y0+j.k)

-We seek an extremum of the function by equating to zero the partial derivatives with respect to p and q
-This system of 2 non-linear equations is re-written  p = F(p,q) ; q = G(p,q)
  and a successive approximation method is used.
 

N.B:

-The routine fails if the denominator(s) of  F or G  equal(s) zero!
-We could also calculate directly the 24 coefficients of the above polynomial ( like "DEXY9" does ), but the program would be huge!
-That's why "DEXY25" is slow...
 

Data Registers:                         ( Registers R01 thru R25 are to be initialized before executing "DEXY25" )

                                      •  R01 = f(x-2,y-2)    •  R06 = f(x-2,y-1)    •  R11 = f(x-2,y0)    •  R16 = f(x-2,y1)    •  R21 = f(x-2,y2)
                                      •  R02 = f(x-1,y-2)    •  R07 = f(x-1,y-1)    •  R12 = f(x-1,y0)    •  R17 = f(x-1,y1)    •  R22 = f(x-1,y2)
                                      •  R03 = f(x0,y-2)     •  R08 = f(x0,y-1)     •  R13 = f(x0,y0)     •  R18 = f(x0,y1)     •  R23 = f(x0,y2)
                                      •  R04 = f(x1,y-2)     •  R09 = f(x1,y-1)     •  R14 = f(x1,y0)     •  R19 = f(x1,y1)     •  R24 = f(x1,y2)
                                      •  R05 = f(x2,y-2)     •  R10 = f(x2,y-1)     •  R15 = f(x2,y0)     •  R20 = f(x2,y1)     •  R25 = f(x2,y2)

Flag:   F10  is cleared                                      R00 = p ,  R26 = q        R27 to R57: temp
Subroutines: /

-Line 41 is a three-byte GTO 04 - not absolutely necessary since this line is executed only once !
-The accuracy is controlled by the display format.
-If you don't have an X-Functions module, replace lines 204-205 by
  RCL 41  STO 46  RCL 42  STO 47  RCL 43  STO 48  RCL 44  STO 49  RCL 45  STO 50
 
 

01  LBL "DEXY25"
  02  STO 27         
  03  RDN
  04  STO 28
  05  RDN
  06  STO 29
  07  X<>Y
  08  STO 30
  09  2.90005
  10  STO 41
  11  2
  12  +
  13  STO 42
  14  1
  15  +
  16  STO 43
  17  4
  18  -
  19  STO 44
  20  31.035
  21  STO 45
  22  XEQ 00
  23  1.9
  24  STO 44
  25  5
  26  +
  27  STO 41
  28  10
  29  +
  30  STO 42
  31  5
  32  +
  33  STO 43
  34  .005
  35  ST+ 45
  36  XEQ 00
  37  CLX
  38  STO 00
  39  STO 26
  40  SF 10
  41  GTO 04
  42  LBL 00
  43  RCL IND 41
  44  RCL IND 42
  45  -
  46  8
  47  *
  48  RCL IND 43
  49  +
  50  RCL IND 44 
  51  -
  52  STO IND 45
  53  ISG 41  
  54  ISG 42
  55  ISG 43
  56  ISG 44
  57  ISG 45
  58  GTO 00
  59  RTN
  60  LBL 01
  61  RCL 26          
  62  ENTER^
  63  ENTER^
  64  X^2
  65  1
  66  -
  67  STO 43
  68  *
  69  STO 41
  70  STO 45
  71  CLX
  72  SIGN
  73  ST+ Z
  74  -
  75  R^
  76  ST* Z
  77  *
  78  STO 42
  79  CLX
  80  +
  81  STO 44
  82  CLX
  83  2
  84  ST+ Z
  85  -
  86  ST* 41
  87  X<>Y
  88  ST* 45 
  89  *
  90  6
  91  *
  92  ST* 43 
  93  1.5
  94  CHS
  95  /
  96  ST* 42
  97  ST* 44
  98  1.9
  99  STO 51          
100  FC? 10
101  RTN
102  RCL 00
103  X^2
104  ST+ X
105  1
106  -
107  STO 46
108  STO 50
109  4
110  -
111  6
112  *
113  STO 48
114  16
115  RCL 00
116  X^2
117  8
118  *
119  -
120  STO 47
121  RCL 00
122  3
123  *
124  ST- 46
125  ST+ 50
126  ST+ X
127  ST+ 47
128  -
129  STO 49          
130  RTN
131  LBL 02
132  RCL IND 51 
133  RCL IND 52
134  *
135  +
136  ISG 51
137  ISG 52
138  GTO 02
139  RCL IND 53
140  *
141  +
142  5
143  ST- 52
144  CLX
145  RCL 54
146  ST- 51
147  CLX
148  ISG 53
149  GTO 02
150  FC? 10
151  RTN
152  LBL 03
153  RCL IND 55
154  RCL IND 56
155  *
156  +
157  ISG 55
158  ISG 56
159  GTO 03
160  X<>Y
161  /
162  ENTER^
163  X<> 00
164  -
165  ABS
166  RTN
167  LBL 04 
168  VIEW 00 
169  XEQ 01
170  46.05
171  STO 52          
172  41.045
173  STO 53
174  STO 56
175  31.9
176  STO 55
177  CLST
178  STO 54
179  XEQ 02
180  STO 57
181  RCL 00
182  X<> 26
183  STO 00
184  XEQ 01
185  5
186  ST- 53
187  ST- 56
188   E5
189  /
190  ST+ 51
191  24
192  STO 54
193  CLST
194  XEQ 02
195  RCL 00
196  X<> 26
197  STO 00
198  X<> 57
199  X<Y?
200  X<>Y
201  RND
202  X#0?
203  GTO 04
204  41.046005
205  REGMOVE 
206  CF 10
207  XEQ 01
208  5
209  ST- 53
210  CLST
211  STO 54          
212  XEQ 02
213  CLX
214  576
215  /
216  RCL 26
217  RCL 30
218  *
219  RCL 28
220  +
221  RCL 00
222  RCL 29
223  *
224  RCL 27
225  +
226  CLD
227  END

 
   ( 397 bytes / SIZE 058 )
 
 

      STACK        INPUTS      OUTPUTS
           T             k         f(x,y)
           Z             h         f(x,y)
           Y             y0            y
           X             x0            x

        f(x,y) = f-extremum

Example:     Estimate the minimum of f(x,y) from these 25 values
 
 

    x \ y     -2       1       4       7     10
     -1     12       9       6       9     13
      1      9       7       4       6      9
      3      9       6       3       4      7
      5     10       8       6       9     11
      7     14      11       8      12     14

 
                12   9     6    9    13                    R01   R06   R11   R16   R21
                 9    7     4    6     9                     R02   R07   R12   R17   R22
    Store     9    6     3    4     7       into        R03   R08   R13   R18   R23       respectively
                10   8     6    9    11                    R04   R09   R14   R19   R24
                14  11    8   12   14                    R05   R10   R15   R20   R25

 We have    x0 = 3      h = 1-(-1) = 3-1 = 5-3 = 7-5  = 2
    and         y0 = 4      k = 1-(-2) = 4-1 = 7-4 = 10-7 = 3     and if we choose  FIX 9

       FIX 9
   3  ENTER^
   2  ENTER^
   4  ENTER^
   3  XEQ "DEXY25"  >>>>  the successive p-values are displayed and finally

                                      x    = 2.508078775                                ( execution time = 3mn38s )
                      RDN        y    = 4.813612947
                      RDN     f(x,y) = 2.682512101  =  f-minimum

Note:

-The 2nd program listed hereunder saves 56 bytes and is slightly faster
-It uses quite a similar method but less data registers are employed.
 

Data Registers:                                           ( Registers R01 thru R25 are to be initialized before executing "EXY25" )

                                      •  R01 = f(x-2,y-2)    •  R06 = f(x-2,y-1)    •  R11 = f(x-2,y0)    •  R16 = f(x-2,y1)    •  R21 = f(x-2,y2)
                                      •  R02 = f(x-1,y-2)    •  R07 = f(x-1,y-1)    •  R12 = f(x-1,y0)    •  R17 = f(x-1,y1)    •  R22 = f(x-1,y2)
                                      •  R03 = f(x0,y-2)     •  R08 = f(x0,y-1)     •  R13 = f(x0,y0)     •  R18 = f(x0,y1)     •  R23 = f(x0,y2)
                                      •  R04 = f(x1,y-2)     •  R09 = f(x1,y-1)     •  R14 = f(x1,y0)     •  R19 = f(x1,y1)     •  R24 = f(x1,y2)
                                      •  R05 = f(x2,y-2)     •  R10 = f(x2,y-1)     •  R15 = f(x2,y0)     •  R20 = f(x2,y1)     •  R25 = f(x2,y2)

Flag:   F01  is cleared                                      R00 = p ,  R26 = q        R27 to R50: temp
Subroutines: /

-Lines 154 and 160 are three-byte GTO's
-The accuracy is controlled by the display format ( line 158 ).
 
 

  01  LBL "EXY25"
  02  STO 27
  03  RDN
  04  STO 28 
  05  RDN
  06  STO 29
  07  X<>Y 
  08  STO 30
  09  CLX
  10  STO 00
  11  STO 26         
  12  STO 31
  13  GTO 10
  14  LBL 01
  15  ENTER
  16  ENTER
  17  X^2
  18  1
  19  -
  20  STO 36
  21  *
  22  STO 34
  23  STO 38
  24  CLX
  25  SIGN
  26  ST+ Z
  27  -
  28  R^
  29  ST* Z
  30  *
  31  STO 35
  32  CLX
  33  +
  34  STO 37         
  35  CLX
  36  2
  37  ST+ Z
  38  -
  39  ST* 34
  40  X<>Y
  41  ST* 38
  42  *
  43  6
  44  *
  45  ST* 36
  46  1.5
  47  /
  48  CHS
  49  ST* 35
  50  ST* 37
  51  RTN
  52  LBL 02
  53  ENTER
  54  X^2
  55  STO T
  56  *
  57  ST+ X
  58  STO Y
  59  6
  60  *
  61  STO 41
  62  X<>Y
  63  ENTER
  64  STO 39         
  65  R^
  66  3
  67  *
  68  1
  69  -
  70  ST- 39
  71  +
  72  STO 43
  73  LASTX
  74  ST+ X
  75  6
  76  -
  77  STO 40
  78  CHS
  79  R^
  80  4
  81  *
  82  ST- 40
  83  -
  84  STO 42
  85  RTN
  86  LBL 10
  87  SF 01
  88  LBL 11
  89  RCL 00 
  90  XEQ 01
  91  RCL 26
  92  XEQ 02
  93  5 E-5
  94  1.9
  95  FC? 01
  96  +
  97  STO 33         
  98  34.038
  99  STO 44
100  39.043
101  STO 45
102  46.9
103  STO 32         
104  CLST
105  LBL 12
106  RCL IND 33
107  RCL IND 44
108  *
109  +
110  ISG 33
111  ISG 44
112  GTO 12
113  STO IND 32
114  5
115  ST- 44
116  CLX
117  RCL IND 45
118  *
119  +
120  24
121  FC? 01
122  ST- 33
123  CLX
124  ISG 32
125  ISG 45
126  GTO 12
127  CLX
128  RCL 46 
129  RCL 50         
130  +
131  RCL 47
132  RCL 49
133  +
134  16
135  *
136  -
137  RCL 48
138  30
139  *
140  +
141  /
142  ENTER
143  X<> 26
144  -
145  ABS
146  RCL 31
147  X<Y?
148  X<>Y
149  STO 31
150  RCL 00
151  X<> 26
152  STO 00 
153  FS?C 01
154  GTO 11
155  VIEW 00
156  CLX
157  X<> 31
158  RND
159  X#0?
160  GTO 10
161  RCL 26         
162  XEQ 01
163  34.039005
164  REGMOVE
165  RCL 00
166  XEQ 01
167  5
168  ST- 45
169  SIGN
170  STO 33
171  CLST
172  LBL 13
173  RCL IND 33
174  RCL IND 44
175  *
176  +
177  ISG 33
178  CLX
179  ISG 44
180  GTO 13
181  RCL IND 45
182  *
183  +
184  5
185  ST- 44
186  CLX
187  ISG 45
188  GTO 13
189  CLX
190  576
191  /
192  RCL 26         
193  RCL 30
194  *
195  RCL 28
196  +
197  RCL 00
198  RCL 29
199  *
200  RCL 27
201  +
202  CLD
203  END

 
       ( 340 bytes / SIZE 051 )
 
 

      STACK        INPUTS      OUTPUTS
           T             k         f(x,y)
           Z             h         f(x,y)
           Y             y0            y
           X             x0            x

        f(x,y) = f-extremum

Example:     The same one: Estimate the minimum of f(x,y) from these 25 values
 
 

    x \ y     -2       1       4       7     10
     -1     12       9       6       9     13
      1      9       7       4       6      9
      3      9       6       3       4      7
      5     10       8       6       9     11
      7     14      11       8      12     14

 
                12   9     6    9    13                    R01   R06   R11   R16   R21
                 9    7     4    6     9                     R02   R07   R12   R17   R22
    Store     9    6     3    4     7       into        R03   R08   R13   R18   R23       respectively
                10   8     6    9    11                    R04   R09   R14   R19   R24
                14  11    8   12   14                    R05   R10   R15   R20   R25

 We have    x0 = 3      h = 1-(-1) = 3-1 = 5-3 = 7-5  = 2
    and         y0 = 4      k = 1-(-2) = 4-1 = 7-4 = 10-7 = 3     and if we choose  FIX 9

       FIX 8
   3  ENTER^
   2  ENTER^
   4  ENTER^
   3  XEQ "EXY25"  >>>>  the successive p-values are displayed and finally

                                      x    = 2.508078774                                ( execution time = 3mn21s )
                      RDN        y    = 4.813612946
                      RDN     f(x,y) = 2.682512101  =  f-minimum
 

Notes:

-Each iteration lasts about 31 seconds.
-You can reduce execution time if you store the different control numbers ( 34.038  39.043 ... ) in data registers ( R51  R52 ... ) in the first part of the routine.
-A few bytes may be saved if you replace lines 146 to 149 by a single instruction:  ST+ 31
 

3°) Three-dimensional Problems
 

     3-a) Program#1
 

-Like "EXY" listed in paragraph 2-a), "EXYZ"  finds a minimum of a given function f.
-Seek a minimum of (-f)  if you want to know a maximum of the function.

-In fact, "EXN" - listed in §4-a) - does the same job !
 

Data Registers:           •  R00 = function name        ( This register is to be initialized before executing "EXYZ" )

                                          R01 = x ;  R02 = y  ;  R03 = z  ; R04 = f(x) ;  R05 = h  ;  R06: temp
Flags: /
Subroutine:   A program which computes f(x,y,z) assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.

-Lines 113-119 are three-byte GTO 01
 
 

  01  LBL "EXYZ"
  02  R^
  03  STO 05
  04  R^
  05  STO 03 
  06  R^
  07  STO 02
  08  R^
  09  STO 01
  10  XEQ IND 00
  11  STO 04
  12  LBL 01 
  13  VIEW 01
  14  4 
  15  STO 06
  16  RCL 03
  17  RCL 02
  18  RCL 01
  19  RCL 05
  20  +
  21  XEQ IND 00
  22  RCL 04
  23  X>Y?
  24  GTO 02
  25  RCL 03
  26  RCL 02
  27  RCL 01
  28  RCL 05
  29  -
  30  XEQ IND 00
  31  RCL 04
  32  X>Y?
  33  GTO 03
  34  DSE 06
  35  GTO 01
  36  LBL 02
  37  X<>Y
  38  STO 04
  39  RCL 05
  40  ST+ 01
  41  GTO 01
  42  LBL 03
  43  X<>Y
  44  STO 04
  45  RCL 05
  46  ST- 01
  47  LBL 01
  48  RCL 03
  49  RCL 02
  50  RCL 05
  51  +
  52  RCL 01
  53  XEQ IND 00
  54  RCL 04
  55  X>Y?
  56  GTO 02
  57  RCL 03
  58  RCL 02
  59  RCL 05
  60  -
  61  RCL 01
  62  XEQ IND 00
  63  RCL 04
  64  X>Y?
  65  GTO 03
  66  DSE 06
  67  GTO 01
  68  LBL 02
  69  X<>Y
  70  STO 04
  71  RCL 05
  72  ST+ 02
  73  GTO 01
  74  LBL 03
  75  X<>Y
  76  STO 04
  77  RCL 05
  78  ST- 02
  79  LBL 01
  80  RCL 03
  81  RCL 05
  82  +
  83  RCL 02
  84  RCL 01
  85  XEQ IND 00
  86  RCL 04
  87  X>Y?
  88  GTO 02
  89  RCL 03
  90  RCL 05
  91  -
  92  RCL 02
  93  RCL 01
  94  XEQ IND 00
  95  RCL 04
  96  X>Y?
  97  GTO 03
  98  DSE 06
  99  GTO 01
100  LBL 02
101  X<>Y
102  STO 04
103  RCL 05
104  ST+ 03
105  GTO 01
106  LBL 03
107  X<>Y
108  STO 04
109  RCL 05
110  ST- 03
111  LBL 01
112  DSE 06
113  GTO 01
114  PI
115  ST/ 05
116  RCL 05
117  RND
118  X#0?
119  GTO 01
120  RCL 04         
121  RCL 03
122  RCL 02
123  RCL 01
124  END

 
           ( 168 bytes / SIZE 007 )
 
 

      STACK        INPUTS      OUTPUTS
           T             h  min(f) = f(x,y,z)
           Z             z0            z
           Y             y0            y
           X             x0            x

 
Example:  Find the minimum of the function defined by   f(x,y,z) = ( x2 + y2 + z2 - 4 )2  + ( x2 - y2 + z - 2 )2 + ( x + y + z2 - 1 )2 + ( x - y - z )2

-Load the short routine:
 
 

 01  LBL "T"
 02  STO 07
 03  X^2
 04  X<>Y
 05  STO 08
 06  X^2
 07  +
 08  X<>Y
 09  STO 09
 10  X^2
 11  +
 12  4
 13  -
 14  X^2
 15  RCL 07
 16  X^2
 17  RCL 08
 18  X^2
 19  -
 20  RCL 09
 21  +
 22  2
 23  -
 24  X^2
 25  +
 26  RCL 07
 27  RCL 08
 28  +
 29  RCL 09
 30  X^2
 31  +
 32  1
 33  -
 34  X^2
 35  +
 36  RCL 07
 37  RCL 08
 38  -
 39  RCL 09
 40  -
 41  X^2
 42  +
 43  RTN
 44  END

-If you choose  x0 = y0 = z0 = h = 1   and  FIX 4

    FIX 4
       1     ENTER^  ENTER^  ENTER^  XEQ "EXYZ"   the HP-41 displays the successive x-approximations and returns eventually:

                  x =  1.1054
                  y = -0.9168
                  z =   2.2630   and  fmin = 1.4344                      ---Execution time = 4m32s---

-In FIX 5 , you'll get

                  x =  1.10542
                  y = -0.91684
                  z =   2.26299   and  fmin = 1.43438

-And with FIX 9:

                  x =  1.105422397
                  y = -0.916837262
                  z =   2.262990778   and  fmin = 1.434376612
 

     3-b) Discrete Data: grid of 125 points
 

-Here, we also use the Lagrange polynomials of degree 4 in 3 directions x , y , z
-We assume that the xi's , yj's & zm's  are equally spaced. ( xi+1 - xi = h  , yj+1 - yj = k  ,  ym+1 - ym = l  )         (  i , j , m = -2 , -1 , 0 , 1 )

   f( x0+p.h , y0+q.k , z0+r.l )  is a function in 3 variables in p , q , r

-We seek an extremum of the function by equating to zero the partial derivatives with respect to p , q and r
-This system of 3 non-linear equations is re-written  p = F(p,q,r) ; q = G(p,q,r) ; r = H(p,q,r)
  and a successive approximation method is used.
 

N.B:

-The routine fails if the denominator(s) of  F , G or H  equal(s) zero !
-This program is very slow and an emulator in turbo mode is highly recommended.
 
 

Data Registers:              R00 thru R36: temp      ( Registers Rbb thru Rbb+124 are to be initialized before executing "EXYZ125" )

                                         R01 = p ,  R02 = q  ,  R03 = r           and   with   bb > 36

                                      •  Rbb = f(x-2,y-2,z-2)    •  Rb+5 = f(x-2,y-1,z-2)    •  Rb+10 = f(x-2,y0,z-2)    •  Rb+15 = f(x-2,y1,z-2)    •  Rb+20 = f(x-2,y2,z-2)
                                      •  Rb+1 = f(x-1,y-2,z-2)    ................................................................................................................    •  Rb+21 = f(x-1,y2,z-2)
                                      •  Rb+2 = f(x0,y-2,z-2)      ...............................................................................................................    •  Rb+22 = f(x0,y2,z-2)
                                      •  Rb+3 = f(x1,y-2,z-2)     ................................................................................................................    •  Rb+23 = f(x1,y2,z-2)
                                      •  Rb+4 = f(x2,y-2,z-2)     .................................................................................................................   •  Rb+24 = f(x2,y2,z-2)

                                                    •  Rb+25 = f(x-2,y-2,z-1)
                                                    •  Rb+26 = f(x-1,y-2,z-1)

                                                                                                        .............. and so on  ...............                                                  •  Rb+124 = f(x2,y2,z2)

                                           with   fp q r = f( x0+p.h , y0+q.k , z0+r.l )
 

Flags:   F01 & F02  are cleared
Subroutines: /
 

-Lines 197-199-205 are three-byte GTO's
 
 

01  LBL "EXYZ125"
  02  "H^K^L=?"
  03  PROMPT
  04  STO 36
  05  RDN
  06  STO 35
  07  X<>Y
  08  STO 34
  09  "X^Y^Z=?"
  10  PROMPT
  11  STO 33
  12  RDN
  13  STO 32
  14  X<>Y
  15  STO 31
  16  "bbb>36=?"
  17  PROMPT
  18  INT
  19  STO 20
  20  CLX
  21  STO 00
  22  STO 01
  23  STO 02
  24  STO 03
  25  GTO 10
  26  LBL 01
  27  ENTER
  28  ENTER
  29  X^2
  30  1
  31  -
  32  STO 06
  33  *
  34  STO 04
  35  STO 08
  36  CLX
  37  SIGN
  38  ST+ Z
  39  -
  40  R^
  41  ST* Z
  42  *
  43  STO 05
  44  CLX
  45  +
  46  STO 07           
  47  CLX
  48  2
  49  ST+ Z
  50  -
  51  ST* 04
  52  X<>Y
  53  ST* 08
  54  *
  55  6
  56  *
  57  ST* 06
  58  1.5
  59  /
  60  CHS
  61  ST* 05
  62  ST* 07
  63  RTN
  64  LBL 02
  65  ENTER
  66  X^2
  67  STO T
  68  *
  69  ST+ X
  70  STO Y
  71  6
  72  *
  73  STO 11
  74  X<>Y
  75  ENTER
  76  STO 09
  77  R^
  78  3
  79  *
  80  1
  81  -
  82  ST- 09
  83  +
  84  STO 13           
  85  LASTX
  86  ST+ X
  87  6
  88  -
  89  STO 10
  90  CHS
  91  R^
  92  4
  93  *
  94  ST- 10
  95  -
  96  STO 12
  97  RTN
  98  LBL 10
  99  SF 01
100  SF 02
101  LBL 11
102  RCL 03
103  XEQ 02
104  RCL 02
105  XEQ 01
106  4.014005
107  REGMOVE
108  RCL 01
109  XEQ 01
110  25
111  FS? 02
112  SQRT
113  FS? 01
114  SIGN
115   E5
116  /
117  .9
118  +
119  RCL 20
120  +
121  STO 21
122  4.008
123  STO 22           
124  5.005
125  +
126  STO 24
127  LASTX
128  +
129  STO 23
130  26.9
131  STO 25
132  CLST
133  STO 19
134  LBL 12
135  RCL IND 21
136  RCL IND 22
137  *
138  +
139  ISG 21
140  ISG 22
141  GTO 12
142  RCL IND 23
143  *
144  +
145  5
146  ST- 22
147  CLX
148  124
149  FC? 01
150  FS? 02
151  CLX
152  ST- 21
153  CLX
154  ISG 23
155  GTO 12
156  X<>Y
157  STO IND 25
158  RCL IND 24
159  *
160  ST+ 19
161  5
162  ST- 23
163  124
164  FS? 02
165  FS? 01
166  CLX
167  ST- 21
168  CLST
169  ISG 25
170  ISG 24
171  GTO 12
172  RCL 19           
173  RCL 26
174  RCL 30
175  +
176  RCL 27
177  RCL 29
178  +
179  16
180  *
181  -
182  RCL 28
183  30
184  *
185  +
186  /
187  ENTER
188  X<> 03
189  -
190  ABS
191  ST+ 00
192  RCL 03
193  X<> 02
194  X<> 01
195  STO 03
196  FS?C 01
197  GTO 11
198  FS?C 02
199  GTO 11
200  VIEW 01
201  CLX
202  X<> 00
203  RND
204  X#0?
205  GTO 10
206  RCL 03
207  XEQ 01
208  4.009005
209  REGMOVE
210  RCL 02           
211  XEQ 01
212  4.014005
213  REGMOVE
214  RCL 01
215  XEQ 01
216  5
217  ST- 24
218  RCL 20
219  STO 21
220  CLST
221  STO 19
222  LBL 13
223  RCL IND 21
224  RCL IND 22
225  *
226  +
227  ISG 21
228  CLX
229  ISG 22
230  GTO 13
231  RCL IND 23
232  *
233  +
234  5
235  ST- 22
236  CLX
237  ISG 23
238  GTO 13
239  CLX
240  RCL IND 24
241  *
242  ST+ 19
243  5
244  ST- 23
245  CLST
246  ISG 24
247  GTO 13
248  13824
249  ST/ 19
250  RCL 03           
251  RCL 36
252  *
253  RCL 33
254  +
255  RCL 02
256  RCL 35
257  *
258  RCL 32
259  +
260  RCL 01
261  RCL 34
262  *
263  RCL 31
264  +
265  RCL 19
266  RDN
267  CLD
268  END

 
     ( 465 bytes / SIZE bbb+125 )
 
 

      STACK       INPUT1        INPUT2      INPUT3      OUTPUTS
           T            /             /           /        f(x,y,z)
           Z            h             x0           /            z
           Y            k             y0           /            y
           X            l             z0         bbb            x

     bbb > 36  ,  bbb = 1st register of the hypermatrix   ->     f(x,y,z) = f-extremum

Example:    The following routine stores in registers R01 thru R125 the values of the function

   f(x,y,z) = 4 + ( x2 - 3 x + 1 )2 + ( y2 - 4 y + 1 )2 + ( z2 - 5 z + 1 )2     for  x , y , z  =  -2 , -1 , 0 , +1 , +2
 
 

 01  LBL "FF"
 02  5
 03  STO M
 04  STO N
 05  STO O
 06  125
 07  STO 00        
 08  LBL 01
 09  RCL M
 10  3
 11  -
 12  ENTER
 13  ENTER
 14  3
 15  -
 16  *
 17  1
 18  +
 19  X^2
 20  RCL N         
 21  3
 22  -
 23  ENTER
 24  ENTER
 25  4
 26  -
 27  *
 28  1
 29  +
 30  X^2
 31  +
 32  RCL O         
 33  3
 34  -
 35  ENTER
 36  ENTER
 37  5
 38  -
 39  *
 40  1
 41  +
 42  X^2
 43  +
 44  4
 45  +
 46  STO IND 00
 47  DSE 00        
 48  CLX
 49  DSE M
 50  GTO 01
 51  5
 52  STO M
 53  DSE N
 54  GTO 01
 55  STO N         
 56  DSE O
 57  GTO 01
 58  CLA
 59  END

 
1°)  XEQ "FF"

2°)  If you choose, for instance, bbb = 41:      1.041125  REGMOVE  places thess 125 numbers into  R41 thru R165 , then

3°)  XEQ "EXYZ125"     >>>>     H^K^L=?

       1  ENTER^
           ENTER^   R/S     >>>>      X^Y^Z=?

       0  ENTER^
           ENTER^   R/S     >>>>      bbb>36=?

     41   FIX 8      R/S

-The HP-41 displays the successive h-approximations and after 11 seconds with V41 - almost 2 hours with HP-41 - you get:

                        x = 0.381966008
          RDN      y = 0.267949193
          RDN      z = 0.208712153
          RDN   fmin = 3.999999999   ( the exact value is of course 4 )

Note:

-In this example,  x , y , z  are roots of the 3 polynomials

      x2 - 3 x + 1  ,  y2 - 4 y + 1  ,  z2 - 5 z + 1
 

4°) N-dimensional Problems
 

     4-a) Program#1
 

-This program finds a minimum of a given function f.
-Seek a minimum of (-f)  if you want to know a maximum of the function.
 

Data Registers:           •  R00 = function name                                   ( R00 and R11 to R10+n  are to be initialized before executing "EXN" )

                                          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

-Line 14 = VIEW 10  if you want to display the f-values   or   VIEW 11  if you prefer the x-values ( in this case, add  CLD  after line 44 )
 
 

01  LBL "EXN"
02  STO 09
03  X<>Y
04  STO 08
05  XEQ IND 00
06  STO 10
07  GTO 01
08  LBL 00
09  X<>Y
10  STO 10        
11  DSE 07
12  GTO 02
13  LBL 01
14  VIEW 10 
15  RCL 09
16  STO 06
17  10.01
18  +
19  STO 07
20  LBL 02
21  RCL 08
22  ST+ IND 07
23  XEQ IND 00
24  RCL 10
25  X>Y?
26  GTO 00
27  RCL 08
28  ST+ X
29  ST- IND 07
30  XEQ IND 00
31  RCL 10
32  X>Y?
33  GTO 00
34  RCL 08
35  ST+ IND 07
36  DSE 06
37  GTO 00
38  PI
39  ST/ 08
40  RCL 08 
41  RND
42  X#0?
43  GTO 01
44  RCL 10        
45  END

 
     ( 74 bytes / SIZE nnn+11 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h             0
           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 f-values are displayed ( the same values may be displayed several times ) and  4mn28s 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 mentioned for "EXY" apply a fortiori to "EXN"
 

-This program may be simplified if we assume that, for instance, n does not exceed 10: we can store x1 , x2 , ...... , xn  into  R01 , R02 , ...... , Rnn
  which is more natural and 6 extra-bytes can be saved.
 

Data Registers:           •  R00 = function name                      ( R00 to Rnn  are to be initialized before executing "EXN2" )

                                      •  R01 = x1  •  R02 = x2  ...........   •  Rnn = xn             R11 = n   R12 = h   R13 = f ( x1 , x2 , ...... , xn )    R14-R15: temp
Flags: /
Subroutine:   A program which computes f ( x1 , x2 , ...... , xn ) assuming x1 is in R01 , x2 is in R02 , ............ , xn is in Rnn
 
 

01  LBL "EXN2"
02  STO 11
03  X<>Y
04  STO 12
05  XEQ IND 00
06  STO 13
07  GTO 01
08  LBL 00
09  X<>Y
10  STO 13
11  DSE 14
12  GTO 02 
13  LBL 01
14  VIEW 13 
15  RCL 11        
16  STO 14
17  STO 15
18  LBL 02
19  RCL 12
20  ST+ IND 14
21  XEQ IND 00
22  RCL 13
23  X>Y?
24  GTO 00
25  RCL 12
26  ST+ X
27  ST- IND 14
28  XEQ IND 00
29  RCL 13
30  X>Y?
31  GTO 00
32  RCL 12
33  ST+ IND 14
34  DSE 15
35  GTO 00
36  PI
37  ST/ 12
38  RCL 12 
39  RND
40  X#0?
41  GTO 01 
42  RCL 13         
43  END

 
  ( 69 bytes / SIZE 016 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h             0
           X             n          min(f)

-The same example is solved in 4mn21s
 

     4-b) Program#2 - Parabolic Interpolation
 

-The following program seeks an extremum of a function of n variables  f ( x1 , x2 , ...... , xn )   provided  n does not exceed 10
 

Data Registers:           •  R00 = function name                                   ( R00 and R01 to Rnn  are to be initialized before executing "PEXN" )

                                      •  R01 = x1  •  R02 = x2  ...........   •  Rnn = xn        R11 = n       R12 = f ( x1 , x2 , ...... , xn )          R13 to R20+3n: temp
Flags: /
Subroutine:   A program which computes f ( x1 , x2 , ...... , xn ) assuming x1 is in R01 , x2 is in R02 , ............ , xn is in Rnn

-The accuracy is controlled by the display setting.
 
 

  01  LBL "PEXN"
  02  STO 11
  03  X<>Y
  04  STO 12
  05  CLX
  06  3
  07  *
  08  20
  09  +
  10  STO 19
  11  RCL 11
  12  LBL 00
  13  RCL IND X 
  14  RCL 12
  15  -
  16  STO IND Z
  17  RDN
  18  DSE Y
  19  DSE X
  20  GTO 00
  21  X<>Y
  22  STO 18
  23  RCL 11
  24  LBL 01
  25  RCL IND X
  26  STO IND Z
  27  RDN
  28  DSE Y
  29  DSE X
  30  GTO 01
  31  X<>Y
  32  STO 17
  33  RCL 11
  34  LBL 02
  35  RCL IND X 
  36  RCL 12
  37  +
  38  STO IND Y
  39  STO IND Z
  40  RDN
  41  DSE Y
  42  DSE X
  43  GTO 02
  44  STO 15
  45  LBL 03
  46  RCL 11
  47  STO 16
  48  CLX
  49  STO 20
  50  VIEW 01
  51  LBL 04
  52  RCL IND 19
  53  STO IND 16
  54  XEQ IND 00
  55  STO 14
  56  RCL IND 18
  57  STO IND 16
  58  XEQ IND 00
  59  STO 13
  60  RCL IND 17
  61  STO IND 16
  62  XEQ IND 00
  63  STO 12
  64  RCL 13
  65  -
  66  RCL IND 17
  67  RCL IND 18
  68  -
  69  X#0?
  70  /
  71  RCL 13
  72  RCL 14
  73  -
  74  RCL IND 18
  75  RCL IND 19
  76  -
  77  X#0?
  78  /
  79  -
  80  RCL IND 19
  81  RCL IND 17
  82  -
  83  X#0?
  84  /
  85  X#0?
  86  STO 15
  87  RCL IND 18
  88  RCL IND 17
  89  -
  90  *
  91  +
  92  RCL 15
  93  ST+ X
  94  X#0?
  95  /
  96  ABS
  97  ST+ 20
  98  LASTX
  99  RCL IND 16
100  +
101  STO IND 16
102  X<> IND 17
103  X<> IND 18
104  STO IND 19
105  DSE 19
106  DSE 18
107  DSE 17
108  DSE 16
109  GTO 04
110  RCL 20
111  RCL 11
112  ST+ 17
113  ST+ 18
114  ST+ 19
115  /
116  RND
117  X#0?
118  GTO 03
119  RCL 15
120  X#0?
121  SIGN
122  RCL 12
123  CLD
124  END

 
   ( 190 bytes / SIZE 3n+21 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h          +/-1
           X             n    extremum(f)

  Y-output = +1 suggests a maximum
  Y-output = -1 suggests a minimum         ( Y = 0 suggests nothing! )

Example:   Let's calculate the distance between the paraboloid  (P)  z = x2 + 2 y2  and the paraboloid  (P')  z' = PI - 2 ( x' - 7 )2 - ( y' - 5 )2

-We load the following routine that computes the square of the distance between one point M(x,y) on the first paraboloid
  and another point M'(x',y') on the second one.

    LBL "T"      ST+ X        7              RCL 04      X^2              +                  +                    assuming      R01 = x        R03 = y
    RCL 01      +                 -              5                RCL 01         RCL 03       RTN                  and           R02 = x'       R04 = y'
    X^2            PI               X^2         -                 RCL 02         RCL 04
    RCL 03      -                 ST+ X     X^2            -                    -
    X^2            RCL 02      +             +                X^2               X^2

    "T"   ASTO 00

    1   STO 01    STO 03     if we take   x =  y = 1
    6   STO 02    STO 04          and       x' = y' = 6    as initial guesses

   n = 4  ( four unknowns )   and with  h = 1  and  FIX 4

   1    ENTER^
   4    XEQ "PEXN"   the successive x-values are displayed and eventually, it yields:

   extremum(f)  =  37.78538053                         ( execution time = 6 minutes )
        RDN          -1  ( a minimum )

>>>   whence the distance between these paraboloids = SQRT(37.78538053) = 6.146981416    ( the last digit should be a "2" )

-We have also:

       x = 1.4552     x' = 6.2724
       y = 0.5197     y' = 3.9606