# hp41programs

approx A Successive Approximation Method for the HP-41

Overview

1°) Real Unknowns
2°) Complex Unknowns

-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 systems 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
-To avoid a possible infinite loop, line 21 might be replaced by   E-8     X<Y?   or   E-8   RCL 00   *   X<Y?

 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? 22  GTO 01        23  CLA 24  END

( 43 bytes / SIZE 2n+1 )

 STACK INPUTS OUTPUTS X / 0

Example1:   Solve the system:

x2 / y + y / z - z2 = 0
x.y.z  -  x  -  y2  = 0
x.ln y - z = 0

-First, let's rewrite this system into the form:

x = z / ln y
y = ( x.y.z - x )1/2
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 )

-For instance, if the system:

10.x +   y    -   z   = 1                                 x = ( 1 - y + z ) / 10
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:

-Occasionally, it will also work for a system like:

2.x + 3.y - 4.z = -2                                  x = ( ( 7.x + 3.x.y - 2.x.z ) / 4 )1/2
3.x + 2.y + 5.z = 30      If it's rewritten     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

-To avoid a possible infinite loop, line 33 might be replaced by   E-8     X<Y?   or   E-8   RCL 00   *   X<Y?

 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? 34  GTO 01      35  CLA 36  END

( 59 bytes / SIZE 3n+1 )

 STACK INPUTS OUTPUTS X / 0

Example:    Find a solution of the system:     z = ( z2 - z' )1/3  ,   z' = ( (z')2 - z )1/4

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 )