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: 4 tabulated values
1-e) 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 |

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
RDN    fmin =  2.9583

1-c) Discrete Data:  4 Tabulated Values

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

Data Registers:               R00 = h                                ( Registers R01 thru R04 are to be initialized before executing "DEX4" )

•  R01 = f(x-1)    •  R02 = f(x0)    •  R03 = f(x1)    •  R04 = f(x2)                  R05 to R08: temp
Flags: /
Subroutines: /

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

( 100 bytes / SIZE 009 )

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

Example1:     Estimate the maximum of f(x) from these 4 values

 x 1 3 5 7 y 2 6 7 4

2  STO 01   6  STO 02   7  STO 03   4  STO 04

h = 2  and  x0 = 3

2  ENTER^
3  XEQ "DEX4"   >>>>     x   =  4.571877795
RDN    fmax =  7.088374732

Example2:     Estimate the maximum of f(x) from these 4 values

 x 1 3 5 7 y 2 6 7 5

2  STO 01   6  STO 02   7  STO 03   5  STO 04

h = 2  and  x0 = 3

2  ENTER^
3  XEQ "DEX4"   >>>>     x   =  4.666666667
RDN    fmax =  7.041666667

Notes:

-In example 2, the 4 points are on a parabola.
-Otherwise, if f is not a plolynomial of degree 2, lines 40-41 & 56 to 60 are superfluous.
-So, in most cases, they could be deleted.

-In fact, if you delete these lines and if you replace 5 by 4.999999999  STO 04 in the 2nd example,
you will get x = 4.666666666 & fmax = 7.041666667  almost exact results !

-This routine solves a quadratic equation to find the extremum of a polynomial of degree 3.
-Lines 51 to 55 choose the solution that is the closest to x0

1-e) 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  XY 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  XY 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

 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