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⁽⁰⁾ is arbitrary.
-This theorem is used in the following program to solve a system of n equations in n unknowns written in the form:

      x₁ = f₁ ( x₁ , ... , x_n )
      x₂ = f₂ ( x₁ , ... , x_n )
       .............................

x_n = f_n ( x₁ , ... , x_n )

-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 = x₁        • Rn+1 = f₁ name
                                  •   R02 = x₂        • Rn+2 = f₂ name          ( x₁ , ... , x_n ) is an initial guess that you have to choose.
                                       ................................................

• Rnn = x_n • R2n = f_n name

Flags: /
Subroutines: The n functions f_i computing f_i ( x₁ , ... , x_n ) in X-register assuming x₁ in R01 ; ...... ; x_n in Rnn

-Actually, this program calculates ( x'₁ , ... , x'_n ) from ( x₁ , ... , x_n ) by the formulae:

      x'₁ = f₁ ( x₁ , x'₂ , ... , x'_n-1 , x'_n )
      x'₂ = f₂ ( x₁ , x₂ , ... , x'_n-1 , x'_n )
       .............................

    x'_n-2 = f_n-1 ( x₁ , x₂ , ... , x'_n-1 , x'_n )
    x'_n-1 = f_n-1 ( x₁ , x₂ , ... , x_n-1 , x'_n )
      x'_n = f_n   ( x₁ , x₂ , ... , x_n-1 , x_n )

-In other words, every new estimation of x_i 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:

      x² / y + y / z - z² = 0
      x.y.z - x - y² = 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 = ( x² / 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

Note:

-In fact, we can write a simple program for each problem without any subroutine:

01 LBL "FX"
02 LBL 01
03 VIEW 01
04 RCL 03
05 RCL 02
06 LN
07 /
08 ENTER

09 X<> 01
10 -
11 ABS
12 STO 00
13 RCL 02
14 RCL 03
15 *
16 RCL 01

17 ST* Y
18 -
19 SQRT
20 ENTER
21 X<> 02
22 -
23 ABS
24 ST+ 00

25 RCL 01
26 X^2
27 RCL 02
28 ST/ Y
29 RCL 03
30 /
31 +
32 SQRT

33 ENTER
34 X<> 03
35 -
36 ABS
37 RCL 00
38 +
39 X#0?
40 GTO 01
41 END

( 56 bytes / SIZE nnn+1 )

STACK	INPUTS	OUTPUTS
X	/	0

-Let's store the initial guesses into R01 to Rnn ( 2 STO 01 STO 02 STO 03 ) and XEQ "FX"

-The successive x-values are displayed 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 = ( y² + z² - z.t + u.v + w² )^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 ) / ( x² + y² + u² + v² + w² ) ; v = ( x.y².z + t.u.w² )^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 = x₁      •   R02 = y₁       • R2n+1 = f₁ name
                                  •   R03 = x₂      •   R04 = y₂       • R2n+2 = f₂ name          ( z₁ = x₁ + i.y₁ , ... , z_n = x_n + i.y_n ) is an initial guess you have to choose.
                                       ....................................................................

• R2n-1 = x_n • R2n = y_n • R3n = f_n name

Flags: /
Subroutines: The n functions f_i computing f_i ( z₁ , ... , z_n ) = X + i.Y in X- and Y-registers assuming z₁ in R01 R02 , ...... , z_n 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 = ( z² - z' )^1/3 , z' = ( (z')² - 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 x₁ 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 )