The Museum of HP Calculators


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

This program is Copyright © 2002-2007 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

Overview

 1°) Linear Systems

    a) Program#1
    b) Program#2
    c) Underdetermined Systems
    d) Overdetermined Systems
    e) Tridiagonal Systems  ( new )
    f)  Pentadiagonal Systems ( new )

 2°) Non-Linear Systems

    a) 1 Equation
    b) 2 Equations in 2 unknowns
    c) 3 Equations in 3 unknowns
    d) N Equations in N unknowns

 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 for 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.
 

  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                    ( synthetic )
  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                  ( a synthetic three-byte GTO )
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        into        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  you key in:   1.02004 ENTER^
                                      CF01               10-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         into          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^
                    10-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         meaning         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         into      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 = 1mn19s )

                                      1  0  0  0  0   5  -10  10  -5   1
                                      0  1  0  0  0 -10  30 -35  19 -4
 The array is now:           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 can't 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
 

  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                    ( or another "small" number to identify the "tiny" elements )
  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               ( synthetic )
  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                  ( a three-byte GTO )
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

 where n is the number of rows
          m  --------------- 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        into        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

 where n is the number of rows
          m  --------------- columns
  and  p is the "small" number you have choosen.
 

     c)  Underdetermined Systems
 

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

-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
 

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                           Lines 34 thru 44 may be replaced by
35  RCL 00
36  +                                    RCL 00        /              1           +
37   E3                                 LASTX       +              +          REGMOVE
38  /                                     +                 RCL 00    E3
39  +                                    E6               X^2          /
40  RCL 00
41  X^2
42  1
43  +
44  XEQ "MCO"
45  RCL 00                       if you prefer to use "LS" instead of "LS2" replace lines 45 to 51 by
46  RCL 00
47  1                                  1               E2         /         XEQ "LS"
48  +                                 RCL 00      /           1         LASTX
49  XEQ "LS2"                 ST+ Y        +          +         FRC
50  X<> Z                         ST* Y       E3      E-7       ISG X
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

 where n is the number of rows
          m  --------------- 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     We 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)  Overdetermined Systems
 

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

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

-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
 

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                       if you prefer to use "LS" instead of "LS2" replace lines 25 to 33 by
26  DSE X
27  RCL 00                       RCL 00       E2         /         XEQ "LS"       INT
28  XEQ "LS2"                 ENTER^      /           1         LASTX           X^2
29  X<> Z                         DSE X        +          +         FRC                STO Y
30  INT                             ST* Y         E3       E-7      ISG X             LASTX
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

 where n is the number of rows
          m  --------------- 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         into           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 or is very "small" ( suggesting tA.A is singular ) , the results are probably meaningless.
 

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

   where  (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  --------------------------  solution
 

Example:

   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      if you choose  bb = 01               1   3   7                 7                          R05   R07   R09                        R19
                     2.z + 4.t + 6.u           = 1      store these 22 elements                   2   4   6            1          into                    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 )

  and           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.
 

     f)  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: /
 

  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                    ( synthetic )
  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                   a three-byte GTO
  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                    ( synthetic )
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

   where  (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  --------------------------  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               into                     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 )

    and      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:   f(x) = 0
 

-This program uses the secant method to solve f(x) = 0
-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!
 

     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,
   the solution x is in the X-register and in R01
   ----------- y -------  Y--------------- R02
                                    Z-register contains | f(x;y) | + | g(x;y) |                       f(x;y) is in R05
                                                                                                            and  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 verify 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 subroutines can 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

                                                                 Z = z                  Z' = f(x;y;z)
>> The stack has to be changed from:       Y = y     into       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,

     the solution   x is in X-register and in R01
                         y ---- Y--------------- R02
                         z -----Z--------------- R03
                           and  T-register contains | f(x;y;z) | + | g(x;y;z) | + | h(x;y;z) |         f(x;y;z) in R07
                                                                                                                            g(x;y;z) in R08
                                                                                                                            h(x;y;z) in 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) ------------  65 to 67 by  RCL 10   X<>16  STO 10  RCL 11  X<>17  STO 11  RCL 12  X<>18  STO 12
                         d) ------------  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.
 

  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              ( a synthetic three-byte GTO )
  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 |
           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
 

Then   "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
     and  | 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 gives 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 verify 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.
  This program calls "LS" 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 )
 

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                    R3n+1 thru Rn2+4n  contains the partial derivatives
                                   R02 = x2      Rn+2 = x'2        R2n+2 = f2  name                    of the n functions and the values of these n functions too,
                                    .....................................................................                         i-e 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
 

  01  LBL "NLS"
  02  CF 01                       replace this line by  SF 01  if you want to solve the successive linear systems without pivoting:
  03  LBL 01                     it's of course faster but also more risky! ( it saves about 41 seconds in the example below )
  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"               the value of n in R00 is lost when "LS" is executed
  79  LASTX                   but then, the control number bb.eeenn
  80  FRC                        is recalled from L-register
  81  ISG X                     and n is stored into R00 again
  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                   ( a synthetic three-byte GTO )
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
 

Then      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 in 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 )
 

3°) A short In/Out routine:

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

01  LBL "IO"
02  CF 22
03  CF 23
04  1
05  LBL 01
06  RCL d           ( or RCLFLAG )
07  FIX 0
08  CF 29
09  " "                  ( one space )
10  ARCL Y
11  STO d          ( or STOFLAG )
12  RDN
13  FS? 02
14  GTO 03
15  "~?"               ( I mean  "append ?" )
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  "~="                  ( "append =" )
47  ARCL IND Y
48  AVIEW
40  1
50  +
51  ISG Y
52  GTO 01
53  CLA
54  END

( 91 bytes )
 

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

-   You key in      CF 02             4    XEQ "IO"    the HP-41  displays     1?
                                                           7   R/S      --------------------     2?
                                                         12   R/S      --------------------     3?
                                                         41   R/S      --------------------     4?
 switch to alpha mode and            "AAA"  R/S      --------------------     5?
                                                   "HP41"  R/S      --------------------     6?
                              simply press    R/S     and   you'll get   4.008     in X-register:  this is the control number  bb.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 every  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!

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

Reference:        "Numerical Analysis" - Francis Scheid - ( McGraw-Hill )  ISBN:  07-055197-9
 
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall