The Museum of HP Calculators
The Museum of HP Calculators
This program is Copyright © 2003-2006 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.
1°) Real Unknowns
2°) Complex Unknowns (
new )
-If F is a contraction mapping on a complete metric space, then the
equation F ( X ) = X has a unique solution
which is the limit of the sequence: X(k+1) =
F(X(k)) where X(0) is arbitrary.
-This theorem is used in the following program to solve a system of
n equations in n unknowns written in the form:
x1 = f1 ( x1
, ... , xn )
x2 = f2 ( x1
, ... , xn )
.............................
xn = fn ( x1 , ... , xn )
-The advantage of this method is its simplicity: there is no need to
solve a n x n linear system like with "NLS" ( cf "linear and non-linear
system for the HP-41" ).
-It could be used to solve large linear or non-linear systems, with
n > 16 , which would be impossible otherwise with an HP-41.
-Unfortunately, it's not so easy to find the good function F that leads
to convergence...
1°) Real Unknowns
Data Registers: • R00 = n = the number of unknowns ( All these registers are to be initialized before executing "FXN" )
• R01 = x1
• Rn+1 = f1 name
• R02 = x2
• Rn+2 = f2 name
( x1 , ... , xn ) is an initial guess that
you have to choose.
................................................
• Rnn = xn • R2n = fn name
Flags: /
Subroutines: The n functions
fi computing fi ( x1
, ... , xn ) in X-register assuming x1
in R01 ; ...... ; xn in Rnn
-Actually, this program calculates ( x'1 , ... , x'n ) from ( x1 , ... , xn ) by the formulae:
x'1 = f1 ( x1
, x'2 , ... , x'n-1 , x'n
)
x'2 = f2 ( x1
, x2 , ... , x'n-1 , x'n )
.............................
x'n-2 = fn-1 ( x1
, x2 , ... , x'n-1 , x'n )
x'n-1 = fn-1 ( x1
, x2 , ... , xn-1 , x'n )
x'n = fn
( x1 , x2 , ... , xn-1 , xn
)
-In other words, every new estimation of xi
replaces the previous one as soon as it is computed.
-Thus, the program is shorter, it requires less registers and convergence
is improved!
Note: As usual, synthetic registers M and N may be replaced
by any unused data registers, but if , for instance, you replace register
N by R99
replace
line 04 by the two lines: CLX STO 99
01 LBL "FXN"
02 LBL 01
03 VIEW 01
04 CLA
05 RCL 00
06 STO M
07 LBL 02
08 RCL 00
09 RCL M
10 +
11 RCL IND X
12 XEQ IND X
13 ENTER^
14 X<> IND M
15 -
16 ABS
17 ST+ N
18 DSE M
19 GTO 02
20 X<> N
21 X#0?
to avoid a possible infinite loop, line 21 might be replaced by
E-8 or
E-8
22 GTO 01
X<Y?
RCL 00
23 CLA
*
24 END
X<Y?
( 43 bytes / SIZE 2n+1 )
x2 / y + y / z - z2 = 0
x = z / ln y
Example1: Solve the system:
x.y.z - x - y2 = 0
First, let's rewrite this system into the form:
y = ( x.y.z - x )1/2
x.ln y - z = 0
z = ( x2 / y + y / z )1/2
and let's call these 3 functions: "X" "Y" "Z"
LBL "X"
RCL 03
RCL 02
LN
/
RTN
LBL "Y"
RCL 02
RCL 03
*
RCL 01
ST* Y
-
SQRT
RTN
LBL "Z"
RCL 01
X^2
RCL 02
ST/ Y
RCL 03
/
+
SQRT
RTN
Then,
3 STO 00 ( there are 3 functions )
and if we choose ( 2 ; 2 ; 2 ) as initial guess,
2 STO 01 STO 02 STO 03
"X" ASTO 04
"Y" ASTO 05
"Z" ASTO 06
XEQ "FXN" displays the successive x-values and finally:
x = 1.894868809 in R01
y = 2.457696043 in R02
z = 1.703912160 in R03
Example2: This program may also be used
to solve n x n linear systems A.X = B , especially when A is a "dominant
diagonal" matrix.
( If n > 16 , it couldn't be solved by a classic linear-system program
because at most 319 registers are available )
10.x + y - z = 1
x = ( 1 - y + z ) / 10
For instance, if the system 2.x + 11.y + 3.z
= 4 is rewritten:
y = ( 4 - 2.x - 3.z ) / 11 the
guess ( 0 ; 0 ; 0 ) assures convergence.
3.x - y - 12. z = 2
z = ( -2 + 3.x - y ) / 12
( The solution is x = 0.040098200
; y = 0.408346972 ; z = -0.190671031
)
Remark:
2.x + 3.y - 4.z = -2
x = ( ( 7.x + 3.x.y - 2.x.z ) / 4 )1/2
-Occasionally, it will also work for a system like:
3.x + 2.y + 5.z = 30 If it's reshaped in
the form y = ( ( -2.y - 2.x.y + 4.y.z ) / 3 )1/2
4.x - 3.y + 2.z = 7
z = ( ( 30.z - 3x.z - 2.y.z ) / 5 )1/2
The initial guess ( 3 ; 3 ; 3 ) leads to the solution!
Exercise: Find a solution of the following system:
x = ( y2 + z2 - z.t + u.v + w2
)1/2 ; y = ( x + y + z + t - u
- v - w ) / ( 2.y ) ; z = ( x.y.(
u + v ) + z.( t - w ) )1/3
t = ( ( x + y + z ).( u + v + w ) )1/2 ;
u = 2.( z + 2.t ) / ( x2 + y2 + u2 + v2
+ w2 ) ; v = ( x.y2.z + t.u.w2
)1/4 ; w = ( x.y.z + 2.t.u.v )1/3
Answer: x = 2.820317813 ; y = 1.992285489
; z = 3.435794297 ; t = 8.547543120 ; u = 0.924440439
; v = 3.672472185 ; w = 4.260625159
2°) Complex Unknowns
z = x + i.y
Data Registers: • R00 = n = the number of unknowns ( All these registers are to be initialized before executing "FZN" )
• R01 = x1 •
R02 = y1 • R2n+1 =
f1 name
• R03 = x2 •
R04 = y2 • R2n+2 =
f2 name
( z1 = x1 + i.y1 , ... , zn
= xn + i.yn ) is an initial guess you
have to choose.
....................................................................
• R2n-1 = xn • R2n = yn • R3n = fn name
Flags: /
Subroutines: The n functions
fi computing fi ( z1
, ... , zn ) = X + i.Y in X- and Y-registers assuming
z1 in R01 R02 , ...... , zn in R2n-1
R2n
01 LBL "FZN"
02 LBL 01
03 VIEW 01
04 CLA
05 RCL 00
06 ST+ X
07 STO M
08 LBL 02
09 RCL 00
10 ST+ X
11 RCL M
12 2
13 /
14 +
15 RCL IND X
16 XEQ IND X
17 X<>Y
18 ENTER^
19 X<> IND M
20 -
21 ABS
22 X<>Y
23 ENTER^
24 DSE M
25 X<>IND M
26 -
27 ABS
28 +
29 ST+ N
30 DSE M
31 GTO 02
32 X<> N
33 X#0?
to avoid a possible infinite loop, line 33 may be replaced by
E-8 or
E-8
34 GTO 01
X<Y?
RCL 00
35 CLA
*
36 END
X<Y?
( 59 bytes / SIZE 3n+1 )
Example: Find a solution of the system: z = ( z2 - z' )1/3 , z' = ( (z')2 - z )1/4
-Load the following routine
LBL "Z"
RCL 02
RCL 01
XEQ "Z^2"
( cf "Complex Functions for the HP-41" )
RCL 03
-
X<>Y
RCL 04
-
X<>Y
3
1/X
XEQ "Z^X"
( cf "Complex Functions for the HP-41" )
RTN
LBL "T"
RCL 04
RCL 03
XEQ "Z^2"
( cf "Complex Functions for the HP-41" )
RCL 01
-
X<>Y
RCL 02
-
X<>Y
4
1/X
XEQ "Z^X"
( cf "Complex Functions for the HP-41" )
RTN
Z ASTO 05 T ASTO 06
2 STO 00 and if we choose z = z' = 1 + i
as initial guesses 1 STO 01 STO 02 STO 03
STO 04
XEQ "FZN" the successive x1 are displayed and finally:
z = R01 + i R02 = 1.038322757 + 0.715596476
i
z' = R03 + i R04 = 1.041713085 - 0.462002405
i
( execution time = 5mn43s )
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall