HP-71B Program to Find a Function Minimum

Using Newton's Method

by Namir Shammas

The following program calculates the minimum point of a multi-variable function using Newton's method. This method is implemented in a basic form without any enhancements.

The program prompts you to enter for each variable:

1. The maximum number of iterations

2. The values for the initial set of variables

3. The values for the tolerances for each variable.

The program displays intermediate results to show the iteration progress and also the following final results:

1. The coordinates of the minimum value.

2. The fine step sizes for each variable.

3. The minimum function value.

4. The number of iterations

The current code finds the minimum for the following function:

f(x1,x2) = x1 - x2 + 2 * x1 ^ 2 + 2 * x1 * x2 + x2 ^ 2

Here is a sample session. The first figure shows the input.:

The second figure shows the output:

Here is the BASIC listing:  

10 ! Basic Newton's Method for Optimization
20 !
30 ! USES JPCROM TO SUPPORT ADVANCED LOOPS AND
40 ! DECISION-MAKING CONSTRUCTS
50 !
60 ! USES MATH ROM TO SUPPORT MAT STATEMENTS
70 !
80 N = 2
90 !
100 DIM X(N), G(N)
110 DIM T0(N)
120 DIM D(N), J(N, N)
130 !
140 DISP "NEWTON'S OPTIMIZATIOM METHOD"
150 INPUT "Enter maximum number iterations? "; I9
160 For I = 1 To N
170 DISP "Enter guess for X(";I;")";
180 INPUT X(I)
190 DISP "Enter tolerance for X(";I;")";
200 INPUT T0(I)
210 Next I
220 !
230 I0 = 0
240 REPEAT
250 I0 = I0 + 1
260 If I0 > I9 Then
270 DISP "Reached maximum iterations limit"
280 LEAVE
290 End If
300 REM
310 CALL GetDrv1(N, X, G)
320 REM
330 ! test if gradient is shallow enough
340 N0 = 0
350 For I = 1 To N
360 N0 = N0 + G(I)^2
370 Next I
380 N0 = Sqr(N0)
390 ! If N0 < E Then LEAVE
400 REM
410 CALL GetDrv2(N, X, J)
420 MAT J = INV(J)
430 MAT G = J * G
440 For I = 1 To N
450 D(I) = G(I)
460 X(I) = X(I) - D(I)
470 Next I
480 REM
490 S0 = 1
500 For I = 1 To N
510 If Abs(D(I)) > T0(I) Then
520 S0 = 0
530 I = N
540 End If
550 Next I
560 REM
570 CALL MyFx(X, N, F)
580 DISP "F = ";F;" ";
590 For I = 1 To N
600 DISP "X=(";I;")=";X(I);" Delta(";I;")=";D(I);" ";
610 Next I
620 DISP
630 REM
640 UNTIL S0 <> 1
650 REM
660 CALL MyFx(X, N, F)
670 DISP "**********FINAL RESULTS************"
680 DISP "Optimum at:"
690 For I = 1 To N
700 DISP "X(";I;")=";X(I) @ PAUSE
710 Next I
720 For I = 1 To N
730 DISP "Delta X(";I;")=";D(I) @ PAUSE
740 Next I
750 DISP "Function value ="; F @ PAUSE
760 DISP "Number of iterations = ";I0
770 !
780 !
790 SUB MyFx(X(), N, R)
800 ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
810 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
820 END SUB
830 !
840 SUB FrstDrv(N, X(), I0, R)
850 X0 = X(I0)
860 h = 0.01 * (1 + Abs(X0))
870 X(I0) = X0 + h
880 CALL MyFx(X, N, F1)
890 X(I0) = X0 - h
900 CALL MyFx(X, N, F2)
910 X(I0) = X0
920 R = (F1 - F2) / 2 / h
930 END SUB
940 !
950 SUB SecDerv(N, X(), I0, J0, R)
960 ! calculate second derivative?
970 If I0 = J0 Then
980 CALL MyFx(X, N, F0)
990 X0 = X(I0)
1000 H1 = 0.01 * (1 + Abs(X0))
1010 X(I0) = X0 + H1
1020 CALL MyFx(X, N, F1)
1030 X(I0) = X0 - H1
1040 CALL MyFx(X, N, F2)
1050 X(I0) = X0
1060 R = (F1 - 2 * F0 + F2) / H1 ^ 2
1070 Else
1080 X0 = X(I0)
1090 Y0 = X(J0)
1100 H1 = 0.01 * (1 + Abs(X0))
1110 H2 = 0.01 * (1 + Abs(Y0))
1120 ! calculate F3
1130 X(I0) = X0 + H1
1140 X(J0) = Y0 + H2
1150 CALL MyFx(X, N, F3)
1160 ! calculate F4
1170 X(I0) = X0 - H1
1180 X(J0) = Y0 - H2
1190 CALL MyFx(X, N, F4)
1200 ! calculate F5
1210 X(I0) = X0 + H1
1220 X(J0) = Y0 - H2
1230 CALL MyFx(X, N, F5)
1240 ! calculate F6
1250 X(I0) = X0 - H1
1260 X(J0) = Y0 + H2
1270 CALL MyFx(X, N, F6)
1280 X(I0) = X0
1290 X(J0) = Y0
1300 R = (F3 - F6 - F5 + F4) / (4 * H1 * H2)
1310 End If
1320 END SUB
1330 !
1340 SUB GetDrv1(N, X(), D1())
1350 For I = 1 To N
1360 CALL FrstDrv(N, X, I, D1(I))
1370 Next I
1380 END SUB
1390 !
1400 SUB GetDrv2(N, X(), D2(,))
1410 For I = 1 To N
1420 For J = 1 To N
1430 CALL SecDerv(N, X, I, J, D2(I, J))
1440 Next J
1450 Next I
1460 END SUB
The above listing uses the JPCROM to employ more structured constructs for loops and decision-making. The above listing uses the REPEAT-UNTIL and LOOP-ENDLOOP loops and also the IF-THEN-ELSE-ENDIF enhanced version of the basic IF statement. These program flow control statements limits or eliminates using the GOTO statements and make the programs easier to read and update. Here is the version of the listing that is indented and void of line numbers. I am including it here since it is easier to read than the above listing although it will not run on either the HP-71B handheld computer or the EMU71 emulator..   
! Basic Newton's Method for Optimization
!
! USES JPCROM TO SUPPORT ADVANCED LOOPS AND 
! DECISION-MAKING CONSTRUCTS
!
! USES MATH ROM TO SUPPOT MAT STATEMENTS
!
N = 2
!
DIM X(N), G(N) 
DIM T0(N) 
DIM D(N), J(N, N) 
!
DISP "NEWTON'S OPTIMIZATIOM METHOD"
INPUT "Enter maximum number iterations? "; I9
For I = 1 To N
  DISP "Enter guess for X(";I;")";
  INPUT X(I)
  DISP "Enter tolerance for X(";I;")";
  INPUT T0(I)  
Next I
!    
I0 = 0
REPEAT
  I0 = I0 + 1
  If I0 > I9 Then
    DISP "Reached maximum iterations limit"
    LEAVE
  End If
  
  CALL GetDrv1(N, X, G)
  
  ! test if gradient is shallow enough
  N0 = 0
  For I = 1 To N
    N0 = N0 + G(I)^2
  Next I  
  N0 = Sqr(N0)
  ! If N0 < E Then LEAVE
  
  CALL GetDrv2(N, X, J)
  MAT J = INV(J)
  MAT G = J * G
  For I = 1 To N
    D(I) = G(I)
    X(I) = X(I) - D(I)
  Next I
  
  S0 = 1
  For I = 1 To N
    If Abs(D(I)) > T0(I) Then
      S0 = 0
      I = N
    End If
  Next I
  
  CALL MyFx(X, N, F)
  DISP "F = ";F;" ";
  For I = 1 To N
    DISP "X=(";I;")=";X(I);" Delta(";I;")=";D(I);" ";
  Next I
  DISP

UNTIL S0 <> 1

CALL MyFx(X, N, F)
DISP "**********FINAL RESULTS************"
DISP "Optimum at:"
For I = 1 To N
  DISP "X(";I;")=";X(I) @ PAUSE
Next I
For I = 1 To N
  DISP "Delta X(";I;")=";D(I) @ PAUSE
Next I
DISP "Function value ="; F @ PAUSE
DISP "Number of iterations = ";I0
! 
!  
SUB MyFx(X(), N, R)
  ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
  R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
END SUB
!
SUB FrstDrv(N, X(), I0, R) 
  X0 = X(I0)
  h = 0.01 * (1 + Abs(X0))
  X(I0) = X0 + h
  CALL MyFx(X, N, F1)
  X(I0) = X0 - h
  CALL MyFx(X, N, F2)
  X(I0) = X0
  R = (F1 - F2) / 2 / h
END SUB
!
SUB SecDerv(N, X(), I0, J0, R) 
  ! calculate second derivative?
  If I0 = J0 Then
    CALL MyFx(X, N, F0)
    X0 = X(I0)
    H1 = 0.01 * (1 + Abs(X0))
    X(I0) = X0 + H1
    CALL MyFx(X, N, F1)
    X(I0) = X0 - H1
    CALL MyFx(X, N, F2)
    X(I0) = X0
    R = (F1 - 2 * F0 + F2) / H1 ^ 2
  Else
    X0 = X(I0)
    Y0 = X(J0)
    H1 = 0.01 * (1 + Abs(X0))
    H2 = 0.01 * (1 + Abs(Y0))
    ! calculate F3
    X(I0) = X0 + H1
    X(J0) = Y0 + H2
    CALL MyFx(X, N, F3)
    ! calculate F4
    X(I0) = X0 - H1
    X(J0) = Y0 - H2
    CALL MyFx(X, N, F4)
    ! calculate F5
    X(I0) = X0 + H1
    X(J0) = Y0 - H2
    CALL MyFx(X, N, F5)
    ! calculate F6
    X(I0) = X0 - H1
    X(J0) = Y0 + H2
    CALL MyFx(X, N, F6)
    X(I0) = X0
    X(J0) = Y0
    R = (F3 - F6 - F5 + F4) / (4 * H1 * H2)
  End If
END SUB
!
SUB GetDrv1(N, X(), D1())  
  For I = 1 To N
    CALL FrstDrv(N, X, I, D1(I))
  Next I
END SUB
!
SUB GetDrv2(N, X(), D2(,))
  For I = 1 To N 
    For J = 1 To N 
      CALL SecDerv(N, X, I, J, D2(I, J))
    Next J
  Next I
END SUB

Here is a version of the first listing that DOES NOT rely on the JPCROM. This version should run on the basic HP-71B with the Math module.

10 ! Basic Newton's Method for Optimization
20 !
60 ! USES MATH ROM TO SUPPOT MAT STATEMENTS
70 !
80 N = 2
90 !
100 DIM X(N), G(N)
110 DIM T0(N)
120 DIM D(N), J(N, N)
130 !
140 DISP "NEWTON'S OPTIMIZATIOM METHOD"
150 INPUT "Enter maximum number iterations? "; I9
160 For I = 1 To N
170 DISP "Enter guess for X(";I;")";
180 INPUT X(I)
190 DISP "Enter tolerance for X(";I;")";
200 INPUT T0(I)
210 Next I
220 !
230 I0 = 0
240 !REPEAT
250 I0 = I0 + 1
260 If I0 <= I9 Then 290
270 DISP "Reached maximum iterations limit"
280 LEAVE
290 !End If
300 REM
310 CALL GetDrv1(N, X, G)
320 REM
330 ! test if gradient is shallow enough
340 N0 = 0
350 For I = 1 To N
360 N0 = N0 + G(I)^2
370 Next I
380 N0 = Sqr(N0)
390 ! If N0 < E Then LEAVE
400 REM
410 CALL GetDrv2(N, X, J)
420 MAT J = INV(J)
430 MAT G = J * G
440 For I = 1 To N
450 D(I) = G(I)
460 X(I) = X(I) - D(I)
470 Next I
480 REM
490 S0 = 1
500 For I = 1 To N
510 If Abs(D(I)) <= T0(I) Then 540
520 S0 = 0
530 I = N
540 !End If
550 Next I
560 REM
570 CALL MyFx(X, N, F)
580 DISP "F = ";F;" ";
590 For I = 1 To N
600 DISP "X=(";I;")=";X(I);" Delta(";I;")=";D(I);" ";
610 Next I
620 DISP
630 REM
640 IF S0 = 1 THEN 240
650 REM
660 CALL MyFx(X, N, F)
670 DISP "**********FINAL RESULTS************"
680 DISP "Optimum at:"
690 For I = 1 To N
700 DISP "X(";I;")=";X(I) @ PAUSE
710 Next I
720 For I = 1 To N
730 DISP "Delta X(";I;")=";D(I) @ PAUSE
740 Next I
750 DISP "Function value ="; F @ PAUSE
760 DISP "Number of iterations = ";I0
770 !
780 !
790 SUB MyFx(X(), N, R)
800 ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
810 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
820 END SUB
830 !
840 SUB FrstDrv(N, X(), I0, R)
850 X0 = X(I0)
860 h = 0.01 * (1 + Abs(X0))
870 X(I0) = X0 + h
880 CALL MyFx(X, N, F1)
890 X(I0) = X0 - h
900 CALL MyFx(X, N, F2)
910 X(I0) = X0
920 R = (F1 - F2) / 2 / h
930 END SUB
940 !
950 SUB SecDerv(N, X(), I0, J0, R)
960 ! calculate second derivative?
970 If I0 <> J0 Then 1070
980 CALL MyFx(X, N, F0)
990 X0 = X(I0)
1000 H1 = 0.01 * (1 + Abs(X0))
1010 X(I0) = X0 + H1
1020 CALL MyFx(X, N, F1)
1030 X(I0) = X0 - H1
1040 CALL MyFx(X, N, F2)
1050 X(I0) = X0
1060 R = (F1 - 2 * F0 + F2) / H1 ^ 2
1065 GOTO 1310
1070 !Else
1080 X0 = X(I0)
1090 Y0 = X(J0)
1100 H1 = 0.01 * (1 + Abs(X0))
1110 H2 = 0.01 * (1 + Abs(Y0))
1120 ! calculate F3
1130 X(I0) = X0 + H1
1140 X(J0) = Y0 + H2
1150 CALL MyFx(X, N, F3)
1160 ! calculate F4
1170 X(I0) = X0 - H1
1180 X(J0) = Y0 - H2
1190 CALL MyFx(X, N, F4)
1200 ! calculate F5
1210 X(I0) = X0 + H1
1220 X(J0) = Y0 - H2
1230 CALL MyFx(X, N, F5)
1240 ! calculate F6
1250 X(I0) = X0 - H1
1260 X(J0) = Y0 + H2
1270 CALL MyFx(X, N, F6)
1280 X(I0) = X0
1290 X(J0) = Y0
1300 R = (F3 - F6 - F5 + F4) / (4 * H1 * H2)
1310 !End If
1320 END SUB
1330 !
1340 SUB GetDrv1(N, X(), D1())
1350 For I = 1 To N
1360 CALL FrstDrv(N, X, I, D1(I))
1370 Next I
1380 END SUB
1390 !
1400 SUB GetDrv2(N, X(), D2(,))
1410 For I = 1 To N
1420 For J = 1 To N
1430 CALL SecDerv(N, X, I, J, D2(I, J))
1440 Next J
1450 Next I
1460 END SUB

BACK

Copyright (c) Namir Shammas. All rights reserved.