hp41programs

Systems Linear and non-Linear Systems for the HP-41
 

Overview

 1°) Linear Systems

    a) Program#1
    b) Program#2
    c) Program#3
    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"
 

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 2n2 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 = n2.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.
 

     c)  Program#3
 

-In this third variant, the matrix is the right-part of the array so that, when the program stops,
  the solution ( x1 , x2 , ............. , xn ) 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 Rn2

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

     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+m2 -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+m2 -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 = [ ai,j ]  is a tridiagonal matrix  i-e  ai,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  a11 = 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  bi '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     = a11    •  Rbb+2 = a12                                                                                                                          •  R3n+bb-2 = Ree-n+1 = b1
                           •  Rbb+1 = a21   •  Rbb+3 = a22    •  Rbb+5 = a23                                                                                                •  R3n+bb-1 = Ree-n+2 = b2
                                                     •  Rbb+4 = a32    •  Rbb+6 = a33         ............
                                                                                •  Rbb+7 = a43         ............                                                                                      .....................
                                                                                                                ............

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

     >>>>  When the program stops, the solutions  x1 , .... , xn  have replaced  b1 , .... , bn  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 = [ ai,j ]  is a pentadiagonal matrix  i-e  ai,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  bi '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     = a11    •  Rbb+3 = a12   •  Rbb+7 = a13                                                                                                     •  R5n+bb-6 = Ree-n+1 = b1
                 •  Rbb+1 = a21   •  Rbb+4 = a22    •  Rbb+8 = a23     •  Rbb+12 = a24                                                                       •  R5n+bb-5 = Ree-n+2 = b2
                 •  Rbb+2 = a32   •  Rbb+5 = a32    •  Rbb+9 = a33     •  Rbb+13 = a34
                                           •  Rbb+6 = a42    •  Rbb+10 = a43   •  Rbb+14 = a44
                                                                      •  Rbb+11 = a53   •  Rbb+15 = a54
                                                                                                  •  Rbb+16 = a64

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

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

         >>>>  When the program stops, the solutions  x1 , .... , xn  have replaced  b1 , .... , bn  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  x0  and  x'0  are 2 initial guesses with  x0 # x'0

Example:   Find a solution near  x = 2 of  x3 - 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:   x1 , x2 , x3   The next approximation x4 is computed by the formula:

   x4 = x3 + 2.C / [ -B ± ( B2 - 4 A.C )1/2 ]     ( where the sign is chosen to produce the largest denominator )   with:

         A = ( x2 - x1 ) f(x3) + ( x1 - x3 ) f(x2) + ( x3 - x2 ) f(x1)
         B = ( x2 - x1 ) ( 2 x3 - x2 - x1 ) f(x3) - ( x1 - x3 )2 f(x2) + ( x3 - x2 )2 f(x1)
         C = ( x3 - x2 ) ( x2 - x1 ) ( x3 - x1 ) f(x1)

-If, however, the discriminant is negative, x4 is computed by  x4 = x3 - 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 = x3      R04 = f(x3)       R07 thru R10: temp
                                          R02 = x2      R05 = f(x2)
                                          R03 = x1      R06 = f(x1)
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 x3 - 44 x2 - 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 x1 & x2 so that  f(x1) . f(x2) < 0 ,  the formulas are:

     x3 = ( x1 + x2)/2
     x4 = x3 + ( x3 - x1 )  f(x3)  { [ f(x3) ]2 -  f(x1) f(x2) } -1/2  sign [  f(x1) -  f(x2) ]
 

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

                                        R01 = x2          R03 = x1         R05 = x3         R07 = x4
                                        R02 = f(x2)      R04 = f(x1)      R06 = f(x3)      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 - x2 ) -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(x1) f(x2) < 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  xi+1 - xi = h  and that there are 5 data points   y(x1) = y1 , ........... , y(x5) = y5

-In this case, the Lagrange interpolation polynomial may be written  y = y3 + a n + b n2 + c n3 + d n4   with  n = ( x - x3 )/h
  and  f(x) = 0  is re-written   n = - (  y3 + b n2 + c n3 + d n4 ) / 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 = y1
                                      •  R02 = y2
                                      •  R03 = y3                            R06 to R09: temp      R10 = x3
                                      •  R04 = y4
                                      •  R05 = y5
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  x2 + y4 = 30

  First, we rewrite the system:

         x.y -  7  = 0
   x2 + y4 - 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.y2 - 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:   f1( x1,... ,xn )  = 0 , .......... ,  fn( x1, ... ,xn ) = 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 fi in the X-register with x1 in R01 , ... , xn in Rnn
- Synthetic registers M N O  and data registers R00 thru Rn2+4n are used.

- The successive x1-values are displayed ( line 04 ) and when the program stops, | f1 | + .... + | fn | is in X-register
  and the solution  ( x1 , .... , xn )  is  in  ( R01 , ..... , Rnn )
 

        d-1)  Program#1
 

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

                                 •  R01 = x1   •   Rn+1 = x'1      •  R2n+1 = f1  name
                                 •  R02 = x2   •   Rn+2 = x'2      •  R2n+2 = f2  name                    R3n+1 thru Rn2+4n  contain the coefficients
                                    .....................................................................                          of the linear system that "LS" solves at every iteration.

                                 •  Rnn = xn   •     R2n = x'n      •  R3n  =  fn  name

        where   (  x1 , .... , xn ) and ( x'1 , .... , x'n ) are the 2 initial guesses which must verify      xi # x'i   for  i = 1 , 2 , ... , n

Flag:   F01
Subroutines:     "LS"  ( synthetic version )  and  n programs that take  x1 in R01 , ..... , xn  in Rnn
                            and compute  fi ( x1 , .... , xn )  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 n2+4n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           X             /       Sum | fi

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

       x1 + x2 + x3 + x4 - 16 = 0
             x1 x2 x3 - 3 x4 = 0
        4x12 - x2 x3 x4 - 40 = 0
          x1 x2 x3  x4 - 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 x1-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  | f1 | + | f2 | + | f3 | + | f4 | = 10-8  in X-register  and the solution in registers  R01 , R02 , R03 , R04

Thus,   x1 =  4.266540475 , x2 =  1.353632236 ,  x3 = 3.548526779 ,  x4 = 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  xi = 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:

     x3 + y2z + t.u - v2 - w2 = 0
     x2y - z.t.u2 + x.v - w3 = 0
     x + y + z + t - u - v - w = 0
     x3 - y.z.t + t.u.v - w2 = 0
     x.y4 - 2y.z3 - t.u.v2w = 0
     x + y.z + t.u - v.w2 = 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  | fi |  =  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 = x1   •   Rn+1 = x'1      •  R2n+1 = f1  name
                                 •  R02 = x2   •   Rn+2 = x'2      •  R2n+2 = f2  name
                                    .....................................................................                                       R3n+1 to Rn2+4n:  temp

                                 •  Rnn = xn   •     R2n = x'n      •  R3n  =  fn  name

        where   (  x1 , .... , xn ) and ( x'1 , .... , x'n ) are the 2 initial guesses which must verify      xi # x'i   for  i = 1 , 2 , ... , n

Flags:  F00 - F01
Subroutines:     "LS3"  and  n routines that take  x1 in R01 , ..... , xn  in Rnn  and compute  fi ( x1 , .... , xn )  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 n2+4n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           X             /       Sum | fi

-The same example gives - in 5mn37s -   | f1 | + | f2 | + | f3 | + | f4 | = 2 10-8  in X-register  and the solution in registers  R01 , R02 , R03 , R04

  namely:   x1 =  4.266540475 , x2 =  1.353632241 ,  x3 = 3.548526764 ,  x4 = 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