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

-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
HTML