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	x₀	x

Example: Find the minimum of the function f(x) = | x² - 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 x₀ = 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 f_min = 10^-9 = R02

-Here, the solution was trivial: x = SQRT(2) , f_min = 0

     y
      |     *
      |      *              *                  *
      |         *     *          *           *
      |             *                *       *
      |                                 *   *
      |                                   *
   --|--------x₁-------------x₂---------------------------- 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 = x₁
-Choose another guess to find x₂

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 x² + b x + c passing through 3 points ( x₁ , f(x₁) ) , ( x₂ , f(x₂) ) , ( x₃ , f(x₃) )
-It calculates the abscissa of the extremum x₄ = -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 = x₃       R04 = f(x₃)        R07 = -a
                                          R02 = x₂       R05 = f(x₂)
                                          R03 = x₁       R06 = f(x₁)
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	x₁	+/-1
Y	x₂	f(x)
X	x₃	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/x⁵ )/[ 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-pointinterpolation formula.
-We assume that the x_i's are equally spaced. ( x_i+1 - x_i = h , i = -1 , 0 )
-So, f(x₀+p.h) - f(x₀) = A.p + B.p²

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

• R01 = f(x_-1) • R02 = f(x₀) • R03 = f(x₁) 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	x₀	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 x₀ = 5

2 ENTER^
5 XEQ "DEX3" >>>> x = 5.3333
RDN f_min = 2.9583

1-c) Discrete Data: 4 Tabulated Values

-DEX4 uses the Lagrange 4-pointinterpolation formula.
-We assume that the x_i's are equally spaced. ( x_i+1 - x_i = h , i = -1 , 0 , 1 )
-So, f(x₀+p.h) - f(x₀) = A.p + B.p² + C.p³

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

• R01 = f(x_-1) • R02 = f(x₀) • R03 = f(x₁) • R04 = f(x₂) 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	x₀	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 x₀ = 3

2 ENTER^
3 XEQ "DEX4" >>>> x = 4.571877795
RDN f_max = 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 x₀ = 3

2 ENTER^
3 XEQ "DEX4" >>>> x = 4.666666667
RDN f_max = 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 & f_max = 7.041666667 almost exact results !

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

1-e) Discrete Data: 5 Tabulated Values

-The following program uses the Lagrange 5-pointinterpolation formula.
-We assume that the x_i's are equally spaced. ( x_i+1 - x_i = h , i = -2 , -1 , 0 , 1 )
-So, f(x₀+p.h) - f(x₀) = A.p + B.p² + C.p³ + D.p⁴

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

• R01 = f(x_-2) • R02 = f(x_-1) • R03 = f(x₀) • R04 = f(x₁) • R05 = f(x₂) 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	x₀	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 x₀ = 5

2 ENTER^
5 XEQ "DEX5" >>>> x = 5.316624790 ( in 6 seconds )
RDN f_min = 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	x₀	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 x₀ = 5

      FIX 9
2 ENTER^
5 XEQ "DEX5" >>>>   x   = 5.316624790       ( in 6 seconds )
                             RDN f_min = 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	y₀	y
X	x₀	x

Example: Find the minimum of the function f(x,y) = | x.y - PI | + | x² - 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 x₀ = y₀ = 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    f_min = 10^-9           = R03

-The solution was again trivial: x = SQRT(2) , y = PI/SQRT(2) , f_min = 0

Warning: To find f(x,y) minimum, we usually solve the system: f_x' = 0 ; f_y' = 0 ( equating the partial derivatives to zero )

-But the solution of this system is not always aminimum of f !
-Likewise, "EXY" may converge to any values ( x₀ , y₀ )
such that r(x) = f ( x , y₀ ) has a minimum for x = x₀
and s(y) = f ( x₀ , y ) has a minimum for y = y₀
-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 = x₃        R03 = x₂          R05 = x₁        R07 = f(x,y)            R10 = -a_x
                                          R02 = y₃        R04 = y₂          R06 = y₁       R08-R09: temp        R11 = -a_y
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	y₀	y
X	x₀	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) )² - ( y - sqrt(3) )² ]

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 x₀ = y₀ = 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 x_i's are equally spaced and that the y_j's are equally spaced. ( x_i+1 - x_i = h & y_j+1 - y_j = k i , j = -1 , 0 )

-So, f(x₀+p.h,y₀+q.k) - f(x₀,y₀) = A.p + B.q + C.p² + D.p.q + E.q² + F.p²q + G.p.q² + H.p²q²

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

       2.A = f₁₀ - f_-10                      4.D = f_-1-1 - f_1-1 - f_-11 + f₁₁           F = D - B - ( f_-1-1 - f_-11 )/2
       2.B = f₀₁ - f_0-1                      2.E =   f_0-1 - 2 f₀₀ + f₀₁                  G = D - A - ( f_-1-1 + f_-11 )/2                             with   f_{i j} = f(x₀+i.h,y₀+j.k)
       2.C = f_-10 - 2 f₀₀ + f₁₀                                                                H = - C - E - f₀₀ + D + ( f_-11 + f_1-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.q² )/( 2.C + 2.F.q + 2.H.q² ) and this system is solved
q = ( - B - D.p - F.p²)/( 2.E + 2.G.p + 2.H.p² ) 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,y₀)    • R07 = f(x_-1,y₁)
                                      • R02 = f(x₀,y_-1)     • R05 = f(x₀,y₀)     • R08 = f(x₀,y₁)           R00 & R11 to R23: temp   ( R11 = p , R12 = q )
                                      • R03 = f(x₁,y_-1)     • R06 = f(x₁,y₀)     • R09 = f(x₁,y₁)
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	y₀	y
X	x₀	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 x₀ = 3 h = 3-1 = 5-3 = 2
y₀ = 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 x_i's are equally spaced and that the y_j's are equally spaced. ( x_i+1 - x_i = h & y_j+1 - y_j = k i , j = -2 , -1 , 0 , 1 )

-We have then

576 f(x₀+p.h,y₀+q.k)   =   (q²-1)q(q-2) [ (p²-1)p(p-2) f_-2-2 - 4.(p-1)p(p²-4) f_-1-2 + 6.(p²-1)(p²-4) f_0-2 - 4.(p+1)p(p²-4) f_1-2 + (p²-1)p(p+2) f_-2-2 ]
                                     - 4.(q-1)q(q²-4) [ (p²-1)p(p-2) f_-2-1 - 4.(p-1)p(p²-4) f_-1-1 + 6.(p²-1)(p²-4) f_0-1 - 4.(p+1)p(p²-4) f_1-1 + (p²-1)p(p+2) f_-2-1 ]
                                     + 6.(q²-1)(q²-4) [ (p²-1)p(p-2) f_-20 - 4.(p-1)p(p²-4) f_-10 + 6.(p²-1)(p²-4) f₀₀ - 4.(p+1)p(p²-4) f₁₀ + (p²-1)p(p+2) f_-20 ]
                                    - 4.(q+1)q(q²-4) [ (p²-1)p(p-2) f_-21 - 4.(p-1)p(p²-4) f_-11 + 6.(p²-1)(p²-4) f₀₁ - 4.(p+1)p(p²-4) f₁₁ + (p²-1)p(p+2) f_-21 ]
                                      + (q²-1)q(q+2) [ (p²-1)p(p-2) f_-22 - 4.(p-1)p(p²-4) f_-12 + 6.(p²-1)(p²-4) f₀₂ - 4.(p+1)p(p²-4) f₁₂ + (p²-1)p(p+2) f_-22 ]

with f_{i j} = f(x₀+i.h,y₀+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,y₀)    • R16 = f(x_-2,y₁)    • R21 = f(x_-2,y₂)
                                      • R02 = f(x_-1,y_-2)    • R07 = f(x_-1,y_-1)    • R12 = f(x_-1,y₀)    • R17 = f(x_-1,y₁)    • R22 = f(x_-1,y₂)
                                      • R03 = f(x₀,y_-2)     • R08 = f(x₀,y_-1)     • R13 = f(x₀,y₀)     • R18 = f(x₀,y₁)     • R23 = f(x₀,y₂)
                                      • R04 = f(x₁,y_-2)     • R09 = f(x₁,y_-1)     • R14 = f(x₁,y₀)     • R19 = f(x₁,y₁)     • R24 = f(x₁,y₂)
                                      • R05 = f(x₂,y_-2)     • R10 = f(x₂,y_-1)     • R15 = f(x₂,y₀)     • R20 = f(x₂,y₁)     • R25 = f(x₂,y₂)

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	y₀	y
X	x₀	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 x₀ = 3 h = 1-(-1) = 3-1 = 5-3 = 7-5 = 2
and y₀ = 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" )

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	y₀	y
X	x₀	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

We have x₀ = 3 h = 1-(-1) = 3-1 = 5-3 = 7-5 = 2
and y₀ = 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	z₀	z
Y	y₀	y
X	x₀	x

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

-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 x₀ = y₀ = z₀ = 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 f_min = 1.4344                      ---Execution time = 4m32s---

-In FIX 5 , you'll get

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

-And with FIX 9:

                  x = 1.105422397
                  y = -0.916837262
                  z =   2.262990778   and f_min = 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 x_i's , y_j's & z_m's are equally spaced. ( x_i+1 - x_i = h , y_j+1 - y_j = k , y_m+1 - y_m = l ) ( i , j , m = -2 , -1 , 0 , 1 )

f( x₀+p.h , y₀+q.k , z₀+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,y₀,z_-2)    • Rb+15 = f(x_-2,y₁,z_-2)    • Rb+20 = f(x_-2,y₂,z_-2)
                                      • Rb+1 = f(x_-1,y_-2,z_-2)    ................................................................................................................    • Rb+21 = f(x_-1,y₂,z_-2)
                                      • Rb+2 = f(x₀,y_-2,z_-2)      ...............................................................................................................    • Rb+22 = f(x₀,y₂,z_-2)
                                      • Rb+3 = f(x₁,y_-2,z_-2)     ................................................................................................................    • Rb+23 = f(x₁,y₂,z_-2)
                                      • Rb+4 = f(x₂,y_-2,z_-2)     .................................................................................................................   • Rb+24 = f(x₂,y₂,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(x₂,y₂,z₂)

with f_{p q r} = f( x₀+p.h , y₀+q.k , z₀+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	x₀	/	z
Y	k	y₀	/	y
X	l	z₀	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 + ( x² - 3 x + 1 )² + ( y² - 4 y + 1 )² + ( z² - 5 z + 1 )² 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   f_min = 3.999999999   ( the exact value is of course 4 )

Note:

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

x² - 3 x + 1 , y² - 4 y + 1 , z² - 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 ( x₁ , x₂ , ...... , x_n )                               • R11 = x₁ • R12 = x₂ ...........   • R10+n = x_n
Flags: /
Subroutine:   A program which computes f ( x₁ , x₂ , ...... , x_n ) assuming x₁ is in R11 , x₂ is in R12 , ............ , x_n 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: x² + y² + z² = 4 ; x² - y² + z = 2 ; x + y + z² = 1 ; x - y - z = 0
by the least-squares method ( the system itself has no solution )

-Let f(x,y,z) = ( x² + y² + z² - 4 )² + ( x² - y² + z - 2 )² + ( x + y + z² - 1 )² + ( x - y - z )²

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 x₁ , x₂ , ...... , x_n 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 = x₁ • R02 = x₂ ........... • Rnn = x_n R11 = n R12 = h R13 = f ( x₁ , x₂ , ...... , x_n ) R14-R15: temp
Flags: /
Subroutine: A program which computes f ( x₁ , x₂ , ...... , x_n ) assuming x₁ is in R01 , x₂ is in R02 , ............ , x_n 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 ( x₁ , x₂ , ...... , x_n ) provided n does not exceed 10

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

• R01 = x₁ • R02 = x₂ ........... • Rnn = x_n R11 = n R12 = f ( x₁ , x₂ , ...... , x_n ) R13 to R20+3n: temp
Flags: /
Subroutine: A program which computes f ( x₁ , x₂ , ...... , x_n ) assuming x₁ is in R01 , x₂ is in R02 , ............ , x_n 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 = x² + 2 y² and the paraboloid (P') z' = PI - 2 ( x' - 7 )² - ( y' - 5 )²

-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

hp41programs

Extrema for the HP-41