hp41programs

Linear and non-Linear Systems for the HP-41
 

Overview

 1°) Linear Systems

    a) Program#1
    b) Program#2
    c1) Program#3
    c2) free42 Program     
    d) Underdetermined Systems

      d-1) Program#1
      d-2) Program#2

    e) Overdetermined Systems

      e-1) Program#1
      e-2) Program#2

    f) Tridiagonal Systems
    g) Pentadiagonal Systems

 2°) Non-Linear Systems

    a) 1 Equation in 1 unknown

      a-1)  Secant Method
      a-2)  Quadratic Interpolation
      a-3)  Ridders' Method
      a-4)  Discrete Data ( 5 points )

    b) 2 Equations in 2 unknowns
    c) 3 Equations in 3 unknowns
    d) N Equations in N unknowns

      d-1) Program#1
      d-2) Program#2

 3°) A short In/Out Routine
 

** To solve non-linear systems of equations ( including with complex unknowns ).
      see also "A Successive Approximation Method for the HP-41"

Latest Update:  paragraph1°)c2)

 

1°) Linear Systems
 

     a) Program#1
 

 "LS" allows you to solve linear systems, including overdetermined and underdetermined systems.
  You can also invert up to a 12x12 matrix.
  It is essentially equivalent to the RREF function of the HP-48 ( a "little" slower, I grant you that! ):
  Its objective is to reduce the matrix on the upper left to a diagonal form with only ones in the diagonal.
  The determinant of this matrix is also computed and stored in register R00.
  ( if there are more rows than columns, R00 is not always equal to the determinant of the upper left matrix because of rows exchanges )

  If flag F01 is clear, Gaussian elimination with partial pivoting is used.
  If flag F01 is set, the pivots are the successive elements of the diagonal:
  It's sometimes useful for matrices like Pascal's matrices of high order.
  They are extremely troublesome and many roundoff errors can occur.
  But if you set flag F01, all the coefficients will be computed with no roundoff error at all, because all the pivots will be ones!

  One advantage of this program is that you can choose the beginning data register - except R00 - ( this feature is used to solve non-linear systems too ):
  You store the first coefficient into Rbb , ... , up to the last one into Ree ( column by column )             ( bb > 00 )

  Then you key in  bbb.eeerr  ENTER^
                                 p          XQ "LS"    and the system will be solved.    ( bbb.eeerr ends up in L-register )

  where r is the number of rows of the matrix
  and  p  is a small non-negative number that you choose to determine tiny elements: if a pivot is smaller than p it will be considered to be zero.

  Don't interrupt "LS" because registers P and Q are used ( there is no risk of crash, but their contents could be altered )
  If you're not familiar to synthetic programming, you can replace status registers M N O P Q by R01 R02 R03 R04 R05
  In that case, you must choose b > 5.

-Line 121 is a three-byte  GTO 01
 
 

01  LBL "LS"    02  STO M   03  X<>Y   04  ENTER^   05  INT   06  LASTX   07  FRC   08  ISG X   09  INT   10  STO O   11  RCLY   12  +   13  DSE X   14   E3   15  /   16  +   17  X<> O          18  .1   19  %   20  +   21  STO Q   22  SIGN   23  STO 00   24  X<>Y   25  STO P    26  LBL 01

27  STO N   28  FS? 01   29  GTO 04   30  INT   31  RCL O   32  FRC   33  +   34  ENTER^   35  ENTER^   36  CLX   37  LBL 02         38  RCL IND Z   39  ABS   40  X>Y?   41  STO Z   42  X>Y?   43  +   44  RDN   45  ISG Z   46  GTO 02   47  RCL N   48  ENTER^   49  FRC   50  R^   51  INT   52  +

53  X=Y?   54  GTO 04   55  LBL 03   56  RCL IND X   57  X<> IND Z   58  STO IND Y   59  ISG Y   60  RDN   61  ISG Y   62  GTO 03   63  RCL 00   64  CHS   65  STO 00   66  LBL 04   67  RCL M   68  RCL IND N   69  ST* 00   70  ABS   71  X<=Y?   72  GTO 09   73  RCL N   74  LASTX   75  LBL 05   76  ST/ IND Y   77  ISG Y   78  GTO 05

79  LBL 06   80  RCL N   81  ENTER^   82  FRC   83  RCL O   84  INT   85  +   86  X=Y?   87  GTO 08   88  RCL IND X   89  SIGN   90  RDN   91  LBL 07   92  RCL IND Y   93  LASTX   94  *   95  ST- IND Y   96  ISG Y   97  RDN   98  ISG Y   99  GTO 07 100  LBL 08 101  ISG O 102  GTO 06 103  RCL Q 104  INT

105  ST- O 106  LBL 09        107  RCL Q 108  ST+ O 109  X^2 110  RCL P 111  INT 112  ENTER^ 113  SIGN 114  ST+ N 115  - 116  + 117  RCL N  118  X>Y? 119  CLX 120  ISG X 121  GTO 01 122  X<> P 123  SIGN 124  RCL 00 125  CLA 126  END

 
   ( 189 bytes / SIZE eee+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	bbb.eeerr	1
X	p	determinant
L	/	bbb.eeerr

 
Example1:     Find the solution of the system:

    2x + 3y + 5z + 4t = 39
   -4x + 2y + z + 3t  = 15
    3x  -  y + 2z + 3t = 19
    5x + 7y - 3z + 2t = 18

 If you choose to store the first element into R01, you have to store these 20 numbers:

          2  3  5  4  39                  R01  R05  R09  R13  R17
        -4  2  1  3  15        in        R02  R06  R10  R14  R18      respectively           ( you can use "IO" at the bottom of this page )
         3 -1  2  3  19                  R03  R07  R11  R15  R19
         5  7 -3  2  18                  R04  R08  R12  R16  R20

The first register is R01, the last register is R20 and there are 4 rows, therefore the control number of the matrix is  1.02004
If you choose p = 10^-7  key in:

                CF 01
  1.02004 ENTER^
      E-7    XEQ "LS"   and  you obtain  840 in X-register and in R00:  it is the determinant of the 4x4 matrix on the left.              ( execution time = 31s )

Registers R01 thru R16 now contains the unit matrix and registers R17 thru R20 contains the solution  x = 1 ; y = 2 ; z = 3 ; t = 4
Thus, the array has been changed into:

         1  0  0  0  1
         0  1  0  0  2              the solution is the last column
         0  0  1  0  3              of the new matrix.
         0  0  0  1  4

Example2:    Solve the system:

          5x + y + z  = 8
          4x - y + 2z = 13
           x + 2y - z  = -5
         7x - 4y + 5z = 31

-Suppose we need to store the first element into  R11 , then we store these 16 numbers

          5   1  1   8                       R11  R15  R19  R23
          4  -1  2  13         in          R12  R16  R20  R24          respectively
          1   2  -1  -5                     R13  R17  R21  R25
          7  -4  5  31                      R14  R18  R22  R26

Here, the control number of this array is   11.02604   and if we take p = 10^-7 again,

              11.02604  ENTER^
                    E-7     XEQ "LS"   gives   -5.4  10^-17  in X-register and in R00 = the determinant of the 4x4 matrix, which is of course 0   ( execution time = 21s )

and registers R11 thru R27 are now:

        1  0   0.333333333     2.333333333
        0  1  -0.666666667   -3.666666669
        0  0      5.10^-10              -2.10^-9
        0  0          0                    4.10^-9

Thus, the system is equivalent to:

      x + z /3  =  7/3              and there is an infinite set of solutions:
      y - 2z /3 = -11/3           z may be chosen at random and x and y are determined by
                0  =  0                 x = -z /3 + 7/3
                0  =  0                 y = 2z /3 - 11/3
 

Note:    If the last number of the initial matrix is 32 instead of 31 the array is changed into:

       1  0   0.333333333   0        i-e         x + z /3 = 0
       0  1  -0.666666667   0                     y - 2z /3 = 0
       0  0        5.10^-10        0                                 0 = 0
       0  0           0              1                                 0 = 1        and there is no solution at all !

Example3:   Invert Pascal's matrix:

      1  1  1   1   1
      1  2  3   4   5
      1  3  6  10 15
      1  4 10 20 35
      1  5 15 35 70

This is equivalent to solve 5 linear systems simultaneously. If you begin at register R01, you have to store the 50 numbers

        1  1  1   1  1    1  0  0  0  0                  R01  R06  R11  R16  R21  R26  R31  R36  R41  R46
        1  2  3   4  5    0  1  0  0  0                  R02  R07  R12  R17  R22  R27  R32  R37  R42  R47
        1  3  6  10 15  0  0  1  0  0         in      R03  R08  R13  R18  R23  R28  R33  R38  R43  R48        ( the right half is the unit matrix )
        1  4 10 20 35  0  0  0  1  0                  R04  R09  R14  R19  R24  R29  R34  R39  R44  R49
        1  5 15 35 70  0  0  0  0  1                  R05  R10  R15  R20  R25  R30  R35  R40  R45  R50

The control number of this rectangular array is 1.05005 and we can choose p = 0 because Pascal's matrix is non-singular:

    1.05005  ENTER^
           0      XEQ "LS"  yields  1 in X-register and in R00  ( the determinant of Pascal's matrices is always 1 )         ( execution time = 79s )

 The array is now:

       1  0  0  0  0   5  -10  10  -5   1
       0  1  0  0  0 -10  30 -35  19 -4
       0  0  1  0  0  10 -35  46 -27  6           ( the unit matrix is now on the left and the inverse matrix is on the right,
       0  0  0  1  0  -5   19 -27  17 -4              in registers R21 thru R50 )
       0  0  0  0  1   1   -4    6   -4   1
 

-This program can be used to invert a n x n matrix, but it requires 2n² registers and thus, it cannot invert a 13x13 matrix.
-See "Matrix Operations for the HP-41" if you want to find the inverse of up to a 17x17 matrix.
 

     b)  Program#2
 

-This variant below performs the same calculations as "LS"  but the coefficients are to be stored from register R01
-If flag F00 is clear the pivots smaller than E-7 ( line 60 ) is regarded as zero
-If flag F00 is set, only zero is regarded as zero
-Like "LS" , CF 01 = partial pivoting ,  SF 01 = no pivoting

-Line 60,  E-7  is a small number that identifies tiny elements.
-Line 108 is a three-byte GTO 01
 
 

01  LBL "LS2"    02  X<>Y   03  .1   04  %   05  +   06  STO N   07  FRC   08  STO O   09  *   10  1   11  STO 00          12  ST+ O   13  LASTX   14  %   15  +   16  +   17  LBL 01   18  STO M   19  FS? 01   20  GTO 04   21  INT   22  RCL O   23  FRC

24  +   25  ENTER^   26  ENTER^   27  CLX   28  LBL 02   29  RCL IND Z   30  ABS   31  X>Y?   32  STO Z   33  X>Y?   34  +   35  RDN   36  ISG Z   37  GTO 02   38  RCL M   39  ENTER^   40  FRC   41  R^   42  INT   43  +   44  X=Y?   45  GTO 04   46  LBL 03

47  RCL IND X   48  X<> IND Z   49  STO IND Y   50  ISG Y   51  RDN   52  ISG Y   53  GTO 03   54  RCL 00   55  CHS   56  STO 00   57  LBL 04   58  CLX   59  FC? 00   60   E-7    61  RCL IND M   62  ST* 00   63  ABS   64  X<=Y?   65  GTO 09   66  RCL M   67  LASTX   68  LBL 05   69  ST/ IND Y

70  ISG Y   71  GTO 05   72  RCL O   73  STO P    74  LBL 06   75  RCL M   76  ENTER^   77  FRC   78  RCL P   79  INT   80  +   81  X=Y?   82  GTO 08   83  RCL IND X   84  SIGN   85  RDN   86  LBL 07   87  RCL IND Y   88  LASTX   89  *   90  ST- IND Y   91  ISG Y   92  RDN

93  ISG Y   94  GTO 07   95  LBL 08   96  ISG P   97  GTO 06   98  LBL 09          99  RCL N 100  ST+ O 101  X^2 102  1 103  RCL M  104  + 105  X>Y? 106  CLX 107  ISG X 108  GTO 01  109  RCL 00 110  CLA 111  END

 
   ( 168 bytes / SIZE n.m+1 )
 
 

STACK	INPUTS	OUTPUTS
Z	/	( n.nnn )2
Y	n	/
X	m	det A

      n is the number of rows
      m is the number of columns

 Z-output is useful to retrieve n
 If n <= m  L-output = n².nm,nn

Example:   Of course, all the examples given after the "LS" listing can be solved ( but the first coefficient must be stored into R01 )  for instance:

          2x + 3y + 5z + 4t = 39
         -4x + 2y + z + 3t  = 15
          3x  -  y + 2z + 3t = 19
          5x + 7y - 3z + 2t = 18

-We store these 20 numbers:

              2  3  5  4  39                  R01  R05  R09  R13  R17
            -4  2  1  3  15        =        R02  R06  R10  R14  R18      respectively
             3 -1  2  3  19                  R03  R07  R11  R15  R19
             5  7 -3  2  18                  R04  R08  R12  R16  R20

-There are 4 rows and 5 columns,

      CF 00   CF 01
      4  ENTER^
      5  XEQ "LS2"  >>>>  Det A = 840 = R00  ( in 30 seconds )

-Registers R01 thru R16 now contains the unit matrix and registers R17 thru R20 contains the solution  x = 1 , y = 2 , z = 3 , t = 4
-Thus, the array has been changed into:

               1  0  0  0  1
               0  1  0  0  2              the solution is the last column
               0  0  1  0  3              of the new matrix.
               0  0  0  1  4

-"LS2" is slightly faster than "LS"
-When the program stops, R00 = det A
 

Note:   If you prefer to choose the "small" number that identifies the tiny elements:

  -Replace register M by register Q
  -Replace line 102 by ENTER^   SIGN
  -Replace lines 58 to 60 by  RCL M
  -Replace line 02 by STO M   X<> Z

-In this case, the inputs must be:
 

STACK	INPUTS	OUTPUTS
Z	n	( n.nnn )2
Y	m	/
X	p	det A

     n = the number of rows
     m = the number of columns
     p = the "small" number you have chosen.
 

     c1)  Program#3
 

-In this third variant, the matrix is the right-part of the array so that, when the program stops,
  the solution ( x₁ , x₂ , ............. , x_n ) is in registers  R01  R02  ..............  Rnn
-Here, the attempt to diagonalization starts by the lower right corner of the matrix.
-Flags F00 & F01 play the same role as above.

-Line 113 is a three-byte GTO 01
 
 

01  LBL "LS3"   02  X<>Y   03  .1   04  %   05  +   06  STO N   07  ST* Y   08  FRC   09  -   10  STO O   11  1   12  STO 00          13  LASTX   14  %   15  RCL Z   16  INT   17  +   18  LBL 01   19  STO M   20  FS? 01   21  GTO 04   22  INT   23  RCL O   24  FRC

25  +   26  ENTER^   27  ENTER^   28  CLX   29  LBL 02   30  RCL IND Z   31  ABS   32  X>Y?   33  STO Z   34  X>Y?   35  +   36  RDN   37  DSE Z   38  GTO 02   39  RCL M   40  ENTER^   41  FRC   42  R^   43  INT   44  +   45  X=Y?   46  GTO 04   47  LBL 03   48  RCL IND X

49  X<> IND Z   50  STO IND Y   51  DSE Y   52  RDN   53  DSE Y   54  GTO 03   55  RCL 00   56  CHS   57  STO 00   58  LBL 04   59  CLX   60  FC? 00   61   E-7   62  RCL IND M   63  ST* 00   64  ABS   65  X<=Y?   66  CLX   67  X=0?   68  STO 00   69  X=0?   70  GTO 09   71  RCL M   72  LASTX

73  LBL 05   74  ST/ IND Y   75  DSE Y   76  GTO 05   77  RCL O   78  STO P   79  LBL 06   80  RCL M   81  ENTER^   82  FRC   83  RCL P   84  INT   85  +   86  X=Y?   87  GTO 08   88  RCL IND X   89  SIGN   90  RDN   91  LBL 07   92  RCL IND Y   93  LASTX   94  *   95  ST- IND Y   96  DSE Y

97  RDN   98  DSE Y   99  GTO 07 100  LBL 08 101  DSE P 102  GTO 06        103  LBL 09 104  RCL N 105  ST- O 106  RCL O 107  RCL M 108  1 109  - 110  X<=Y? 111  CLX 112  DSE X 113  GTO 01 114  RCL 00 115  CLA 116  END

 
   ( 175 bytes / SIZE  n.m+1 )
 
 

STACK	INPUTS	OUTPUTS
T	/	n.nnn
Z	/	/
Y	n	/
X	m	det A

     n  =  number of rows
     m = number of columns

 T-output is useful to retrieve n

Example:

          2x + 3y + 5z + 4t = 39                                       39 = 2x + 3y + 5z + 4t
         -4x + 2y + z + 3t  = 15           is re-written          15 = -4x + 2y + z + 3t
          3x  -  y + 2z + 3t = 19                                       19 = 3x  -  y + 2z + 3t
          5x + 7y - 3z + 2t = 18                                       18 = 5x + 7y - 3z + 2t

  and we store these 20 numbers:

         39    2  3  5  4                        R01  R05  R09  R13  R17
         15  -4  2  1  3         in            R02  R06  R10  R14  R18      respectively
         19   3 -1  2  3                        R03  R07  R11  R15  R19
         18   5  7 -3  2                        R04  R08  R12  R16  R20

-There are 4 rows and 5 columns,

      CF 00   CF 01
      4  ENTER^
      5  XEQ "LS3"  >>>>  Det A = 840 = R00  ( in 30 seconds )

-Registers R05 thru R16 now contain the unit matrix and registers R01 thru R04 contain the solution  x = 1 , y = 2 , z = 3 , t = 4
-Thus, the array has been changed into:

           1  1  0  0  0
           2  0  1  0  0                the solution is the first column
           3  0  0  1  0                of the new matrix.
           4  0  0  0  1

-When the program stops, R00 = det A
-If you have to invert a matrix, the inverse will be in registers R01 thru Rn²

-Moreover, the programs "LS-" & "LS+" below become shorter if they call "LS3" instead of "LS2"

       c2)  Free42 Program

-Here is a version which can solve much larger linear systems with free42:

Data Registers:             R00 = Det M          R01 .... R09: temp                     ( Registers Rbb thru Ree are to be initialized before executing "LS42" )

                                    •  Rbb     = a₁₁    .....................    •  Ree-n+1 = a_1,m
                                    •  Rb+1   = a₂₁    .....................    •  Ree-n+2 = a_2,m
                                        .............................................

                                    •  Rb+n-1 = a_n1   ......................  •  Ree     =   a_n,m

Flags: /
Subroutines:  /

-Line 128 is a three-byte GTO 07
-  STO X  may be replaced by a TEXT0  or another NOP instruction
-With free42, line 63 ( E-7 ) could be replaced by a smaller number ( perhaps  E-20 )

01 LBL "LS42"  02 STO 08  03 X<>Y  04 STO 02  05 STO 04  06 *  07 X<>Y  08 STO 09  09 +  10 1  11 STO 00  12 -  13 STO 03  14 RCL 04  15  E5  16 /  17 +  18 LBL 07  19 STO 01            20 INT  21 ENTER  22 ENTER  23 RCL 04  24 STO 05  25 CLX  26 LBL 08  27 RCL IND Z  28 ABS

29 X>Y?  30 STO Z  31 X>Y?  32 +  33 RDN  34 DSE Z  35 STO X  36 DSE 05  37 GTO 08  38 RCL 01  39 ENTER  40 FRC  41 R^  42 INT  43 +  44 X=Y?  45 GTO 00  46 RCL 08  47 STO 05            48 RDN  49 LBL 09  50 RCL IND X  51 X<> IND Z  52 STO IND Y  53 DSE Y  54 RDN  55 DSE Y  56 STO X

57 DSE 05  58 GTO 09  59 RCL 00  60 CHS  61 STO 00  62 LBL 00  63  E-7  64 RCL IND 01  65 ST* 00  66 ABS  67 X<=Y?  68 CLX  69 X=0?  70 STO 00  71 X=0?  72 GTO 00  73 RCL 08            74 STO 05  75 RCL 01  76 LASTX  77 LBL 10  78 ST/ IND Y  79 DSE Y  80 STO X  81 DSE 05  82 GTO 10  83 RCL 02  84 STO 05

85 RCL 03  86 STO 06  87 LBL 11  88 RCL 01  89 ENTER  90 FRC  91 RCL 06  92 +  93 X=Y?  94 GTO 13  95 RCL 08            96 STO 07  97 CLX  98 RCL IND Y  99 SIGN 100 RDN 101 LBL 12 102 RCL IND Y 103 LASTX 104 * 105 ST- IND Y 106 DSE Y 107 RDN 108 DSE Y 109 STO X 110 DSE 07 111 GTO 12 112 LBL 13

113 DSE 06 114 DSE 05 115 GTO 11 116 LBL 00 117 RCL 02 118 ST- 03 119 RCL 01         120 1 121 - 122 DSE X 123 STO X 124 DSE 08 125 FS? 30 126 GTO 14 127 DSE 04 128 GTO 07 129 LBL 14 130 RCL 00 131 RCL 02 132 1 133 - 134  E3 135 / 136 RCL 09 137 .1 138 % 139 + 140 + 141 END

 
     ( 199 bytes / SIZE eee+1 )
 
 

STACK	INPUTS	OUTPUTS
Z	bbb > 9	det A
Y	n	det A
X	m	bbb.e'e'e'(1st column)

     n  =  number of rows
     m = number of columns

Example:

          2x + 3y + 5z + 4t = 39                                       39 = 2x + 3y + 5z + 4t
         -4x + 2y + z + 3t  = 15           is re-written          15 = -4x + 2y + z + 3t
          3x  -  y + 2z + 3t = 19                                       19 = 3x  -  y + 2z + 3t
          5x + 7y - 3z + 2t = 18                                       18 = 5x + 7y - 3z + 2t

  and if you choose bbb = 11,  store these 20 numbers:

         39    2  3  5  4                        R11  R15  R19  R23  R27
         15  -4  2  1  3         in            R12  R16  R20  R24  R28      respectively
         19   3 -1  2  3                        R13  R17  R21  R25  R29
         18   5  7 -3  2                        R14  R18  R22  R26  R30

-There are 4 rows and 5 columns,

      11  ENTER^
       4   ENTER^
       5  XEQ "LS42"  >>>>   11.014
                                 X<>Y  Det A = 840 = R00

-Registers R15 thru R30 now contain the unit matrix and registers R11 thru R14 contain the solution  x = 1 , y = 2 , z = 3 , t = 4
-Thus, the array has been changed into:

           1  1  0  0  0
           2  0  1  0  0                the solution is the first column
           3  0  0  1  0                of the new matrix.
           4  0  0  0  1

-When the program stops, R00 = det A
-If you want to invert a nxn matrix, the inverse will be in registers Rbb thru Rbb-1+n²   

     d)  Underdetermined Systems
 

-"LS" & "LS2" can find the solution(s) of an underdetermined system A.X = B ( if any )
-But the following programs compute the Moore-Penrose solution given by the formula:   X = ^tA ( A.^tA ) ^-1 B
 

        d-1)  Program#1
 

-Store the coefficients of the system in column order , starting with register R01
-Flags F00 & F01 are used like with "LS2"
- "LS2" cf paragraph b) above
  "TRN" transpose of a matrix                                     ( cf "Matrix Operations for the HP-41" )
  "M*"  multiplication of 2 matrices
  "TM*" multiplication of the transpose of a matrix by another one
  "MCO" copying a matrix
   are called as subroutines

-Lines 34 thru 44 may be replaced by

    RCL 00        /              1           +
    LASTX       +              +          REGMOVE
    +                 RCL 00    E3
    E6               X^2          /

-If you prefer to use "LS" instead of "LS2" replace lines 45 to 51 by

    1               E2         /         XEQ "LS"
    RCL 00      /           1         LASTX
    ST+ Y        +          +         FRC
    ST* Y       E3      E-7       ISG X
 
 

01  LBL "LS-"  02  STO 00        03  STO Z 04  DSE X 05  X<>Y 06  ST* Y 07  ST* Z 08   E2 09  / 10  + 11   E3 12  / 13  1 14  + 15  X<>Y

16  2 17  + 18  XEQ "TRN" 19  ENTER^ 20  STO Z 21  1 22  - 23  RDN 24  STO IND T 25  LASTX 26  XEQ "TM*" 27  X<>Y 28  ENTER^ 29  X<> 00        30  DSE X

31  * 32  1 33  + 34  STO Y 35  RCL 00 36  +  37   E3  38  /  39  +  40  RCL 00 41  X^2 42  1 43  + 44  XEQ "MCO" 45  RCL 00

46  RCL 00 47  1  48  +  49  XEQ "LS2" 50  X<> Z  51  SQRT 52  INT 53  STO 01 54  X^2 55  STO Y 56  LASTX 57  + 58  1 59  + 60  RCL IND X

61  STO T 62  SIGN 63  - 64  E3 65  / 66  + 67  ISG X 68  RCL 01        69   E5 70  / 71  + 72  1 73  XEQ "M*"  74  END

 
   ( 120 bytes / SIZE 2n.m-n+2 )
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	m	1.rrr,rr

     n = number of rows
     m = number of columns      ( n < m-1 )
    1.rrr,rr  is the control number of the solution ( in fact r = m-1 )

Example:    Let's find the Moore-Penrose solution of the following system:

      2x + 3y + 7z + 4t = 1                2  3   7   4  1           R01  R04  R07  R10  R13
      3x + 2y - 5z + 8t = 4     store    3  2  -5  8  4   into   R02  R05  R08  R11  R14   respectively
      4x + 5y + 6z + t  = 7                 4  5   6   1  7           R03  R06  R09  R12  R15

-There are 3 rows and 5 columns so:

      3  ENTER^
      5  XEQ "LS-"  >>>>   1.00404  ( in 49 seconds )

-The solution is in registers   R01 R02 R03 R04  namely:    x = 1.095088161 , y = 1.075776658 , z = -0.389378673 , t = -0.422963897

-When the program stops, R00 = det ( A.^tA )   [ in this example, R00 = 128628 ]
-If R00 = 0 or is very "small" ( suggesting A.^tA is singular ) , the results are probably meaningless.
 

        d-2)  Program#2
 

-This variant also calls   "TRN"  "M*"  "TM*"  "MCO"   as subroutines    ( cf "Matrix Operations for the HP-41" )
  but it uses "LS3" instead of "LS" or "LS2"
 
 

01  LBL "LS3-"  02  STO 00 03   E-2 04  + 05   E3 06   / 07  RCL Y 08  ST* Y 09  + 10  RCL 00        11  RCL Z 12  * 13  1

14  ST+ Z 15  + 16  XEQ "TRN" 17  ENTER^ 18  STO 00 19  X<> Z 20  FRC 21  ISG X 22  INT 23  1 24  + 25  XEQ "TM*"  26  LASTX

27  ENTER^ 28  ENTER^ 29  DSE X 30  STO T 31  * 32  1 33  + 34  X<> 00 35  STO IND 00 36  RDN 37  XEQ "LS3" 38  R^ 39  INT

40  ENTER^ 41  ENTER^ 42  X^2 43  + 44  1 45  + 46  RCL IND X  47  ENTER^ 48  INT 49  DSE X 50  R^ 51  1 52  %

53  + 54   E3 55  / 56  1 57  + 58  X<>Y 59  XEQ "M*" 60  1 61  XEQ "MCO" 62  END

 
   ( 106 bytes / SIZE 2n.m-n+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	m	1.rrr,rr

    n = number of rows
   m = number of columns      ( n < m-1 )
   1.rrr,rr  is the control number of the solution ( in fact r = m-1 )
 

Example:    Let's find the Moore-Penrose solution of the same system:

      1 = 2x + 3y + 7z + 4t                 1  2  3   7   4           R01  R04  R07  R10  R13
      4 = 3x + 2y - 5z + 8t      store    4  3  2  -5  8   into   R02  R05  R08  R11  R14   respectively
      7 = 4x + 5y + 6z + t                   7  4  5   6   1           R03  R06  R09  R12  R15

-There are 3 rows and 5 columns so:

      3  ENTER^
      5  XEQ "LS3-"  >>>>   1.00404  ( in 47 seconds )

-The solution is in registers   R01 R02 R03 R04  namely:    x = 1.095088162 , y = 1.075776659 , z = -0.389378674 , t = -0.422963896

-When the program stops, R00 = det ( A.^tA )   [ in this example, R00 = 128628 ]
-If R00 = 0 ( suggesting A.^tA is singular ) , the results are probably meaningless.
 

     e)  Overdetermined Systems
 

-An overdetermined system  A.X = B  has often no solution at all!
-But the following programs calculate the least-squares solution:
-It minimizes || A.X - B ||

Formula:   X = ( ^tA.A ) ^-1 ( ^tA.B )
 

        e-1)  Program#1
 

-Store the coefficients of the system in column order , starting with register R01
-Flags F00 & F01 are used like with "LS2"
- "LS2"  cf paragraph b) above
  "TM*" multiplication of the transpose of a matrix by another one             ( cf "Matrix Operations for the HP-41" )
  "MCO" copying a matrix
   are called as subroutines

-If you prefer to use "LS" instead of "LS2" replace lines 25 to 33 by

     RCL 00       E2         /         XEQ "LS"       INT
     ENTER^      /           1         LASTX           X^2
     DSE X        +          +         FRC                STO Y
     ST* Y         E3       E-7      ISG X             LASTX
 
 

01  LBL "LS+"  02  STO 00        03  RCL Y 04  ST* Y 05   E2 06  / 07  + 08  RCL X 09  RCL Z

10  - 11   E3 12  ST/ Z 13  / 14  1 15  ST+ Z 16  + 17  X<> Z 18  RCL 00

19  * 20  1 21  + 22  XEQ "TM*" 23  1 24  XEQ "MCO" 25  RCL 00  26  DSE X 27  RCL 00

28  XEQ "LS2"    29  X<> Z  30  INT  31  ENTER^ 32  ENTER^ 33  SQRT 34  + 35   E3 36  /

37  + 38  1 39  ST+ Y 40  XEQ "MCO" 41  END

 
    ( 78 bytes / SIZE n.m+m² -m+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	n	/
X	m	1.rrr,01

     n = number of rows
     m = number of columns      ( n >= m )
     1.rrr,01 is the control number of the solution ( in fact r = m-1 )

Example:       The system:

     5x + y + z  =   8             has no solution.
     4x - y + 2z =  13
       x + 2y - z =  -5
     7x - 4y + 5z =  32
     2x + 5y - 9z = -20

-We have to store these 20 coefficients into registers R01 thru R20 ( in column order )

       5   1   1   8                       R01  R06  R11  R16
       4  -1  2  13                      R02  R07  R12  R17          respectively
       1   2  -1 -5          =          R03  R08  R13  R18
       7  -4  5  32                      R04  R09  R14  R19
       2   5 -9 -20                      R05  R10  R15  R20

-There are 5 rows and 4 columns so

  5  ENTER^
  4  XEQ "LS+"  >>>>   1.00301  ( in 54 seconds )

-The "least-squares solution" is in registers    R01 R02 R03  namely:   x = 2.071207430  ,  y = -3.201238385  ,  z = 0.904024772

-When the program stops, R00 = det ( ^tA.A )   [ In this example  R00 = 55233 ]
-If det ( ^tA.A )  = 0 ( suggesting ^tA.A is singular ) , the results are probably meaningless.
 

        e-2)  Program#2
 

-This variant also calls  "TM*" and "MCO"   as subroutines    ( cf "Matrix Operations for the HP-41" )
  but it uses "LS3" instead of "LS" or "LS2"
 
 

01  LBL "LS3+"  02  STO 00        03   E3 04  / 05   E-5 06  + 07  1

08  + 09  RCL Y  10  * 11  ENTER^      12  FRC 13  1 14  ST+ Z

15  + 16  RCL 00 17  R^ 18  * 19  1 20  + 21  XEQ "TM*"

22  1 23  XEQ "MCO" 24  RCL 00 25  DSE X 26  RCL 00 27  XEQ "LS3" 28  R^

29  FRC             30  ISG X 31  END

 
    ( 59 bytes / SIZE n.m+m² -m+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	n	det
X	m	1.rrr

     n = number of rows
     m = number of columns      ( n >= m )
    1.rrr  is the control number of the solution ( in fact r = m-1 )

Example:      The system:

       8  =  5x + y + z              still has no solution.
      13 =  4x - y + 2z
      -5 =   x + 2y - z
      32 = 7x - 4y + 5z
     -20 =  2x + 5y - 9z

-We store these 20 coefficients into registers R01 thru R20   ( in column order )

          8    5   1   1                     R01  R06  R11  R16
         13   4  -1  2                     R02  R07  R12  R17
         -5   1   2  -1         =         R03  R08  R13  R18          respectively
         32   7  -4  5                     R04  R09  R14  R19
        -20   2   5  -9                    R05  R10  R15  R20

-There are 5 rows and 4 columns so

  5  ENTER^
  4  XEQ "LS3+"  >>>>   1.003  ( in 49 seconds )

-The "least-squares solution" is in registers    R01 R02 R03  namely:   x = 2.071207430  ,  y = -3.201238398  ,  z = 0.904024763

-When the program stops, Y = R00 = det ( ^tA.A )   [ In this example  R00 = 55233 ]
-If det ( ^tA.A )  = 0 or is very "small" ( suggesting ^tA.A is singular ) , the results are probably meaningless.
 
 

     f)  Tridiagonal Systems
 

-The following program solves a  nxn  linear system  A.X = B  where  A = [ a_i,j ]  is a tridiagonal matrix  i-e  a_i,j = 0  if  | i - j | > 1
-In other words, only the main diagonal and its 2 nearest neighbors have non-zero elements.

-Gaussian elimination is used without pivoting, so the method can fail, for instance if  a₁₁ = 0
  but practically, most of the problems leading to tridiagonal systems do not have these pitfalls.

-The components of the 3 diagonals are to be stored in column order into contiguous registers, starting with a register Rbb of your own choosing ( with bb  > 00 )
  followed by the  b_i 's ( the n elements of B )

-The determinant  det A  is also computed and stored in R00.
 

Data Registers:    R00 = det A  at the end                                              ( Registers Rbb thru Ree are to be initialized before executing "3DLS" )
 

                           •  Rbb     = a₁₁    •  Rbb+2 = a₁₂                                                                                                                          •  R3n+bb-2 = Ree-n+1 = b₁
                           •  Rbb+1 = a₂₁   •  Rbb+3 = a₂₂    •  Rbb+5 = a₂₃                                                                                                •  R3n+bb-1 = Ree-n+2 = b₂
                                                     •  Rbb+4 = a₃₂    •  Rbb+6 = a₃₃         ............
                                                                                •  Rbb+7 = a₄₃         ............                                                                                      .....................
                                                                                                                ............

                                                                                 •  R3n+bb-10 = a_n-3,n-2
                                                                                 •  R3n+bb-9 = a_n-2,n-2      •  R3n+bb-7 = a_n-2,n-1
                                                                                 •  R3n+bb-8 = a_n-1,n-2      •  R3n+bb-6 = a_n-1,n-1    •  R3n+bb-4 = a_n-1,n
                                                                                                                          •  R3n+bb-5 = a_n,n-1      •  R3n+bb-3 = a_n,n            •  R4n+bb-3 = Ree = b_n
 

     >>>>  When the program stops, the solutions  x₁ , .... , x_n  have replaced  b₁ , .... , b_n  in registers  R3n+bb-2  thru  R4n+bb-3

Flags: /
Subroutines: /
 
 

01  LBL "3DLS" 02  INT 03  LASTX 04  FRC 05   E3 06  * 07  RCL X 08  3 09  + 10  R^ 11  - 12  4 13  / 14  STO M 15  - 16  1

17  STO 00 18  + 19  STO N 20   E3 21  / 22  + 23  STO O 24  LBL 01 25  RCL O 26  RCL IND X 27  ST* 00 28  ST/ IND N 29  ISG O 30  ISG O 31  FS? 30 32  GTO 02

33  ST/ IND O 34  CLX 35  SIGN 36  + 37  RCL IND N 38  RCL IND Y 39  * 40  ISG N 41  CLX 42  ST- IND N 43  RCL IND O 44  LASTX 45  * 46  ISG O 47  ST- IND O 48  GTO 01

49  LBL 02 50  CLX 51  SIGN 52  - 53  INT 54  3 E-5 55  + 56  RCL M 57  DSE X 58  STO O 59  X<>Y 60  LBL 03 61  RCL IND X 62  RCL IND N 63  * 64  DSE N

65  ST- IND N  66  DSE Y 67  X<>Y 68  DSE O 69  GTO 03 70  RCL N 71  RCL M 72  + 73  DSE X 74   E3 75  ST/ Y 76  RDN 77  ST+ N 78  X<> N 79  CLA 80  END

 
   ( 132bytes / SIZE 4n+d-2 )
 
 

STACK	INPUTS	OUTPUTS
X	(bbb.eee)system	(bbb.eee)solution
L	/	n

      (bbb.eee)_system  is the control number of the system  ,  eee =  4n + bbb - 3  , bb > 00 ,  n = number of unknowns ,  n > 1
      (bbb.eee)_solution is the control number of the solution

Example:                                                                                   With  bb = 01,  store                                            into

   2.x + 5.y                                     = 2                                       2   5                           2                 R01   R03                                           R17
   3.x + 7.y + 4.z                            = 4                                       3   7   4                      4                 R02   R04   R06                                  R18
               y + 3.z + 7.t                    = 7                                            1   3   7                 7                          R05   R07   R09                         R19
                     2.z + 4.t + 6.u           = 1                                                 2   4   6            1                                    R08   R10   R12                R20
                              8.t +  u   + 7.v  = 5                                                      8   1   7       5                                             R11   R13   R15       R21
                                      9.u + 4.v  = 6                                                           9   4       6                                                       R14   R16       R22

-Then,   1.022  XEQ "3DLS"  >>>>   17.022    ( in 9 seconds )

        R17 = x = -16.47810408            R20 = t = -0.480422463
        R18 = y =     6.991241632          R21 = u =  0.112313241
        R19 = z =     1.123905204          R22 = v =  1.247295209

-We have also  R00 = det A = 3882

Notes:

-Synthetic registers M N O may be replaced by R01 R02 R03 but in this case, choose bb > 03
-The execution time is approximately proportional to n.
-With bb = 01, this program can solve a tridiagonal system of 80 equations in 80 unknowns in about  118 seconds,  provided it is executed from extended memory.
-There is a risk of OUT OF RANGE  line 27 if the determinant exceeds 9.999999999 E99 - especially for large n-values. Set flag F24 in this case.
-If you don't want to calculate det A , delete lines 27 and 17 , register R00 will be unused.
 

     g)  Pentadiagonal Systems
 

-"5DLS" hereunder solves a  nxn  linear system  A.X = B  where  A = [ a_i,j ]  is a pentadiagonal matrix  i-e  a_i,j = 0  if  | i - j | > 2
-Only the main diagonal and its 4 nearest neighbors have non-zero elements.

-Gaussian elimination is used without pivoting.

-The components of the 5 diagonals are to be stored in column order into contiguous registers, starting with a register Rbb of your own choosing ( with bb  > 00 )
  followed by the  b_i 's ( the n elements of B )

-The determinant  det A  is computed and stored in R00.
 

Data Registers:    R00 = det A  at the end                                                 ( Registers Rbb thru Ree are to be initialized before executing "5DLS" )
 

                 •  Rbb     = a₁₁    •  Rbb+3 = a₁₂   •  Rbb+7 = a₁₃                                                                                                     •  R5n+bb-6 = Ree-n+1 = b₁
                 •  Rbb+1 = a₂₁   •  Rbb+4 = a₂₂    •  Rbb+8 = a₂₃     •  Rbb+12 = a₂₄                                                                       •  R5n+bb-5 = Ree-n+2 = b₂
                 •  Rbb+2 = a₃₂   •  Rbb+5 = a₃₂    •  Rbb+9 = a₃₃     •  Rbb+13 = a₃₄
                                           •  Rbb+6 = a₄₂    •  Rbb+10 = a₄₃   •  Rbb+14 = a₄₄
                                                                      •  Rbb+11 = a₅₃   •  Rbb+15 = a₅₄
                                                                                                  •  Rbb+16 = a₆₄

                                                                                               .......................................                                                                 ...................................
 

                               •  R5n+bb-23 = a_n-5,n-3
                               •  R5n+bb-22 = a_n-4,n-3     •  R5n+bb-18 = a_n-4,n-2
                               •  R5n+bb-21 = a_n-3,n-3     •  R5n+bb-17 = a_n-3,n-2   •  R5n+bb-13 = a_n-3,n-1
                               •  R5n+bb-20 = a_n-2,n-3     •  R5n+bb-16 = a_n-2,n-2  •  R5n+bb-12 = a_n-2,n-1    •  R5n+bb-9 = a_n-2,n
                               •  R5n+bb-19 = a_n-1,n-3     •  R5n+bb-15 = a_n-1,n-2  •  R5n+bb-11 = a_n-1,n-1    •  R5n+bb-8 = a_n-1,n
                                                                         •  R5n+bb-14 = a_n,n-2     •  R5n+bb-10 = a_n,n-1      •  R5n+bb-7 = a_n,n              •  R6n+bb-7 = Ree = b_n
 

         >>>>  When the program stops, the solutions  x₁ , .... , x_n  have replaced  b₁ , .... , b_n  in registers  R5n+bb-6  thru  R6n+bb-7

Flags:  F06 & F07  are cleared by this routine.
Subroutines: /

-Line 96 is a three-byte  GTO 01
 
 

01  LBL "5DLS"   02  SF 06   03  CF 07   04  INT   05  LASTX         06  FRC   07   E3   08  *   09  RCL X   10  7   11  +   12  R^   13  -   14  6   15  /   16  STO M   17  -   18  STO N   19   E3   20  /   21  +   22  STO O   23  SIGN   24  STO 00   25  ST+ N   26  LBL 01   27  RCL N   28  RCL M   29  -   30  7   31  -   32  RCL O

33  X>Y?   34  SF 07   35  STO P    36  RCL IND X   37  ST* 00   38  ST/ IND N   39  ISG O   40  ISG O   41  ISG O   42  FC?C 06   43  ISG O   44  ST/ IND O   45  ISG P   46  RCL IND O   47  STO Z   48  RCL IND P   49  STO Q   50  *   51  ISG O   52  ST- IND O   53  X<> L   54  RCL IND N   55  *   56  ISG N   57  CLX   58  ST- IND N   59  ISG O   60  FS? 30   61  GTO 02   62  ISG P   63  X<> L   64  RCL IND P

65  R^   66  *   67  ST- IND O   68  ISG O   69  FC? 07   70  ISG O   71  R^   72  ST/ IND O   73  X<> Q   74  RCL IND O   75  *   76  ISG O   77  ST- IND O   78  X<> L   79  RCL IND P   80  *   81  ISG O   82  ST- IND O   83  R^   84  LASTX   85  *   86  ISG N   87  CLX   88  ST- IND N   89  DSE N   90  5   91  FS? 07   92  4   93  ST- O   94  FS? 07   95  SF 06   96  GTO 01

97  LBL 02   98  RCL O   99  INT 100  DSE X 101  RCL IND X 102  ST* 00 103  ST/ IND N 104  CLX 105  5 E-5 106  + 107  STO O 108  1 109  - 110  RCL IND X 111  RCL IND N 112  * 113  DSE N 114  ST- IND N 115  CLX 116  SIGN 117  - 118  RCL M 119  2 120  - 121  STO P  122  X<>Y 123  DSE O 124  LBL 03 125  ISG N 126  CLX 127  RCL IND X 128  RCL IND N

129  * 130  DSE N 131  RCL IND N 132  RCL IND O 133  * 134  + 135  DSE N 136  ST- IND N 137  CLX 138  SIGN 139  FS?C 07 140  ST+ Y 141  DSE O 142  CLX 143  DSE Y 144  X<>Y 145  DSE P 146  GTO 03 147  RCL N 148  RCL M 149  + 150  DSE X 151   E3 152  ST/ Y 153  RDN 154  ST+ N 155  X<> N 156  CLA 157  END

 
    ( 265 bytes / SIZE 6n+bb-6 )
 
 

STACK	INPUTS	OUTPUTS
X	(bbb.eee)system	(bbb.eee)solution
L	/	n

     (bbb.eee)_system  is the control number of the system  ,  eee =  6n + bbb - 7  ,  bb > 00 ,  n = number of unknowns ,  n > 2
     (bbb.eee)_solution is the control number of the solution

Example:

  7 x + 3 y + 4 z                                          =  1
     x + 8 y + 6 z +  t                                    =  2
  3 x + 2 y + 9 z + 2 t + 3 u                         =  3
           4 y + 4 z + 8 t + 5 u + 2 v                =  4
                    2 z + 3 t + 9 u + 3 v +  w        =  5
                             2 t + 3 u + 7 v + 2 w      =  6
                                        u + 6 v + 8 w      =  7

-With  bb = 07  store these 36 numbers

    7   3   4                           1                             R07   R10   R14                                             R36
    1   8   6   1                      2                             R08   R11   R15   R19                                    R37
    3   2   9   2   3                 3                             R09   R12   R16   R20   R24                           R38
         4   4   8   5   2            4               in                     R13   R17   R21   R25   R29                 R39                respectively
              2   3   9   3   1       5                                                R18   R22   R26   R30   R33        R40
                   2   3   7   2       6                                                          R23   R27   R31   R34        R41
                        1   6   8       7                                                                   R28   R32   R35        R42

  7.042  XEQ "5DLS"  >>>>     36.042    ( in 23 seconds )

      R36 = x = -0.023088805         R39 = t = 0.038788731       R42 = w = 0.364890723
      R37 = y =   0.069105035         R40 = u = 0.235429731
      R38 = z =   0.238576632         R41 = v = 0.640907414

-We have also  R00 = det A = 607264

Notes:

-Synthetic registers M N O P Q may be replaced by R01 R02 R03 R04 R05 but in this case, choose bb > 05
-"5DLS" cannot solve a 2x2 linear system: there must be at least 5 "diagonals".
-The execution time is roughly proportional to n.
-This program can solve a pentadiagonal system of 54 equations in 54 unknowns ( in 197 seconds ),  provided it is executed from extended memory.
-There is a risk of OUT OF RANGE  line 37 or 102 if the determinant exceeds 9.999999999 E99 - especially for large n-values. Set flag F24 in this case.
-If you don't want to calculate det A , delete lines 102 , 37 and 24 , register R00 will be unused.
 

2°) Non-Linear Systems
 

     a) One Equation in one Unknown:   f(x) = 0
 

       a-1) Secant Method
 

-Although it's obviously not a system of equations, it is logical to begin with this short routine.
 

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

                                          R01 = x , R02 = x' , R03 = f(x)
Flags: /
Subroutine:  A program which computes f(x) assuming x is in X-register upon entry.
 
 

01  LBL "SOLVE" 02  STO 01 03  X<>Y 04  STO 02 05  XEQ IND 00 06  STO 03 07  LBL 01 08  VIEW 01 09  RCL 01 10  XEQ IND 00 11  ENTER^ 12  ENTER^ 13  X<> 03 14  - 15  X#0? 16  / 17  RCL 02 18  RCL 01 19  STO 02 20  STO T 21  - 22  * 23  + 24  STO 01 25  X#Y? 26  GTO 01 27  END

 
( 43 bytes / SIZE 004 )
 
 

STACK	INPUTS	OUTPUTS
T	/	x
Z	/	x
Y	x'0	x
X	x0	x

  where  x₀  and  x'₀  are 2 initial guesses with  x₀ # x'₀

Example:   Find a solution near  x = 2 of  x³ - 4x + 1 = 0

   LBL"T"
   ENTER^
   X^2
   4
   -
   *
   1
   +
   RTN      "T"  ASTO 00  and if we choose  2 and 3 as initial guesses:

   2  ENTER^
   3  XEQ "SOLVE"  >>>>  the successive x-values are displayed and eventually:

   x = 1.860805853 = X-register = R01 = R02  ( f(x) ~ 10 ^-10 = R03 )

-The other solutions are  x' = 0.254101688  and  x" = -2.114907542   ( the last decimal should be a 1 )
 

Note:

-Always check f(x) in R03 because the iterations will also stop if  f(x) = f(x') , thus producing a possibly wrong result!
 

       a-2) Quadratic Interpolation
 

-The next step in interpolation methods is to use a parabola instead of a straight-line.
-Three starting-values are needed:   x₁ , x₂ , x₃   The next approximation x₄ is computed by the formula:

   x₄ = x₃ + 2.C / [ -B ± ( B² - 4 A.C )^1/2 ]     ( where the sign is chosen to produce the largest denominator )   with:

         A = ( x₂ - x₁ ) f(x₃) + ( x₁ - x₃ ) f(x₂) + ( x₃ - x₂ ) f(x₁)
         B = ( x₂ - x₁ ) ( 2 x₃ - x₂ - x₁ ) f(x₃) - ( x₁ - x₃ )² f(x₂) + ( x₃ - x₂ )² f(x₁)
         C = ( x₃ - x₂ ) ( x₂ - x₁ ) ( x₃ - x₁ ) f(x₁)

-If, however, the discriminant is negative, x₄ is computed by  x₄ = x₃ - B/(2A)
-This facilitates the determination of double-roots.

-And sometimes , the algorithm leads to a local extremum.
 
 

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

                                          R01 = x₃      R04 = f(x₃)       R07 thru R10: temp
                                          R02 = x₂      R05 = f(x₂)
                                          R03 = x₁      R06 = f(x₁)
Flags: /
Subroutine:  A program which computes f(x) assuming x is in X-register upon entry.
 
 

01  LBL "SOLVE2" 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  LBL 01 13  VIEW 01 14  RCL 01 15  XEQ IND 00 16  STO 04 17  RCL 02 18  RCL 03

19  - 20  * 21  ENTER^ 22  STO 07 23  RCL 01  24  RCL 03 25  - 26  STO 08 27  STO 10            28  ST* Z 29  RCL 05 30  * 31  ST* 08 32  - 33  RCL 02 34  RCL 01 35  - 36  STO 09

37  ST- 10 38  ST* Z 39  RCL 06 40  * 41  ST* 09 42  - 43  STO 06  44  * 45  RCL 07            46  RCL 10 47  * 48  RCL 08 49  - 50  RCL 09 51  + 52  2 53  / 54  STO 10

55  X^2 56  + 57  X<0? 58  GTO 02 59  SQRT 60  RCL 10  61  SIGN 62  * 63  RCL 10            64  + 65  GTO 03 66  LBL 02 67  RCL 10 68  CHS 69  RCL 06 70  LBL 03 71  X#0? 72  /

73  RCL 01 74  + 75  X<> 01 76  X<> 02 77  STO 03  78  RCL 04            79  X<> 05 80  STO 06 81  RCL 01 82  RCL 02 83  X#Y? 84  GTO 01 85  RCL 04 86  X<>Y 87  END

 
   ( 113 bytes / SIZE 011 )
 
 

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

 
Example:    f(x) = 12 x³ - 44 x² - 5 x + 100

  LBL "T"       12       -       -        +
  ENTER^      *        *       *      RTN
  STO Z         44       5     100

  "T"  ASTO 00

-If the starting-values are  -3  -2  -1

   -3  ENTER^
   -2  ENTER^
   -1  XEQ "SOLVE2"  >>>>  the successive approximations are displayed and finally:

     x = -1.333333333   RDN   f(x) = 0.000000050     ( in 17 seconds )

-If the starting-values are 1  2  3

    1  ENTER^
    2  ENTER^
    3     R/S    >>>>  the successive approximations are displayed and:

    x = 2.499981841   RDN    f(x) = 0     ( in 22 seconds )

-The double-root is actually  x = 2.5  and - as usual in such cases - the result is not very accurate.

Note:

-The rate of convergence is better than that of the secant method.
-Actually, it's the same as that of Newton's method: a quadratic convergence.
 

       a-3) Ridders' Method
 

-If a function is defined for a < x < b , interpolation methods may produce an approximation outside the interval [ a , b ] thus causing a DATA ERROR.
-Ridders' method avoids this pitfall: if a root is bracketed beween a and b , the successive approximations remain bracketed.
-Starting with 2 initial guesses x₁ & x₂ so that  f(x₁) . f(x₂) < 0 ,  the formulas are:

     x₃ = ( x₁ + x₂)/2
     x₄ = x₃ + ( x₃ - x₁ )  f(x₃)  { [ f(x₃) ]² -  f(x₁) f(x₂) } ^-1/2  sign [  f(x₁) -  f(x₂) ]
 

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

                                        R01 = x₂          R03 = x₁         R05 = x₃         R07 = x₄
                                        R02 = f(x₂)      R04 = f(x₁)      R06 = f(x₃)      R08 = E-8
Flags: /
Subroutine:  A program that computes f(x) assuming x is in X-register upon entry.

-Line 88  may be replaced by the M-Code routine  X#Y??  to avoid a possible infinite loop  ( cf "A few M-Code routines for the HP-41" )
 
 

01  LBL "RIDDERS" 02  STO 01              03  X<>Y 04  STO 03 05  XEQ IND 00 06  STO 04 07   E-8  08  STO 08 09  RCL 01 10  XEQ IND 00 11  STO 02 12  LBL 01 13  VIEW 01 14  RCL 01 15  RCL 03 16  + 17  2 18  / 19  STO 05 20  XEQ IND 00

21  ABS 22  RCL 08 23  X>Y? 24  GTO 05              25  LASTX 26  STO 06 27  RCL 04 28  RCL 02 29  - 30  SIGN 31  * 32  RCL 06 33  X^2 34  RCL 02 35  RCL 04 36  * 37  - 38  SQRT 39  / 40  RCL 05

41  RCL 03 42  - 43  * 44  RCL 05             45  + 46  STO 07 47  XEQ IND 00  48  ABS 49  RCL 08 50  X>Y? 51  GTO 06 52  RCL 06 53  LASTX 54  * 55  X<0? 56  GTO 02 57  RCL 04 58  LASTX 59  * 60  X<0?

61  GTO 03 62  LASTX 63  X<> 02 64  STO 04             65  RCL 07 66  X<> 01 67  STO 03  68  GTO 04 69  LBL 02 70  RCL 06 71  STO 04 72  RCL 05 73  STO 03 74  LBL 03 75  LASTX 76  STO 02 77  RCL 07 78  STO 01 79  LBL 04 80  RCL 01

81  RCL 03 82  X#Y?  83  GTO 01             84  GTO 08 85  LBL 05 86  RCL 05 87  GTO 07  88  LBL 06 89  RCL 07 90  LBL 07 91  STO 01 92  LASTX 93  STO 02 94  LBL 08 95  RCL 02 96  RCL 01 97  CLD 98  END

 
   ( 127 bytes / SIZE 009 )
 
 

STACK	INPUTS	OUTPUTS
Y	x1	f(x) ~ 0
X	x2	x

 
Example:    Find a solution of the equation  ( 1 - x² ) ^-1/2 - 4 = 0  between  0  and  0.99

  LBL "T"    LASTX    1/X       RTN
  X^2          -               4
  SIGN       SQRT       -

  'T'  ASTO 00

    0     ENTER^
  0.99  XEQ "RIDDERS"  >>>>  the successive x-approximations are displayed and finally:

    x = 0.968245837   RDN   f(x) ~ 0.000000003         --- execution time = 18 seconds ---
 

Note:

-The starting-values should verify  f(x₁) f(x₂) < 0
-Otherwise, you'll probably get a DATA ERROR message line 38 or 39
-However, it may also work ... sometimes!
 

       a-4)  Discrete Data ( 5 points )
 

-When the function is only known by tabulated values, one can use "SOLVE" & "LAGR" ( cf "Lagrange Interpolation for the HP-41" )  to solve  f(x) = 0
-The following routine is much less general:
-We assume that the abscissas are equally spaced  x_i+1 - x_i = h  and that there are 5 data points   y(x₁) = y₁ , ........... , y(x₅) = y₅

-In this case, the Lagrange interpolation polynomial may be written  y = y₃ + a n + b n² + c n³ + d n⁴   with  n = ( x - x₃ )/h
  and  f(x) = 0  is re-written   n = - (  y₃ + b n² + c n³ + d n⁴ ) / a
-So, this method may fail! For instance when a = 0
 

Data Registers:              R00 = h                             ( Registers R00 thru R05 are to be initialized before executing "SLVD" )

                                      •  R01 = y₁
                                      •  R02 = y₂
                                      •  R03 = y₃                            R06 to R09: temp      R10 = x₃
                                      •  R04 = y₄
                                      •  R05 = y₅
Flags: /
Subroutines: /
 
 

01  LBL "SLVD" 02  STO 10        03  X<>Y 04  STO 00 05  RCL 02 06  RCL 04 07  - 08  STO 06 09  LASTX 10  RCL 02 11  + 12  RCL 03 13  ST+ X 14  -

15  12 16  ST* 06 17  * 18  RCL 01        19  RCL 05 20  + 21  RCL Y 22  3 23  / 24  - 25  RCL 03 26  ST+ X 27  - 28  STO 09

29  - 30  STO 07 31  RCL 05 32  RCL 06        33  6 34  / 35  + 36  RCL 01 37  - 38  ST+ X 39  ST+ 06 40  STO 08 41  CLX 42  LBL 01

43  ENTER^ 44  ENTER^ 45  VIEW X 46  RCL 09        47  * 48  RCL 08 49  + 50  * 51  RCL 07 52  + 53  * 54  * 55  RCL 03 56  24

57  * 58  + 59  RCL 06        60  / 61  X#Y? 62  GTO 01 63  RCL 00 64  * 65  RCL 10 66  + 67  CLD 68  END

 
   ( 87 bytes / SIZE 011 )
 
 

STACK	INPUTS	OUTPUTS
Y	h	/
X	x3	x

 
Example:    Solve f(x) = 0  from the following data
 
 

x	-3	-1	1	3	5
y	-7	-3	-1	4	10

    -7   STO 01
    -3   STO 02
    -1   STO 03
     4    STO 04
    10   STO 05

     2   ENTER^
     1   XEQ "SLVD"  the successive n-values are displayed ( the first one is always 0 cf line 41 ) and eventually, it yields:

     x = 1.534588700   --- execution time = 17 seconds ---

-With the starting-values 1 and 2 , "SOLVE" & "LAGR" produce the same result in 1mn51s

Notes:

-To avoid an infinite loop, you can replace line 61 by the M-code routine  X#Y??  ( cf "A few M-Code Routines for the HP-41" )

-Or replace lines 61-62 by  X<>Y  RND  X<>Y  RND  X#Y?  GTO 01  LASTX
 and replace lines 41-42 by  CLX  +  LBL 01  LASTX

-The accuracy will be controlled by the display setting  (  FIX 9 gives the same result in the example above )
 

     b) 2 Equations in 2 Unknowns:   f(x,y) = 0 ; g(x,y) = 0
 
 

- "SXY" uses the Newton's iterative method - more exactly a generalization of the secant method to approximate the partial derivatives.
   Registers R00 thru R11 are used.
- It requires 2 initial guesses ( x , y ) and ( x' , y' ) which are to be stored in R01 , R02 and R03 , R04 respectively.
  You must choose x # x'  and y # y'   ( very important )

- You also have to load a subroutine that takes y in Y-register and x in X-register
   and calculates f(x;y) in Y-register and g(x;y) in X-register.

>>>   Thus, the stack has to be changed from:                  Y = y      into       Y' = f(x,y)             by your subroutine.
                                                                                        X = x                  X' = g(x,y)
   R11 and greater are available for this job.
- Then, you store the name of this subroutine into R00   ( global label of 6 characters maximum )
   and XEQ "SXY"

- The successive x-values are displayed  ( line 03 ) and when the program stops,

    x is in X-register and in R01
    y is in Y-register and in R02
    Z-register contains | f(x,y) | + | g(x,y) |

         f(x,y) is in R05
        g(x,y) is in R06

- To find another solution, re-initialize R01 thru R04 and R/S
 

Warning:

-The program stops when the approximate Jacobian determinant equals zero.
-This happens when   x = x' or y = y' but it may also happen by a stroke of bad luck,
   for instance if x converges much more quickly than y to the solution.

 >>>  That's why it's always wise to check the value of  | f | +| g |
 
 

01  LBL "SXY" 02  LBL 01 03  VIEW 01 04  RCL 02 05  RCL 03 06  XEQ IND 00 07  STO 08 08  X<>Y 09  STO 07 10  RCL 04 11  RCL 01 12  XEQ IND 00 13  STO 10 14  X<>Y 15  STO 09

16  RCL 02 17  RCL 01 18  XEQ IND 00 19  STO 06         20  ST- 08 21  ST- 10 22  X<>Y 23  STO 05 24  ST- 07 25  ST- 09 26  RCL 07 27  RCL 10 28  * 29  RCL 08 30  RCL 09

31  * 32  - 33  STO 11         34  X=0? 35  GTO 02  36  RCL 06 37  RCL 09 38  * 39  RCL 05 40  RCL 10 41  * 42  - 43  RCL 11 44  / 45  RCL 03

46  RCL 01 47  STO 03         48  - 49  * 50  ST+ 01 51  RCL 05  52  RCL 08 53  * 54  RCL 06 55  RCL 07 56  * 57  - 58  RCL 11 59  / 60  RCL 04

61  RCL 02 62  STO 04          63  - 64  * 65  ST+ 02 66  GTO 01 67  LBL 02 68  RCL 05  69  ABS 70  RCL 06 71  ABS 72  + 73  RCL 02 74  RCL 01 75  END

 
   ( 95 bytes / SIZE 012 )
 
 

STACK	INPUTS	OUTPUTS
Z	/	| f | + | g |
Y	/	y
X	/	x

 
Example:   Find x and y such that   x.y = 7 and  x² + y⁴ = 30

  First, we rewrite the system:

         x.y -  7  = 0
   x² + y⁴ - 30 = 0

 The subroutine may be ( let's call it "T" ):

01  LBL "T"
02  RCL Y
03  RCL Y
04  *
05  7
06  -
07  X<>Y
08  X^2
09  R^
10  X^2
11  X^2
12  +
13  30
14  -
15  RTN
 

 "T"  ASTO 00     and  with ( 2 ; 2 ) ( 3 ; 3 ) as initial guesses
  2  STO 01  STO 02
  3  STO 03  STO 04
  XEQ "SXY"              the successive x-values are displayed:

            2.000000000
            3.458333333
            3.369648334
            3.368074217
            3.368200504
            3.368200265
            3.368200265

 and the program stops with

      X = 3.368200265 = x                              ( execution time = 26s )
      Y = 2.078261222 = y
      Z = 0.000000011 = | f(x;y) | + | g(x;y) |    ( generally a small number if it's a "true" solution )
 
 

     c) 3 Equations in 3 Unknowns:   f(x,y,z) = 0 ; g(x,y,z) = 0 ; h(x,y,z) = 0
 

 - The same method leads to the following program for solving a system of 3 non-linear equations.
    Registers R00 thru R20 are used.
 - Two initial guesses  ( x , y , z )  and  ( x' , y' , z' ) are to be stored into  R01 , R02 , R03 and R04 , R05 , R06 respectively   ( x#x' , y#y' , z#z' )
 - You write a subroutine that takes   x in X-register ; y in Y-register ; z in Z-register
    and  computes  f(x,y,z) in Z-register ,  g(x,y,z) in  Y-register  ,  h(x,y,z)  in X-register

    The stack has to be change from:

             Z = z                Z' = f(x,y,z)
             Y = y     to       Y'= g(x,y,z)           ( registers R16 and greater are available for that purpose )
             X = x               X'= h(x,y,z)

 - Then, you store the name of this subroutine into R00 and XEQ "SXYZ"
    The successive x-values are displayed and when the program stops,

       x = X-register = R01
       y = Y-register = R02
       z = Z-register = R03

  and   T-register contains | f(x,y,z) | + | g(x,y,z) | + | h(x,y,z) |

       f(x,y,z) = R07
       g(x,y,z) = R08
       h(x,y,z) = R09

>> The recommendations for "SXY" still apply for "SXYZ"
 

Note:   The REGSWAP function is used in the following listing, but if you don't have an X-Functions module, you can:

     a) delete line 91
     b) replace lines 77 to 79 by  RCL 13   X<>16  STO 13  RCL 14  X<>17  STO 14  RCL 15  X<>18  STO 15
     c) replace lines  65 to 67 by  RCL 10   X<>16  STO 10  RCL 11  X<>17  STO 11  RCL 12  X<>18  STO 12
     d) replace lines  52 to 53 by  RCL 07   X<>16  STO 07  RCL 08  X<>17  STO 08  RCL 09  X<>18  STO 09
     e) delete line 49

In this case, SIZE 020 is enough.

-Line 89 is a three-byte  GTO 01
 
 

01  LBL  "SXYZ"   02  LBL 01   03  VIEW 01   04  RCL 03   05  RCL 02   06  RCL 04   07  XEQ IND 00   08  STO 09   09  RDN   10  STO 08   11  X<>Y   12  STO 07   13  RCL 03   14  RCL 05   15  RCL 01   16  XEQ IND 00   17  STO 12   18  RDN   19  STO 11   20  X<>Y   21  STO 10   22  RCL 06   23  RCL 02   24  RCL 01   25  XEQ IND 00   26  STO 15   27  RDN

28  STO 14   29  X<>Y   30  STO 13   31  RCL 03   32  RCL 02   33  RCL 01   34  XEQ IND 00   35  STO 18           36  ST- 09   37  ST- 12   38  ST- 15   39  RDN   40  STO 17   41  ST- 08   42  ST- 11   43  ST- 14   44  X<>Y   45  STO 16   46  ST- 07   47  ST- 10   48  ST- 13   49  CLX   50  XEQ 02   51  STO 19   52  7.016003   53  STO 20   54  XEQ 02

55  RCL 19   56  X=0?   57  GTO 03   58  /   59  RCL 04   60  RCL 01   61  STO 04           62  -   63  *   64  ST- 01   65  RCL 20    66  3   67  +   68  XEQ 02   69  RCL 19   70  /   71  RCL 05   72  RCL 02   73  STO 05   74  -   75  *   76  ST+ 02   77  RCL 20   78  6   79  +   80  XEQ 02   81  RCL 19

82  /   83  RCL 06   84  RCL 03   85  STO 06   86  -   87  *   88  ST- 03   89  GTO 01          90  LBL 02   91  REGSWAP    92  RCL 11   93  RCL 15   94  *   95  RCL 12   96  RCL 14   97  *   98  -   99  RCL 07 100  * 101  RCL 10 102  RCL 15 103  * 104  RCL 12 105  RCL 13 106  * 107  - 108  RCL 08

109  * 110  - 111  RCL 10 112  RCL 14 113  * 114  RCL 11 115  RCL 13         116  * 117  - 118  RCL 09  119  * 120  + 121  RTN 122  LBL 03 123  RCL 07 124  ABS 125  RCL 08 126  ABS 127  + 128  RCL 09 129  ABS 130  + 131  RCL 03 132  RCL 02 133  RCL 01 134  END

 
   ( 189 bytes / SIZE 021 )
 
 

STACK	INPUTS	OUTPUTS
T	/	| f | + | g | + | h |
Z	/	z
Y	/	y
X	/	x

 
Example:   Find a solution of the system:

      x.y² - z/y = 0
      x - y - z   = 0
      ln x + y.z = 0

 The subroutines for calculating  f , g , h  may be:

 01  LBL "T"
 02  STO 16
 03  X<>Y
 04  STO 17
 05  X^2
 06  *
 07  X<>Y
 08  STO 18
 09  RCL 17
 10  /
 11  -
 12  RCL 16
 13  RCL 17
 14  -
 15  RCL 18
 16  -
 17  RCL 17
 18  LASTX
 19  *
 20  RCL 16
 21  LN
 22  +
 23  RTN
 

   "T"  ASTO 00                                  and, if you choose  ( 2,2,2 ) and (1,1,1) as initial approximations:

    2  STO 01  STO 02  STO 03
    1  STO 04  STO 05  STO 06
    XEQ "SXYZ"

The successive x-values are displayed:

        2.000000000
        1.742625585
        0.896281748
        0.827584649
        0.865942155
        0.865417190
        0.865408830
        0.865408832
        0.865408832

Finally, it yields:

         x = 0.865408832 in X-register
         y = 0.639295476 in Y-register
         z = 0.226113356 in Z-register
        | f | + | g | + | h | = 10^-10        in T-register       ( the whole execution time is about 1mn37s )

More precisely,

         f(x,y,z) = -10^-10  in R07
         g(x,y,z) =   0       in R08
         h(x,y,z) =   0       in R09

>> Note that if you choose:

        1 STO 01 STO 02 STO 03
        2 STO 04 STO 05 STO 06  as initial guesses,

 the program stops after only 2 iterations and returns x = 1 ; y = 0.777777778 ; z = 0.222222222
 but you realize it's a completely wrong result because  T-register contains  0.492... ( not really very small! )

>> Good guesses are not always easy to find but in any case, always check the values of   f , g , h !
 

     d) N Equations in N Unknowns:   f₁( x₁,... ,x_n )  = 0 , .......... ,  f_n( x₁, ... ,x_n ) = 0
 

- Now, we deal with the general case of a system of nxn non-linear equations.
  These programs call "LS" , "LS2" or "LS3" as a subroutine ( synthetic version required ).
- Once again, the same ( quasi- ) Newton's method is used:
  Each iteration solves a linear system of n equations in n unknowns and that's why "LS" is needed.
- Unlike "SXY" and "SXYZ"  it's of course impossible to take the n variables from the stack and calculate the n functions in the stack if n > 4,
  therefore you'll have to key in n different subroutines for computing the f_i in the X-register with x₁ in R01 , ... , x_n in Rnn
- Synthetic registers M N O  and data registers R00 thru Rn²+4n are used.

- The successive x₁-values are displayed ( line 04 ) and when the program stops, | f₁ | + .... + | f_n | is in X-register
  and the solution  ( x₁ , .... , x_n )  is  in  ( R01 , ..... , Rnn )
 

        d-1)  Program#1
 

Data registers:        •  R00 = n                             ( Registers R00 thru R3n are to be initialized before executing "NLS" )

                                 •  R01 = x₁   •   Rn+1 = x'₁      •  R2n+1 = f₁  name
                                 •  R02 = x₂   •   Rn+2 = x'₂      •  R2n+2 = f₂  name                    R3n+1 thru Rn²+4n  contain the coefficients
                                    .....................................................................                          of the linear system that "LS" solves at every iteration.

                                 •  Rnn = x_n   •     R2n = x'_n      •  R3n  =  f_n  name

        where   (  x₁ , .... , x_n ) and ( x'₁ , .... , x'_n ) are the 2 initial guesses which must verify      x_i # x'_i   for  i = 1 , 2 , ... , n

Flag:   F01
Subroutines:     "LS"  ( synthetic version )  and  n programs that take  x₁ in R01 , ..... , x_n  in Rnn
                            and compute  f_i ( x₁ , .... , x_n )  in the X-register   i = 1 , ..... , n

-Replace line 02 by  SF 01  if you want to solve the successive linear systems without pivoting:
  it's of course faster but also more risky! ( it saves about 41 seconds in the example below )
-Line 107 is a three-byte GTO 01
 
 

01  LBL "NLS"   02  CF 01    03  LBL 01    04  VIEW 01   05  4.00301   06  RCL 00   07  ST+ Y   08  *   09  STO M   10  3.002   11  LASTX   12  *   13  STO N   14  LBL 02   15  RCL IND N   16  XEQ IND X   17  LBL 03   18  STO IND M   19  DSE M   20  GTO 03   21  RCL 00   22  ENTER^   23  X^2   24  +   25  DSE X

26  ST+ M   27  DSE N   28  GTO 02   29  3   30  RCL 00   31  ST+ Y   32  *   33  STO M   34  3.002   35  LASTX   36  *   37  STO N   38  LASTX   39  .1   40  %   41  +   42  -   43  STO O   44  LBL 04   45  RCL O   46  RCL 00   47  -   48  RCL IND X   49  X<> IND O   50  STO IND Y

51  LBL 05   52  RCL IND N   53  XEQ IND X   54  ST- IND M   55  DSE M   56  DSE N   57  GTO 05   58  RCL O   59  RCL 00   60  ST+ N   61  -   62  RCL IND X   63  X<> IND O   64  STO IND Y   65  DSE O   66  GTO 04   67  2.01   68  RCL 00   69  ST+ Y   70  *   71   E3   72  /   73  RCL N   74  +   75  1

76  +   77  0   78  XEQ "LS"    79  LASTX    80  FRC    81  ISG X    82  INT   83  STO 00   84  STO M   85  ST+ X   86  STO N   87  ST+ X   88  RCL 00   89  X^2   90  +   91  STO O   92  X<>Y   93  X=0?   94  GTO 07   95  LBL 06   96  RCL IND M   97  ENTER^   98  X<> IND N   99  - 100  RCL IND O

101  * 102  ST- IND M 103  DSE O 104  DSE N 105  DSE M 106  GTO 06 107  GTO 01  108  LBL 07 109  CLA 110  3.002 111  RCL 00 112  * 113  STO N 114  LBL 08 115  RCL IND N 116  XEQ IND X 117  ABS 118  ST+ M 119  DSE N 120  GTO 08 121  X<> M 122  CLA 123  CLD 124  END

 
   ( 219 bytes / SIZE n²+4n+1 )
 
 

STACK	INPUTS	OUTPUTS
X	/	Sum | fi |

 
Example:    Find a solution near ( 4 , 1 , 3 , 6 ) of the system:

       x₁ + x₂ + x₃ + x₄ - 16 = 0
             x₁ x₂ x₃ - 3 x₄ = 0
        4x₁² - x₂ x₃ x₄ - 40 = 0
          x₁ x₂ x₃  x₄ - 140 = 0

-Let's write the 4 subroutines, "F1" "F2" "F3" "F4"  ( any global labels, 6 characters maximum ):
 

 01  LBL "F1"
 02  RCL 01
 03  RCL 02
 04  RCL 03
 05  RCL 04
 06  +
 07  +
 08  +
 09  16
 10  -
 11  RTN
 12  LBL "F2"
 13  RCL 01
 14  RCL 02
 15  RCL 03
 16  *
 17  *
 18  RCL 04
 19  3
 20  *
 21  -
 22  RTN
 23  LBL "F3"
 24  RCL 01
 25  ST+ X
 26  X^2
 27  RCL 02
 28  RCL 03
 29  RCL 04
 30  *
 31  *
 32  -
 33  40
 34  -
 35  RTN
 36  LBL "F4"
 37  RCL 01
 38  RCL 02
 39  RCL 03
 40  RCL 04
 41  *
 42  *
 43  *
 44  140
 45  -
 46  RTN
 

  SIZE 033     ( or greater )
  4  STO 00    ( the number of unknowns )    and  if     ( 4 , 1 , 3 , 6 )  and  ( 4.1 , 1.1 , 3.1 , 6.1 )  are your 2 initial guesses:

     4  STO 01       4.1   STO 05     F1  ASTO 09
     1  STO 02       1.1   STO 06     F2  ASTO 10
     3  STO 03       3.1   STO 07     F3  ASTO 11
     6  STO 04       6.1   STO 08     F4  ASTO 12

   XEQ "NLS"    the successive x₁-approximations are displayed:

            4.000000000
            4.298102981      ( each iteration lasts about 58 seconds )
            4.266193620
            4.266545269
            4.266540470
            4.266540475
            4.266540475

and the program stops with  | f₁ | + | f₂ | + | f₃ | + | f₄ | = 10^-8  in X-register  and the solution in registers  R01 , R02 , R03 , R04

Thus,   x₁ =  4.266540475 , x₂ =  1.353632236 ,  x₃ = 3.548526779 ,  x₄ = 6.831300511      ( the whole execution time is about 6mn33s )

>>>>  This system has several solutions. Another one is  ( 4.266540475 , -0.2552286465 , 18.81998868 , -6.831300511 )
 

Notes:
 

-I've tested this program to solve a 7x7 non-linear system of ( relatively ) simple equations: every iteration requires about 3 minutes ( see the exercise below ).
-Although I think it's rare, it might happen that this program runs for ever because of successive x-values that are very close together:
-In order to avoid this infinite loop:

    Add  E-8 ( or another "small" number )  X<Y?  after line 106
    Add  ABS  +  after line 102
    Add  CLX  after line 94

-This modification also reduces execution time:
  In the example above, the last ( unuseful ) iteration is avoided.

-There is another way to avoid the unuseful iteration ( when  x_i = x'_i for some i )

    Replace lines 97 to 107 by

    ENTER^  ENTER^  ENTER^  X<> IND N  -  RCL IND O  *  -  STO IND M  X=Y?  SF 01  DSE O  DSE N  DSE M  GTO 06  FC?C 01  GTO 01

-The same warnings mentioned for "SXY" and "SXYZ" still apply to "NLS"

*** If you want to use "LS2" instead of "LS", it's quite possible, but a little more complicated. For instance:

    Replace lines 67 to 91 by

    CLX  RCL 00  X^2  E3  /  +  E3  /  RCL M  +  1  +  REGSWAP  RCL 00  RCL 00  LASTX  +  XEQ "LS2"
    RCL Z  INT  ENTER^  SQRT  STO 00  STO M  STO N  ST+ N  +  STO O  E3  /  ISG X  RCL 00  3  *
    ST+ O  +  E3  /  1  +  REGSWAP

    and add   SF 00   after line 01

-Actually, this variant is slightly faster: in the example above, execution time = 6m22s instead of 6m33s
 

Exercise:     Find a solution of the following system of 7 equations in 7 unknowns:

     x³ + y²z + t.u - v² - w² = 0
     x²y - z.t.u² + x.v - w³ = 0
     x + y + z + t - u - v - w = 0
     x³ - y.z.t + t.u.v - w² = 0
     x.y⁴ - 2y.z³ - t.u.v²w = 0
     x + y.z + t.u - v.w² = 0
     x.y - y.z.t.u.v + w - 1 = 0

Answer:    With the initial guesses   x = y = z = t = u = v = w = 1  &  x' = y' = z' = t' = u' = v' = w' = 2  ,  "NLS" yields after about 27 minutes:

        x = 1.200271274         u = 1.280018254
        y = 1.548046396         v = 1.698155509
        z = 1.011876841         w = 1.463999879
        t = 0.681979132

-With these values,  Sum _i=1,....,7  | f_i |  =  4 10 ^-9        ( Simply press XEQ 07 if you ever lose this result )
 

        d-2)  Program#2
 

-Use the following variant if you prefer "LS3" to "LS" or "LS2"
 

Data registers:        •  R00 = n                             ( Registers R00 thru R3n are to be initialized before executing "NLS3" )

                                 •  R01 = x₁   •   Rn+1 = x'₁      •  R2n+1 = f₁  name
                                 •  R02 = x₂   •   Rn+2 = x'₂      •  R2n+2 = f₂  name
                                    .....................................................................                                       R3n+1 to Rn²+4n:  temp

                                 •  Rnn = x_n   •     R2n = x'_n      •  R3n  =  f_n  name

        where   (  x₁ , .... , x_n ) and ( x'₁ , .... , x'_n ) are the 2 initial guesses which must verify      x_i # x'_i   for  i = 1 , 2 , ... , n

Flags:  F00 - F01
Subroutines:     "LS3"  and  n routines that take  x₁ in R01 , ..... , x_n  in Rnn  and compute  f_i ( x₁ , .... , x_n )  in X-register   i = 1 , ..... , n

-Line 127 is a three-byte  GTO 01
 
 

01  LBL "NLS3"   02  CF 01   03  LBL 01   04  VIEW 01   05  4.00301   06  RCL 00   07  ST+ Y   08  *   09  STO M   10  3.002   11  LASTX   12  *   13  STO N   14  STO O   15  LBL 02   16  RCL IND O   17  XEQ IND X   18  LBL 03   19  STO IND M   20  DSE M   21  GTO 03   22  RCL 00   23  ENTER^   24  X^2   25  +   26  DSE X   27  ST+ M   28  DSE O   29  GTO 02

30  4   31  RCL 00   32  ST+ Y   33  *   34  STO M   35  RCL 00   36  LASTX   37  .1   38  %   39  +   40  +   41  STO O   42  LBL 04   43  RCL O   44  RCL 00   45  -   46  RCL IND X   47  X<> IND O   48  STO IND Y   49  LBL 05   50  RCL IND N   51  XEQ IND X   52  ST- IND M   53  DSE M   54  DSE N   55  GTO 05   56  RCL O   57  RCL 00   58  ST+ N

59  -   60  RCL IND X   61  X<> IND O   62  STO IND Y   63  DSE O   64  GTO 04   65  CLX   66  RCL 00   67  ST- M   68  X^2   69   E3   70  /   71  +   72   E3   73  /   74  RCL M   75  +   76  1   77  +   78  REGSWAP   79  RCL 00   80  RCL 00   81  LASTX   82  +   83  XEQ "LS3"   84  R^   85  INT   86  STO 00   87  STO M

88  STO N   89  ST+ N   90  STO O   91  ENTER^   92  X^2   93  +   94   E3   95  /   96  ISG X   97  RCL 00   98  3   99  * 100  ST+ O 101  + 102   E3 103  / 104  1 105  + 106  REGSWAP 107  X<>Y 108  X=0? 109  GTO 07 110  CLX  111  LBL 06 112  RCL IND M 113  ENTER^ 114  X<> IND N 115  - 116  RCL IND O

117  * 118  ST- IND M 119  ABS 120  + 121  DSE O 122  DSE N 123  DSE M 124  GTO 06 125   E-8 126  X<Y? 127  GTO 01  128  LBL 07 129  CLA 130  3.002 131  RCL 00 132  * 133  STO N 134  LBL 08 135  RCL IND N 136  XEQ IND X 137  ABS 138  ST+ M 139  DSE N 140  GTO 08 141  X<> M 142  CLA 143  CLD 144  END

 
   ( 241 bytes / SIZE n²+4n+1 )
 
 

STACK	INPUTS	OUTPUTS
X	/	Sum | fi |

-The same example gives - in 5mn37s -   | f₁ | + | f₂ | + | f₃ | + | f₄ | = 2 10^-8  in X-register  and the solution in registers  R01 , R02 , R03 , R04

  namely:   x₁ =  4.266540475 , x₂ =  1.353632241 ,  x₃ = 3.548526764 ,  x₄ = 6.831300511
 

3°) A short In/Out routine:
 

 If flag F02 is clear, "IO" stores data.
 If flag F02 is set, "IO" displays data.

 RCL d  may be replaced by  RCLFLAG
 STO d  may be replaced by  STO FLAG

-The append character is denoted  ~
 
 

01  LBL "IO"  02  CF 22 03  CF 23 04  1 05  LBL 01       06  RCL d  07  FIX 0 08  CF 29 09  " "  10  ARCL Y 11  STO d

12  RDN 13  FS? 02 14  GTO 03 15  "~?"  16  PROMPT  17  FC?C 22 18  FS? 23 19  FS? 30 20  GTO 02       21  FS? 23 22  ENTER^

23  FS?C 23 24  ASTO X 25  STO IND Z  26  CLX 27  SIGN 28  + 29  ISG Y 30  CLA 31  GTO 01       32  LBL 02 33  DSE X

34  X<>Y 35  - 36  CHS 37  DSE L 38  LASTX        39   E3 40  / 41  + 42  CLA 43  AOFF 44  RTN

45  LBL 03 46  "~="  47  ARCL IND Y 48  AVIEW 40  1 50  + 51  ISG Y 52  GTO 01 53  CLA 54  END

 
   ( 91 bytes )
 
 

STACK	CF02 INPUTS	OUTPUTS
X	bbb	bbb.eee

STACK	SF02 INPUTS	OUTPUTS
X	bbb.eee	/

 
Example:        Suppose you have to store  7 , 12 , 41 , "AAA" , "HP41"   into registers  R04 , R05 , ... etc ...

-Key in      CF 02

     4    XEQ "IO"    the HP-41  displays     1?
     7        R/S          the HP-41  displays     2?
    12       R/S          the HP-41  displays     3?
    41       R/S          the HP-41  displays     4?       switch to alpha mode and
 "AAA"   R/S          the HP-41  displays     5?
 "HP41"  R/S          the HP-41  displays     6?

 simply press    R/S     and   you'll get   4.008     in X-register:  this is the control number  bbb.eee   of your data in registers R04 thru R08

-To obtain the control number needed for "LS", simply add  r.10^-5 to the previous result

   ( r = the number of rows of the matrix = the number of equations of the system )

>>>  If you store a number like "PI" set flag F22 before pressing the R/S key

-If you want to view the contents of registers R04 to R08  press  SF 02    4.008    XEQ "IO"
   and the HP-41 will display successively:

        1 = 7
        2 = 12
        3 = 41
        4 = AAA
        5 = HP41

 ( if you set flag F21 , the calculator will stop at each  AVIEW )

>>>You can also use "IO" to store data into registers  R04 , R07 , R10 , ... etc ...

   Key in   CF 02   4.00003  XEQ "IO"  and so on , but in this case, simply ignore the last prompt because the control number would be completely wrong!

-Likewise, to view registers   R04 , R07 , R10 , R13 , R16   for instance , key in   SF 02   4.01603  XEQ "IO"   ...
 

References:

[1]  Francis Scheid - "Numerical Analysis" - ( McGraw-Hill )  ISBN:  07-055197-9
[2]  Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4