hp41programs

Complex Complex Functions for the HP-41

Overview

1°) Complex Operations
2°) Complex Equations

   a) Secant Method
   b) Quadratic Interpolation
   c) Newton's Method
   d) A Better Method

1°) Complex Operations

-The global labels "Z+Z" "Z-Z" "Z*Z" "Z/Z" mean addition, subtraction, multiplication and division.
"Z^2" "SQRTZ" "LNZ" "E^Z" "1/Z" mean square, square root, logarithm, exponential, reciprocal.
"SINZ" "COSZ" "TANZ" "ASINZ" "ACOSZ" "ATANZ" correspond to sine, cosine, tangent and their inverses.
"SHZ" "CHZ" "THZ" "ASHZ" "ACHZ" "ATHZ" correspond to hyperbolic sine, hyperbolic cosine, hyperbolic tangent and their inverses:
-I've used the French notations ( sh z instead of sinh z ... etc ... ) which have the advantage of saving several bytes and keystrokes.
"Z^X" and "Z^Z" stand for raising a complex number to a real power and a complex power respectively.

-These programs use no data register and no external subroutine.
-However, ArcCosh z clear synthetic registers M &N - which can be replaced by R00 & R01.
-And, for instance, "TANZ" calls "THZ" as a subroutine since tan z = i.tanh (-i.z) and similarly for many other complex functions.
-In these cases, local labels have been inserted after the global labels ( for example, XEQ 03 LBL 03 instead of XEQ "THZ",
which saves 1 byte and executes faster ) and so on ...
-Lines 144 & 156 are only useful te restore the DEG mode: they may be deleted if RAD is your favorite angular mode.

-To get the magnitude of a complex number, simply press R-P.

Formulae:

    Sinh z = ( e^z - e^-z )/2                           ArcSinh z = Ln [ z + ( z² + 1 )^1/2 ]
    Cosh z = ( e^z + e^-z )/2                          ArcCosh z = Ln [ z + ( z + 1 )^1/2 ( z - 1 )^1/2 ]
    Tanh z = ( e^2z - 1 )/( e^2z + 1 )               ArcTanh z = [ Ln ( 1 + z ) - Ln ( 1 - z ) ]/2

    Sin z = -i Sinh iz                                  ArcSin z = -i ArcSinh iz
    Cos z = Cosh iz                                  ArcCos z = PI/2 - ArcSin z
    Tan z = i Tanh -iz                                ArcTan z = i ArcTanh -iz

Data Registers: /
Flags: /
Subroutines: /

01 LBL "SINZ"
02 X<>Y
03 CHS
04 XEQ 01
05 CHS
06 X<>Y
07 RTN
08 LBL "TANZ"
09 CHS
10 X<>Y
11 XEQ 03
12 X<>Y
13 CHS
14 RTN
15 LBL "ASINZ"
16 LBL 04
17 X<>Y
18 CHS
19 XEQ 05
20 CHS
21 X<>Y
22 RTN
23 LBL "ATANZ"
24 CHS
25 X<>Y
26 XEQ 06
27 X<>Y
28 CHS
29 RTN
30 LBL "SHZ"
31 LBL 01
32 XEQ 00
33 XEQ 01
34 GTO 10
35 LBL "COSZ"
36 X<>Y
37 CHS

38 LBL "CHZ"
39 XEQ 00
40 XEQ 02
41 GTO 10
42 LBL "ATHZ"
43 LBL 06
44 RCL X
45 1
46 ST- Z
47 +
48 R^
49 X<>Y
50 XEQ 12
51 R^
52 CHS
53 R^
54 CHS
55 XEQ 12
56 XEQ 01
57 LBL 10
58 2
59 ST/ Z
60 /
61 RTN
62 LBL "THZ"
63 LBL 03
64 2
65 ST* Z
66 *
67 XEQ 08
68 RCL X
69 1
70 ST- Z
71 +
72 R^
73 X<>Y
74 LBL "Z/Z"

75 R-P
76 R^
77 R^
78 R-P
79 R^
80 ST- Z
81 X<> T
82 /
83 P-R
84 RTN
85 LBL 00
86 RCL Y
87 RCL Y
88 XEQ 08
89 R^
90 CHS
91 R^
92 CHS
93 GTO 08
94 LBL "ACOSZ"
95 XEQ 04
96 X<>Y
97 CHS
98 PI
99 2
100 /
101 R^
102 -
103 RTN
104 LBL "ACHZ"
105 STO M
106 X<>Y
107 STO N
108 X<>Y
109 1
110 +
111 XEQ 07

112 RCL M
113 1
114 -
115 RCL N
116 X<>Y
117 XEQ 07
118 XEQ 03
119 X<>Y
120 0
121 X<> N
122 +
123 X<>Y
124 0
125 X<> M
126 +
127 GTO 12
128 LBL "ASHZ"
129 LBL 05
130 RCL Y
131 RCL Y
132 XEQ 05
133 SIGN
134 ST* X
135 ST+ L
136 X<> L
137 XEQ 07
138 XEQ 02
139 LBL "LNZ"
140 LBL 12
141 RAD
142 R-P
143 LN
144 DEG
145 RTN
146 LBL "Z^Z"
147 R^
148 R^

149 XEQ 12
150 XEQ 03
151 LBL "E^Z"
152 LBL 08
153 RAD
154 E^X
155 P-R
156 DEG
157 RTN
158 LBL "Z^X"
159 RDN
160 R-P
161 R^
162 ST* Z
163 Y^X
164 P-R
165 RTN
166 LBL "Z-Z"
167 LBL 01
168 CHS
169 X<>Y
170 CHS
171 X<>Y
172 LBL "Z+Z"
173 LBL 02
174 X<>Y
175 ST+ T
176 RDN
177 +
178 RTN
179 LBL "Z*Z"
180 LBL 03
181 R-P
182 R^
183 R^
184 R-P
185 X<>Y

186 ST+ T
187 RDN
188 *
189 P-R
190 RTN
191 LBL "1/Z"
192 R-P
193 1/X
194 X<>Y
195 CHS
196 X<>Y
197 P-R
198 RTN
199 LBL "Z^2"
200 LBL 05
201 R-P
202 X^2
203 X<>Y
204 ST+ X
205 X<>Y
206 P-R
207 RTN
208 LBL "SQRTZ"
209 LBL 07
210 R-P
211 SQRT
212 SIGN
213 ST+ X
214 ST/ Y
215 X<> L
216 P-R
217 END

( 448 bytes / SIZE 000 )

1°) First case: Function of one complex number


STACK	INPUTS	OUTPUTS
Y	y	v
X	x	u

where z = x+i.y ; f(z) = u+i.v

Example: z = 2+3.i

-Calculate z² ; z^1/2 ; 1/z ; exp(z) ; Ln(z) ; sin z ; cos z ; tan z ; Asin z ; Acos z ; Atan z ; Sinh z ; Cosh z ; Tanh z ; Asinh z ; Acosh z ; Atanh z

3 ENTER^   2   XEQ "Z^2"        >>>      -5      X<>Y      12        whence    (2+3.i)²   =   -5+12.i
3 ENTER^   2   XEQ "SQRTZ" >>>   1.6741 X<>Y   0.8960   whence    (2+3.i)^1/2 =   1.6741+0.8960.i
3 ENTER^   2   XEQ "1/Z"         >>>   0.1538 X<>Y -0.2308   whence   1/(2+3.i) =    0.1538-0.2308.i
3 ENTER^   2   XEQ "E^Z"        >>> -7.3151 X<>Y   1.0427   whence   exp(2+3.i) = -7.3151+1.0427.i
3 ENTER^   2   XEQ "LNZ"       >>>   1.2825 X<>Y    0.9828   whence   Ln(2+3.i) =   1.2825+0.9828.i
3 ENTER^   2   XEQ "SINZ"     >>>    9.1545 X<>Y   -4.1689 whence   Sin(2+3.i) =   9.1545-4.1689.i
3 ENTER^   2   XEQ "COSZ"    >>>   -4.1896 X<>Y -9.1092 whence   Cos(2+3.i) = -4.1896-9.1092.i
3 ENTER^   2   XEQ "TANZ"    >>>   -0.0038 X<>Y   1.0032 whence    Tan(2+3.i) = -0.0038+1.0032.i
3 ENTER^   2   XEQ "ASINZ"   >>>    0.5707 X<>Y   1.9834 whence   Asin(2+3.i) =   0.5707+1.9834.i
3 ENTER^   2   XEQ "ACOSZ" >>>   1.0001 X<>Y -1.9834 whence   Acos(2+3.i) = 1.0001-1.9834.i
3 ENTER^   2   XEQ "ATANZ" >>>   1.4099 X<>Y   0.2291 whence    Atan(2+3.i) = 1.4099+0.2291.i
3 ENTER^   2   XEQ "SHZ"        >>> -3.5906 X<>Y 0.5309 whence    Sinh(2+3.i) = -3.5906+0.5309.i
3 ENTER^   2   XEQ "CHZ"       >>>   -3.7245 X<>Y 0.5118 whence    Cosh(2+3.i) = -3.7245+0.5118.i
3 ENTER^   2   XEQ "THZ"       >>>    0.9654   X<>Y -0.0099 whence   Tanh(2+3.i) = 0.9654-0.0099.i
3 ENTER^   2   XEQ "ASHZ"    >>>    1.9686   X<>Y   0.9647 whence   Asinh(2+3.i) = 1.9686+0.9647.i
3 ENTER^   2   XEQ "ACHZ"   >>>     1.9834   X<>Y 1.0001 whence    Acosh(2+3.i) = 1.9834+1.0001.i
3 ENTER^   2   XEQ "ATHZ"    >>>    0.1469   X<>Y   1.3390 whence    Atanh(2+3.i) = 0.1469+1.3390.i

Note: "Z^2" "SQRTZ" "1/Z" "E^Z" and "LNZ" save registers Z and T ( this feature is important when they are called as subroutines by the other functions ).

2°) Second case: raising a complex number to a real power


STACK	INPUTS	OUTPUTS
Z	y	/
Y	x	v
X	r	u

where z = x+i.y ; z^r = u+i.v

Example: Evaluate (2+3.i)^1/5

3 ENTER^ 2 ENTER^ 5 1/X XEQ "Z^X" >>> 1.2675 X<>Y 0.2524 whence (2+3.i)^1/5 = 1.2675+0.2524.i

Note: "Z^X" saves T-register in registers Z and T.

3°) Third case: Function of two complex numbers.


STACK	INPUTS	OUTPUTS
T	y	/
Z	x	/
Y	y'	v
X	x'	u

where z = x+i.y ; z' = x'+i.y' ; f(z;z') = u+i.v

Example: z = 2+3.i ; z' = 4+7.i Compute z+z' ; z-z' ; z.z' ; z/z' ; z^z'

3 ENTER^   2 ENTER^   7 ENTER^   4 XEQ "Z+Z"    gives      6       X<>Y       10      whence    z+z' = 6+10.i
3 ENTER^   2 ENTER^   7 ENTER^   4 XEQ "Z-Z"     gives     -2       X<>Y       -4      whence    z-z' = -2-4.i
3 ENTER^   2 ENTER^   7 ENTER^   4 XEQ "Z*Z"    gives   -13       X<>Y       26      whence    z.z' = -13+26.i
3 ENTER^   2 ENTER^   7 ENTER^   4 XEQ "Z/Z"     gives   0.4462 X<>Y -0.0308   whence   z/z' = 0.4462-0.0308.i
3 ENTER^   2 ENTER^   7 ENTER^   4 XEQ "Z^Z"    gives 0.1638   X<>Y   0.0583    whence    z^z' = 0.1638+0.0583.i

Notes:

-This program may of course be simplified if you can use the M-code functions:

1/Z Z*Z Z/Z Z^2 SQRTZ which are listed in "A few M-code Routines for the HP-41"

-Slightly different formulae may also be used to compute ArcCosh ArcTanh and ArcCos:

    Sinh z = ( e^z - e^-z )/2                           ArcSinh z = Ln [ z + ( z² + 1 )^1/2 ]
    Cosh z = ( e^z + e^-z )/2                          ArcCosh z = Ln [ z + ( z² - 1 )^1/2 ]
    Tanh z = ( e^2z - 1 )/( e^2z + 1 )               ArcTanh z = [ Ln (( 1 + z ) / ( 1 - z )) ]/2

    Sin z = -i Sinh iz                                  ArcSin z = -i ArcSinh iz
    Cos z = Cosh iz                                  ArcCos z = -i ArcCosh z
    Tan z = i Tanh -iz                                ArcTan z = i ArcTanh -iz

-These definitions are possible and easier to program since no synthetic register is necessary.
-But they don't always give the principal value of some inverse hyperbolic functions.
-Anyway, here is this simpler program:

Data Registers: /
Flags: /
Subroutines: /

01 LBL "SINZ"
02 X<>Y
03 CHS
04 XEQ 01
05 CHS
06 X<>Y
07 RTN
08 LBL "TANZ"
09 CHS
10 X<>Y
11 XEQ 03
12 X<>Y
13 CHS
14 RTN
15 LBL "ASINZ"
16 X<>Y
17 CHS
18 XEQ 05
19 CHS
20 X<>Y
21 RTN
22 LBL "ATANZ"
23 CHS
24 X<>Y
25 XEQ 06
26 X<>Y
27 CHS
28 RTN
29 LBL "SHZ"
30 LBL 01
31 XEQ 00
32 XEQ 01
33 GTO 10
34 LBL "COSZ"
35 X<>Y
36 CHS
37 LBL "CHZ"
38 XEQ 00
39 XEQ 02

40 GTO 10
41 LBL "ATHZ"
42 LBL 06
43 RCL X
44 1
45 ST+ Z
46 X<>Y
47 -
48 R^
49 CHS
50 X<>Y
51 XEQ 07
52 XEQ 12
53 LBL 10
54 2
55 ST/ Z
56 /
57 RTN
58 LBL "THZ"
59 LBL 03
60 2
61 ST* Z
62 *
63 XEQ 08
64 RCL X
65 1
66 ST- Z
67 +
68 R^
69 X<>Y
70 LBL "Z/Z"
71 LBL 07
72 R-P
73 R^
74 R^
75 R-P
76 R^
77 ST- Z
78 X<> T

79 /
80 P-R
81 RTN
82 LBL 00
83 RCL Y
84 RCL Y
85 XEQ 08
86 R^
87 CHS
88 R^
89 CHS
90 GTO 08
91 LBL "ACOSZ"
92 XEQ 04
93 CHS
94 X<>Y
95 RTN
96 LBL "ACHZ"
97 LBL 04
98 XEQ 13
99 CHS
100 GTO 09
101 LBL "ASHZ"
102 LBL 05
103 XEQ 13
104 LBL 09
105 ST+ L
106 X<> L
107 XEQ 11
108 XEQ 02
109 LBL "LNZ"
110 LBL 12
111 RAD
112 R-P
113 LN
114 DEG
115 RTN
116 LBL 13
117 RCL Y

118 RCL Y
119 XEQ 05
120 SIGN
121 ST* X
122 RTN
123 LBL "Z^Z"
124 R^
125 R^
126 XEQ 12
127 XEQ 03
128 LBL "E^Z"
129 LBL 08
130 RAD
131 E^X
132 P-R
133 DEG
134 RTN
135 LBL "Z^X"
136 RDN
137 R-P
138 R^
139 ST* Z
140 Y^X
141 P-R
142 RTN
143 LBL "Z-Z"
144 LBL 01
145 CHS
146 X<>Y
147 CHS
148 X<>Y
149 LBL "Z+Z"
150 LBL 02
151 X<>Y
152 ST+ T
153 RDN
154 +
155 RTN
156 LBL "Z*Z"

157 LBL 03
158 R-P
159 R^
160 R^
161 R-P
162 X<>Y
163 ST+ T
164 RDN
165 *
166 P-R
167 RTN
168 LBL "1/Z"
169 R-P
170 1/X
171 X<>Y
172 CHS
173 X<>Y
174 P-R
175 RTN
176 LBL "Z^2"
177 LBL 05
178 R-P
179 X^2
180 X<>Y
181 ST+ X
182 X<>Y
183 P-R
184 RTN
185 LBL "SQRTZ"
186 LBL 11
187 R-P
188 SQRT
189 SIGN
190 ST+ X
191 ST/ Y
192 X<> L
193 P-R
194 END

( 415 bytes / SIZE 000 )

-Same instructions and - in most cases - same results...

2°) Complex Equations: f(z) = 0

a) The "Secant" Method

-This program uses an iterative algorithm.
-Starting with 2 distinct initial guesses z₀ = x₀+i.y₀ ; z₁ = x₁+i.y₁ , it calculates z_k+1 = z_k - f(z_k).(z_k - z_k-1) / [ f(z_k) - f(z_k-1) ] ( k = 1 ; 2 ; 3 ; ....... )
-The process is repeated until 2 consecutive approximations are equal.

Data Registers: • R00 = Function name ( Global label, maximum of 6 characters. This register is to be initialized before executing "ROOTZ" )

                                R01 = x   R03 = a   R05 = x'
                                 R02 = y   R04 = b   R06 = y'   where f(z) = f(x+i.y) = a+i.b
Flag: /
Subroutine: A program which computes f(z) = a+i.b assuming x is in X-register and y is in Y-register upon entry
                       In other words, the stack      Y = y           Y = b
                        must be changed from          X = x    to    X = a      where z = x+i.y   and   f(z) = f(x+i.y) = a+i.b

01 LBL "ROOTZ"
02 STO 01
03 X<>Y
04 STO 02
05 R^
06 STO 06
07 R^
08 STO 05
09 XEQ IND 00
10 STO 03
11 X<>Y
12 STO 04

13 LBL 01
14 VIEW 01
15 RCL 02
16 RCL 01
17 XEQ IND 00
18 RCL 04
19 RCL 03
20 R^
21 STO 04
22 ST- Z
23 R^
24 STO 03

25 ST- Z
26 R-P
27 R^
28 R^
29 R-P
30 X<>Y
31 ST- T
32 RDN
33 X#0?
34 /
35 RCL 02
36 ENTER^

37 X<> 06
38 -
39 RCL 01
40 STO L
41 X<> 05
42 ST- L
43 X<> L
44 R-P
45 X<>Y
46 ST+ T
47 RDN
48 *

49 P-R
50 ST+ 01
51 X<>Y
52 ST+ 02
53 E-8
54 LASTX
55 X>Y?
56 GTO 01
57 RCL 02
58 RCL 01
59 END

( 86 bytes / SIZE 007 )


STACK	INPUTS	OUTPUTS
T	y₀	/
Z	x₀	/
Y	y₁	y
X	x₁	x

Example: Find a root of sinh z + z² + Pi = 0. First, we load a routine to compute f(z) , for instance:

LBL "T"
STO 07
X<>Y
STO 08
X<>Y
XEQ "SHZ"
RCL 08
RCL 07
XEQ "Z^2"
XEQ "Z+Z"
PI
+
RTN

And if we choose z₀ = 0 ; z₁ = 1+i

   T ASTO 00
   0 ENTER^   ENTER^
   1 ENTER^   XEQ "ROOTZ" the succesive x-values are displayed

and we get -0.278189856 ( execution time = 77 seconds )
X<>Y 1.812880366

Whence -0.278189856 + 1.812880366.i is a solution. ( There are many other solutions like 4.419361546 + 4.697881722.i ... etc ... )

Note:

-The iterations stop when z_k+1 = z_k or when f(z_k+1) = f(z_k) Therefore, check f(z) in registers R03 & R04

b) Quadratic Interpolation

-Starting with 3 initial guesses: z₀ = x₀+i.y₀ ; z₁ = x₁+i.y₁ ; z₂ = x₂+i.y₂ , this program calculates

   A = ( z₁-z₀ ).f(z₂) + ( z₀-z₂ ).f(z₁) + ( z₂-z₁ ).f(z₀)
   B = ( z₁-z₀ ).( 2.z₂-z₁-z₀ ).f(z₂) - ( z₀-z₂ )².f(z₁) + ( z₂-z₁ )².f(z₀)
   C = ( z₂-z₁ ).( z₁-z₀ ).( z₂-z₀ ).f(z₂)

and z₃ = z₂ - 1 / [ B/(2.C) + E.( B²/(4C²) - A/C )^1/2 ] where E = +1 or -1 ( the sign is choosen to produce the larger denominator )

-The process is repeated until | z_k+1-z_k | is smaller than 10^-8 or C = 0 ( check f(z) in R03 & R04 )

Data Registers: • R00 = Function name ( Global label, maximum of 6 characters. This register is to be initialized before executing "RTZ2" )

R01 = x R03 = a R05 thru R20: temp
R02 = y R04 = b where f(z) = f(x+i.y) = a+i.b
Flag: /
Subroutines: "Z+Z" "Z*Z" "Z^2" "SQRTZ"

and a program which computes f(z) = a+i.b assuming x is in X-register and y is in Y-register upon entry

The stack must Y = y Y = b where z = x+i.y
be changed from X = x to X = a and f(z) = f(x+i.y) = a+i.b

( Registers R03 R04 R13 R14 ... etc ... may be used by this subroutine )

-Line 58 is a three-byte GTO 02
-Line 170 is a three-byte GTO 01

01 LBL "RTZ2"
02 STO 01
03 RDN
04 STO 02
05 X<>Y
06 STO 05
07 STO 06
08 STO 09
09 STO 10
10 ST- T
11 -
12 R^
13 XEQ IND 00
14 STO 07
15 X<>Y
16 STO 08
17 RCL 02
18 RCL 01
19 RCL 05
20 ST+ X
21 ST- Z
22 -
23 XEQ IND 00
24 STO 11
25 X<>Y
26 STO 12
27 LBL 01
28 VIEW 01
29 RCL 02
30 RCL 01
31 XEQ IND 00
32 STO 03
33 X<>Y
34 STO 04
35 X<>Y

36 RCL 10
37 RCL 09
38 XEQ "Z*Z"
39 STO 13
40 X<>Y
41 STO 14
42 RCL 06
43 RCL 10
44 +
45 STO 20
46 RCL 05
47 RCL 09
48 +
49 STO 19
50 RCL 14
51 RCL 13
52 XEQ "Z*Z"
53 RCL 06
54 RCL 05
55 XEQ "Z*Z"
56 R-P
57 X=0?
58 GTO 02
59 CHS
60 STO 15
61 X<>Y
62 STO 16
63 RCL 06
64 RCL 20
65 +
66 RCL 05
67 RCL 19
68 +
69 RCL 14
70 RCL 13

71 XEQ "Z*Z"
72 STO 17
73 X<>Y
74 STO 18
75 RCL 20
76 RCL 19
77 RCL 08
78 RCL 07
79 XEQ "Z*Z"
80 ST- 13
81 X<>Y
82 ST- 14
83 X<>Y
84 RCL 20
85 RCL 19
86 XEQ "Z*Z"
87 ST- 17
88 X<>Y
89 ST- 18
90 RCL 06
91 RCL 05
92 RCL 12
93 RCL 11
94 XEQ "Z*Z"
95 ST+ 13
96 X<>Y
97 ST+ 14
98 X<>Y
99 RCL 06
100 RCL 05
101 XEQ "Z*Z"
102 ST+ 17
103 X<>Y
104 ST+ 18
105 RCL 14

106 RCL 13
107 R-P
108 STO 13
109 X<>Y
110 STO 14
111 RCL 18
112 RCL 17
113 R-P
114 RCL 16
115 ST- 14
116 ST- Z
117 X<> 15
118 ST/ 13
119 ST+ X
120 /
121 P-R
122 STO 15
123 X<>Y
124 STO 16
125 X<>Y
126 XEQ "Z^2"
127 RCL 14
128 RCL 13
129 P-R
130 XEQ "Z+Z"
131 XEQ "SQRTZ"
132 RCL 16
133 RCL Z
134 ST- 16
135 +
136 RCL 15
137 RCL Z
138 ST- 15
139 +
140 R-P

141 X<>Y
142 RCL 16
143 RCL 15
144 R-P
145 R^
146 X<Y?
147 RCL Y
148 X#0?
149 1/X
150 R^
151 CHS
152 X<>Y
153 P-R
154 ST+ 01
155 X<> 05
156 STO 09
157 X<>Y
158 ST+ 02
159 X<> 06
160 STO 10
161 RCL 03
162 X<> 07
163 STO 11
164 RCL 04
165 X<> 08
166 STO 12
167 E-8
168 LASTX
169 X>Y?
170 GTO 01
171 LBL 02
172 RCL 02
173 RCL 01
174 END

( 281 bytes / SIZE 021 )


STACK	INPUTS	OUTPUTS
Z	h	/
Y	y₀	y
X	x₀	x

( The 3 initial guesses are actually z₀ = x₀+i.y₀ ; z₀-(1+i).h ; z₀-2.(1+i).h )

Example: With the same equation sinh z + z² + Pi = 0. We load a routine to compute f(z)

LBL "T"
STO 03
X<>Y
STO 04
X<>Y
XEQ "SHZ"
RCL 04
RCL 03
XEQ "Z^2"
XEQ "Z+Z"
PI
+
RTN

And if we choose z₀ = 1+i and h = 1

   T ASTO 00
   1 ENTER^   ENTER^
       XEQ "RTZ2" the succesive x-values are displayed and we get   -0.278189857                  ( execution time = 156 seconds )
                                                                                          X<>Y      1.812880366

Whence z = -0.278189857 + 1.812880366.i

Note: In this example, "ROOTZ" is faster than "RTZ2" but if f is a very complicated function, the advantage is clearly with "RTZ2"

c) Newton's Method

-Only one initial guess z₀ = x₀+i.y₀ is needed but we have to compute the first derivative f '(z)

Formula: z_k+1 = z_k - f(z_k)/f '(z_k)

Data Registers: • R00 = f name • R03 = f ' name ( these 2 registers are to be initialized before executing "NWTZ" )

R01 = x
R02 = y R04 = | f '(z) | R05: temp
Flag: /
Subroutines:

-A program which computes f(z) assuming x is in X-register and y is in Y-register upon entry
-a program which computes f '(z) assuming x is in X-register and y is in Y-register upon entry

01 LBL "NWTZ"
02 STO 01
03 X<>Y
04 STO 02
05 LBL 01
06 VIEW 01
07 RCL 02
08 RCL 01

09 XEQ IND 03
10 R-P
11 STO 04
12 X<>Y
13 STO 05
14 RCL 02
15 RCL 01
16 XEQ IND 00

17 R-P
18 ENTER^
19 X<> 04
20 X#0?
21 /
22 X<>Y
23 RCL 05
24 -

25 X<>Y
26 P-R
27 ST- 01
28 X<>Y
29 ST- 02
30 E-8
31 LASTX
32 X>Y?

33 GTO 01
34 RCL 04
35 RCL 02
36 RCL 01
37 END

( 55 bytes / SIZE 006 )


STACK	INPUTS	OUTPUTS
Z	/	\| f(z) \|
Y	y₀	y
X	x₀	x

-The output in Z-register should be a "small" number if z = x+i.y is a "true" root

Example: sinh z + z² + Pi = 0. We have to load 2 routines to compute f(z) and f ' (z) = cosh z + 2.z for instance:

LBL "T"
STO 07
X<>Y
STO 08
X<>Y
XEQ "SHZ"
RCL 08
RCL 07
XEQ "Z^2"
XEQ "Z+Z"
PI
+
RTN
LBL "U"
STO 07
X<>Y
STO 08
X<>Y
XEQ "CHZ"
RCL 08
ST+ X
RCL 07
ST+ X
XEQ "Z+Z"
RTN

Then    T ASTO 00   U ASTO 03   and with z₀ = 1+i
           1 ENTER^    XEQ "NWTZ" the successive x-approximations are displayed
                                      and we get -0.278189857
                                              RDN    1.812880366     whence z = -0.278189857 + 1.812880366.i      ( execution time = 54 seconds )

Note: | f(z) | in Z-register actually concerns the next to last approximation.

d) A Better method

-Starting with one initial guess z₀ = x₀+i.y₀ , we now use the first and second derivatives f '(z) and f "(z) , the recursion formula is:
z_k+1 = z_k - 1 / [ f '(z_k)/f(z_k) + E.( ( f(z_k)/(2.f '(z_k)) )² - f "(z_k) /(2.f (z_k)) )^1/2 ] where E = +1 or -1 ( the sign is choosen to produce the larger denominator )

Data Registers: • R00 = f name • R03 = f ' name • R04 = f " name ( these 3 registers are to be initialized before executing "RTZ4" )

R01 = x
R02 = y R05 = | f '(z) | R06 thru R09: temp
Flag: /
Subroutines:

-A program which computes f(z)    assuming x is in X-register and y is in Y-register upon entry
-A program which computes f '(z)   assuming x is in X-register and y is in Y-register upon entry
-A program which computes f "(z)   assuming x is in X-register and y is in Y-register upon entry

01 LBL "RTZ4"
02 STO 01
03 X<>Y
04 STO 02
05 LBL 01
06 VIEW 01
07 RCL 02
08 RCL 01
09 XEQ IND 03
10 R-P
11 STO 06
12 X<>Y
13 STO 07
14 RCL 02
15 RCL 01
16 XEQ IND 04
17 R-P
18 STO 08
19 X<>Y

20 STO 09
21 RCL 02
22 RCL 01
23 XEQ IND 00
24 R-P
25 STO 05
26 X=0?
27 GTO 02
28 ST+ X
29 ST/ 06
30 ST/ 08
31 X<>Y
32 ST- 07
33 ST- 09
34 RCL 07
35 ST+ X
36 RCL 06
37 X^2
38 P-R

39 RCL 09
40 RCL 08
41 P-R
42 ST- Z
43 RDN
44 ST- Z
45 RDN
46 R-P
47 SQRT
48 X<>Y
49 2
50 /
51 X<>Y
52 P-R
53 RCL 07
54 RCL 06
55 P-R
56 STO 06
57 RDN

58 STO 07
59 RCL Z
60 ST- 07
61 +
62 RCL 06
63 RCL Z
64 ST- 06
65 +
66 R-P
67 X<>Y
68 RCL 07
69 RCL 06
70 R-P
71 R^
72 X<Y?
73 RCL Y
74 X#0?
75 1/X
76 R^

77 CHS
78 X<>Y
79 P-R
80 ST- 01
81 X<>Y
82 ST- 02
83 E-8
84 LASTX
85 X>Y?
86 GTO 01
87 LBL 02
88 RCL 05
89 RCL 02
90 RCL 01
91 END

( 123 bytes / SIZE 010 )


STACK	INPUTS	OUTPUTS
Z	/	\| f(z) \|
Y	y₀	y
X	x₀	x

-The output in Z-register should be a "small" number if z = x+i.y is a "true" root

Example: sinh z + z² + Pi = 0. We have to load 3 routines to compute f(z) ; f ' (z) = cosh z + 2.z and f "(z) = sinh z + 2 , for instance:

   LBL "T"
   STO 10
   X<>Y
   STO 11
   X<>Y
   XEQ "SHZ"
   RCL 11
   RCL 10
   XEQ "Z^2"
   XEQ "Z+Z"
   PI
   +
   RTN
   LBL "U"
   STO 10
   X<>Y
   STO 11
   X<>Y
   XEQ "CHZ"
   RCL 11
   ST+ X
   RCL 10
   ST+ X
   XEQ "Z+Z"
   RTN
   LBL "V"
   XEQ "SHZ"
   2
   +
   RTN

Then    T ASTO 00   U ASTO 03   V ASTO 04 and with z₀ = 1+i
           1 ENTER^    XEQ "RTZ4"    the successive x-approximations are displayed
                                      and we get -0.278189857
                                              RDN    1.812880366     whence z = -0.278189857 + 1.812880366.i      ( execution time = 70 seconds )

Notes:

| f(z) | in Z-register actually concerns the next to last approximation.
Among these 4 root-finding programs, "RTZ4" has the best rate of convergence.