The Museum of HP Calculators

# A Successive Approximation Method for the HP-41

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.

## Overview

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

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 )