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

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