The Museum of HP Calculators
The Museum of HP Calculators
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
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