ZDERIVE41

Complex Derivatives & Riemannian Manifolds for the HP-41

Overview

0°) An M-Code Routine
1°) Complex Operations & Elementary Functions
2°) First Derivatives of a Function of N variables ( N < 5 )
3°) Second Derivatives of a Function of N variables ( N < 5 )
4°) Third Derivatives of a Function of N variables ( N < 5 )
5°) Hessian Matrix
6°) Laplacian of a Function of N variables ( N < 5 )
7°) Biharmonic Operator
8°) Curl, Divergence, Gradient & Laplacian ( 2 variants )
9°) Gaussian Curvature & Mean Curvature of a Surface
10°) Curvature & Torsion of a Parametric 3D-Curve
11°) Curvature & Torsion of a Parametric 4D-Curve
12°) Complex Linear Systems
13°) Complex Riemannian Manifolds - Initialization - Metric Tensor & Christoffel Symbols
14°) Recalling Gij & Cij^k
15°) Curvature Tensor ( 0-4 )
16°) Curvature Tensor ( 1-3 )
17°) Ricci Tensor
18°) Scalar Riemannian Curvature
19°) Einstein Tensor
20°) Weyl Tensor
21°) Covariant Derivative of a Vector Field
22°) Curl of a Covariant Vector Field
23°) Divergence of a Contravariant Vector Field
24°) Laplacian & Gradient of a Scalar Field
25°) Covariant Vector <> Contravariant Vector
26°) Eigenvalues & Eigenvector of a Complex Matrix
27°) Complex Polar <> Rectangular Conversion
28°) Complex Spherical <> Rectangular Conversion
29°) Dot Product & Hermitian Product
30°) E^Z-1
31°) Complex Einstein's Functions
32°) Complex Planck Radiation Function

-Most of these programs are in the module: ZDERIVE41 that may be found here.

Notes:

-Most of the formulas to estimate derivatives are of order 10.

-With formulas of order 10, the execution time may be prohibitive to compute all tensors..
-So, you can also choose 4th-order methods with"ZFD" and "ZSD", to get faster - but less accurate - results.

-However, a good emulator like V41 or a 41CL is recommended.

Notations:

-In the comments below, I've used Einstein summation convention:

-When an index appears as an upper index & a lower index in a monomial and is not otherwise defined,
it means summation of that term over all the values of the index.

-For example, in a 3 dimensional space,

a_k b^k = a₁ b¹ + a₂ b² + a₃ b³

         g_ij aⁱ a^j = g₁₁ a¹ a¹ + g₁₂ a¹ a² + g₁₃ a¹ a³
                    + g₂₁ a² a¹ + g₂₂ a² a² + g₂₃ a² a³
                    + g₃₁ a³ a¹ + g₃₂ a³ a² + g₃₃ a³ a³

-If the tensor g_ij is symmetric ( g_ij = g_ji ) the last expression may be simplified and becomes:

g_ij aⁱ a^j = g₁₁ a¹ a¹ + g₂₂ a² a² + g₃₃ a³ a³ + 2 ( g₁₂ a¹ a² + g₁₃ a¹ a³ + g₂₃ a² a³ )

-The partial derivative ¶_m g_ij means ¶g_ij / ¶x^m
-Likewise, ¶_km g_ij = ¶g_ij / ¶x^k¶x^m

0°) GIJ name in Alpha register

-In the ZDERIVE41 module, an M-code routine is hidden in the header -ZDERIVE41, it compares the contents of X- and Y-registers i & j

-The function name "GIJ" is placed in the alpha register where I = Inf ( i , j ) & J = Sup ( i , j )

-So, if X-register = i and Y-register = j , the program lines ( with FIX 0 CF 29 ):

X>Y?
X<>Y
"G"
ARCL X
ARCL Y

are replaced by

-ZDERIVE41

0B1 "1"
034 "4"
005 "E"
016 "V"
009 "I"
012 "R"
005 "E"
004 "D"
01A "Z"
02D "-"
04E C=0
168 M=C
1A8 N=C
1E8 O=C
228 P=C
09C PT=5
110 LD@PT- 4
1D0 LD@PT- 7
01C PT=3
0D0 LD@PT- 3
010 LD@PT- 0
0D0 LD@PT- 3
0EE B<>C ALL
0B8 C=Y
10E A=C ALL
0F8 C=X
30E ?A<C ALL
017 JC+02
0AE A<>C ALL
06E A<>B ALL
37C RCR 12
39C PT=0
102 A=C @PT
0EE B<>C ALL
21C PT=2
0FC RCR 10
102 A=C @PT
0AE A<>C ALL
168 M=C
3E0 RTN

1°) Complex Operations & Elementary Functions

Z+Z Z-Z Z*Z Z/Z 1/Z Z^2 SQRTZ ( cf "A Few M-Code Routines for the HP-41" )
and E^Z LNZ are M-Code routines that work like built-in functions except that there is no check for alpha data or overflow or underflow.

E^Z & LNZ calculate exponential & logarithm without changing the angular mode.

09A "Z"
00E "N"
00C "L"
2A0 SETDEC
130 LDI S&X
090 090h
358 ST=C XP
144 CLRF 6
284 CLRF 7
0B8 C=Y
10E A=C ALL
0F8 C=X
128 L=C
070 N=C ALL
125 ?NCXQ
074 1D49
0F0 C<>N ALL
0A8 Y=C
0B0 C=N ALL
3C4 ST=0
115 C=
06C Ln C
0E8 X=C
3E0 RTN

09A "Z"
01E "^"
005 "E"
2A0 SETDEC
3C4 ST=0
0F8 C=X
128 L=C
029 C=
068 E^C
070 N=C ALL
130 LDI S&X
090 090h
358 ST=C XP
144 CLRF 6
284 CLRF 7
0B8 C=Y
10E A=C ALL
0B0 C=N ALL
1D5 ?NCXQ
078 1E75
0E8 X=C
0F0 C<>N ALL
0A8 Y=C
3E0 RTN

-The other functions are focal programs

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

01 LBL "SINZ"
02 X<>Y
03 CHS
04 XROM "SHZ"
05 CHS
06 X<>Y
07 RTN
08 LBL "ACOSZ"
09 CF 07
10 GTO 00
11 LBL "ASINZ"
12 SF 07
13 LBL 00
14 X<>Y
15 CHS
16 XROM "ASHZ"
17 CHS
18 X<>Y
19 FS?C 07
20 RTN
21 X<>Y
22 CHS
23 PI
24 2
25 /
26 R^
27 -

28 RTN
29 LBL "SHZ"
30 SF 07
31 GTO 00
32 LBL "COSZ"
33 X<>Y
34 CHS
35 LBL "CHZ"
36 CF 07
37 LBL 00
38 RCL Y
39 RCL Y
40 E^Z
41 R^
42 CHS
43 R^
44 CHS
45 E^Z
46 FS? 07
47 Z-Z
48 FC?C 07
49 Z+Z
50 2
51 ST/ Z
52 /
53 RTN
54 LBL "ATANZ"

55 CHS
56 X<>Y
57 CF 07
58 GTO 00
59 LBL "ATHZ"
60 SF 07
61 LBL 00
62 RCL X
63 1
64 ST- Z
65 +
66 R^
67 X<>Y
68 LNZ
69 R^
70 CHS
71 R^
72 CHS
73 LNZ
74 Z+Z
75 2
76 ST/ Z
77 /
78 FS?C 07
79 RTN
80 X<>Y
81 CHS

82 RTN
83 LBL "TANZ"
84 CHS
85 X<>Y
86 CF 07
87 GTO 00
88 LBL "THZ"
89 SF 07
90 LBL 00
91 2
92 ST* Z
93 *
94 E^Z
95 RCL X
96 1
97 ST- Z
98 +
99 R^
100 X<>Y
101 Z/Z
102 FS?C 07
103 RTN
104 X<>Y
105 CHS
106 RTN
107 LBL "ASHZ"
108 RCL Y

109 RCL Y
110 Z^2
111 SIGN
112 ST* X
113 ST+ L
114 X<> L
115 SQRTZ
116 Z+Z
117 LNZ
118 RTN
119 LBL "ACHZ"
120 ENTER
121 ENTER
122 STO Q
123 CLX
124 SIGN
125 ST- Z
126 +
127 R^
128 X<>Y
129 SQRTZ
130 R^
131 STO P
132 R^
133 SQRTZ
134 Z*Z
135 ENTER

136 CLX
137 X<> P
138 RCL Q
139 Z+Z
140 LNZ
141 RTN
142 LBL "Z^Z"
143 R^
144 R^
145 E^Z
146 Z*Z
147 E^Z
148 RTN
149 LBL "Z^X"
150 RDN
151 R-P
152 R^
153 ST* Z
154 X=0?
155 X#Y?
156 GTO 00
157 SIGN
158 X<>Y
159 LBL 00
160 Y^X
161 P-R
162 END

( 321 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

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	y
Z	x	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:

"Z^Z" does not save registers Z & T
"ACHZ" uses synthetic registers P & Q

2°) First Derivatives of a Function of N variables ( N < 5 )

-If h > 0 "ZFD" employs a formula of order 10

df/dx = (1/2520.h).[ 2100.( f₁ - f_-1 ) - 600.( f₂ - f_-2 ) + 150.( f₃ - f_-3 ) - 25.( f₄ - f_-4 ) + 2.( f₅ - f_-5 ) ] with f_k = f(x+k.h)

-If h < 0 "ZFD" employs a formula of order 4 ( faster ) i-e

df/dx = (1/12h).[ f(x-2h) - 8.f(x-h) + 8.f(x+h) - f(x+2h) ]

>>> In both cases if the function returns an alpha data, 0 is returned without calculating the whole formula.

Data Registers: • R00 = Function name ( Registers R00 thru R2n are to be initialized before executing "ZFD" )

• R01-R02 = x₁ , • R03-R04 = x₂ , .......... , • R2n-1-R2n = x_n n < 5

R10 = h
R11-R12 = ¶f / ¶x_i R13 ........ R17: temp

Flags: /
Subroutine: A program which computes f(x₁, ..... ,x_n) assuming x₁ , ......... , x_n are in R01 thru R2n

-The following listing combines "ZFD" "ZSD" and "ZLAPN"

01 LBL "ZFD"
02 LBL 01
03 ST+ X
04 DSE X
05 STO 14
06 X<>Y
07 STO 10
08 STO 15
09 RCL IND Y
10 STO 13
11 RCL 15
12 -
13 STO IND 14
14 XEQ IND 00
15 SIGN
16 X#0?
17 GTO 00
18 RCL 13
19 STO IND 14
20 CLX
21 RTN
22 LBL 00
23 X<> L
24 XEQ 06
25 RCL 10
26 X>0?
27 GTO 00
28 CLX
29 8
30 ST* Z
31 *
32 STO 11
33 X<>Y
34 STO 12
35 XEQ 03
36 ST- 11
37 X<>Y
38 ST- 12
39 12
40 LBL 13
41 RCL 10
42 *

43 ST/ 11
44 ST/ 12
45 RCL 13
46 STO IND 14
47 RCL 12
48 RCL 11
49 RTN
50 GTO 01
51 LBL 00
52 CLX
53 21
54 ST* Z
55 *
56 STO 11
57 X<>Y
58 STO 12
59 XEQ 03
60 6
61 ST* Z
62 *
63 ST- 11
64 X<>Y
65 ST- 12
66 XEQ 03
67 1.5
68 ST* Z
69 *
70 ST+ 11
71 X<>Y
72 ST+ 12
73 XEQ 03
74 4
75 ST/ Z
76 /
77 ST- 11
78 X<>Y
79 ST- 12
80 XEQ 03
81 50
82 ST/ Z
83 /
84 ST+ 11

85 X<>Y
86 ST+ 12
87 25.2
88 GTO 13
89 LBL 03
90 RCL 10
91 ST+ 15
92 RCL 13
93 RCL 15
94 -
95 STO IND 14
96 XEQ IND 00
97 LBL 06
98 STO 16
99 X<>Y
100 STO 17
101 RCL 13
102 RCL 15
103 +
104 STO IND 14
105 XEQ IND 00
106 X<>Y
107 RCL 17
108 -
109 X<>Y
110 RCL 16
111 -
112 RTN
113 LBL "ZSD"
114 LBL 02
115 CF 10
116 X=Y?
117 SF 10
118 ST+ X
119 DSE X
120 STO 16
121 RDN
122 ST+ X
123 DSE X
124 STO 15
125 X<>Y
126 STO 10

127 RCL IND 15
128 STO 13
129 RCL IND 16
130 STO 14
131 XEQ IND 00
132 SIGN
133 ENTER
134 X=0?
135 RTN
136 LASTX
137 ST+ X
138 STO 17
139 R^
140 ST+ X
141 STO 18
142 CLX
143 STO 21
144 RCL 10
145 X>0?
146 GTO 00
147 XEQ 03
148 16
149 ST* Z
150 *
151 STO 11
152 X<>Y
153 STO 12
154 XEQ 03
155 ST- 11
156 X<>Y
157 ST- 12
158 12
159 LBL 14
160 FC? 10
161 ST+ X
162 FC?C 10
163 CHS
164 RCL 10
165 X^2
166 *
167 ST/ 11
168 ST/ 12

169 RCL 13
170 STO IND 15
171 RCL 12
172 RCL 11
173 RTN
174 GTO 02
175 LBL 00
176 XEQ 03
177 42
178 ST* Z
179 *
180 STO 11
181 X<>Y
182 STO 12
183 XEQ 03
184 6
185 ST* Z
186 *
187 ST- 11
188 X<>Y
189 ST- 12
190 XEQ 03
191 ST+ 11
192 X<>Y
193 ST+ 12
194 XEQ 03
195 8
196 ST/ Z
197 /
198 ST- 11
199 X<>Y
200 ST- 12
201 XEQ 03
202 125
203 ST/ Z
204 /
205 ST+ 11
206 X<>Y
207 ST+ 12
208 25.2
209 GTO 14
210 LBL "ZLAPN"

211 LBL 04
212 STO 24
213 X<>Y
214 STO 10
215 CLX
216 STO 22
217 STO 23
218 LBL 05
219 RCL 10
220 RCL 24
221 ENTER
222 XEQ 02
223 ST+ 22
224 X<>Y
225 ST+ 23
226 DSE 24
227 GTO 05
228 RCL 23
229 RCL 22
230 RTN
231 GTO 04
232 LBL 03
233 RCL 10
234 ST+ 21
235 RCL 13
236 RCL 21
237 +
238 STO IND 15
239 XEQ IND 00
240 STO 19
241 X<>Y
242 STO 20
243 RCL 13
244 RCL 21
245 -
246 STO IND 15
247 XEQ IND 00
248 RCL 17
249 -
250 ST+ 19
251 CLX
252 RCL 18

253 -
254 ST+ 20
255 RCL 20
256 RCL 19
257 FS? 10
258 RTN
259 RCL 14
260 RCL 21
261 -
262 STO IND 16
263 XEQ IND 00
264 ST- 19
265 X<>Y
266 ST- 20
267 RCL 13
268 STO IND 15
269 XEQ IND 00
270 ST+ 19
271 X<>Y
272 ST+ 20
273 RCL 14
274 RCL 21
275 +
276 STO IND 16
277 XEQ IND 00
278 ST+ 19
279 X<>Y
280 ST+ 20
281 RCL 21
282 ST+ IND 15
283 XEQ IND 00
284 ST- 19
285 X<>Y
286 ST- 20
287 RCL 14
288 STO IND 16
289 RCL 20
290 RCL 19
291 RTN
292 END

( 485 bytes / SIZE 018-022-025 )

STACK	INPUTS	OUTPUTS
Y	h	Im( ¶f / ¶x_i )
X	i	Re( ¶f / ¶x_i )

Example: f(x) = exp(-x²).Ln(y²+z) x = y = z = 2 + i

01 LBL "T"
02 RCL 04
03 RCL 03
04 Z^2
05 RCL 06
06 RCL 05
07 Z+Z
08 LNZ
09 RCL 02
10 RCL 01
11 Z^2
12 X<>Y
13 CHS
14 X<>Y
15 CHS
16 E^Z
17 Z*Z
18 RTN

"T"    ASTO 00
   2     STO 01   STO 03   STO 05
   1     STO 02   STO 04   STO 06

-With h = 0.1

0.1   ENTER^
    1    XEQ "ZFD"   >>>>   0.469272437                     ---Execution time = 28s---
                               X<>Y   -0.006070241

Whence ¶f / ¶x₁ = 0.469272437 - 0.006070241 i

-Likewise,

0.1 ENTER^
2 R/S >>>> f '_x2 = ¶f / ¶x₂ = -0.011990004 + 0.029115986 i ---Execution time = 30s---

0.1 ENTER^
3 R/S >>>> f '_x3 = ¶f / ¶x₃ = 0.000513598 +0.007022198 i ---Execution time = 31s---

Notes:

-Though it's relatively fast, it may take a long time when this routine is called many times as it is the case in the Riemaniann section.

-So, if you choose h < 0 , "ZFD" will use a formula of order 4 , often - not always - less accurate but faster

-For instance, with h = -0.01

0.01 CHS ENTER^
1 R/S >>>> f '_x1 = ¶f / ¶x₁ = 0.469272574 - 0.006070254 i ---Execution time = 11s---

3°) Second Derivatives of a Function of N variables ( N < 5 )

-If h > 0 "ZSD" employs a formula of order 10

d²f/dx² = (1/25200.h²).[ -73766 f₀ + 42000.( f₁ + f_-1 ) - 6000.( f₂ + f_-2 ) + 1000.( f₃ + f_-3 ) - 125.( f₄ + f_-4 ) + 8.( f₅ + f_-5 ) ] with f_k = f(x+k.h)

and for the mixed derivative:

          f "_xy = ¶²f / ¶x¶y = (1/50400.h²).[ 73766 f₀₀ + 42000.( f₁₁ + f_-1-1 - f₁₀ - f_-10 - f₀₁ - f_0-1 ) + 6000.( - f₂₂ - f_-2-2 + f₂₀ + f_-20 + f₀₂ + f_0-2 )
                                                         + 1000.( f₃₃ + f_-3-3 - f₃₀ - f_-30 - f₀₃ - f_0-3 ) + 125.( - f₄₄ - f_-4-4 + f₄₀ + f_-40 + f₀₄ + f_0-4 )
                                                               + 8.( f₅₅ + f_-5-5 - f₅₀ - f_-50 - f₀₅ - f_0-5 ) ]

where f_km = f ( x + k.h , y + m.h )

-If h < 0 "ZSD" employs a formula of order 4 ( faster ) i-e

d²f/dx² = (1/12h²).[ - f(x-2h) + 16.f(x-h) - 30.f(x) +16.f(x+h) - f(x+2h) ]

and, for the crossed derivative

f "_xy = ¶²f / ¶x¶y = (1/24.h²).[ 30 f₀₀ + 16.( f₁₁ + f_-1-1 - f₁₀ - f_-10 - f₀₁ - f_0-1 ) + ( - f₂₂ - f_-2-2 + f₂₀ + f_-20 + f₀₂ + f_0-2 ) ]

>>> In both cases, if the function returns an alpha data, 0 is returned without calculating the whole formula.

Data Registers: • R00 = Function name ( Registers R00 thru R2n are to be initialized before executing "ZSD" )

• R01-R02 = x₁ , • R03-R04 = x₂ , .......... , • R2n-1-R2n = x_n n < 5

                                           R10 = h
                                           R11-R12 = ¶²f / ¶x_i¶x_j        R13 ........ R21: temp
Flag: F10
Subroutine:   A program which computes f(x₁, ..... ,x_n) assuming x₁ , ............ , x_n are in R01 thru R2n

STACK	INPUTS	OUTPUTS
Z	h	/
Y	j	Im ( ¶²f / ¶x_i¶x_j )
X	i	Re ( ¶²f / ¶x_i¶x_j )

Example: f(x) = exp(-x²).Ln(y²+z) x = y = z = 2 + i

01 LBL "T"
02 RCL 04
03 RCL 03
04 Z^2
05 RCL 06
06 RCL 05
07 Z+Z
08 LNZ
09 RCL 02
10 RCL 01
11 Z^2
12 X<>Y
13 CHS
14 X<>Y
15 CHS
16 E^Z
17 Z*Z
18 RTN

"T"    ASTO 00
   2     STO 01   STO 03   STO 05
   1     STO 02   STO 04   STO 06

-With h = 0.1

0.1 ENTER^
   1    ENTER^
         XEQ "ZSD" >>>> -1.702735593                   ---Execution time = 31s---
                             X<>Y -1.010546610

Whence f "_xx = ¶²f / ¶x² = -1.702735593 - 1.010546610 i

Likewise,

0.1   ENTER^
   1    ENTER^
   2       R/S     >>>>   f "_xy = ¶²f / ¶x¶y = 0.106191979 - 0.092483912 i                         ---Execution time = 86s---

Notes:

-Like "ZFD" the formula is of order 10 if h > 0 and of order 4 if h < 0

-With h = -0.01

-0.01 ENTER^
1 ENTER^ R/S >>>> f "_xx = ¶²f / ¶x² = -1.702735267 - 1.010547033 i ---Execution time = 14s---

   -0.01   ENTER^
       1      ENTER^
       2         R/S     >>>>   f "_xy = ¶²f / ¶x¶y = 0.106191683 - 0.092483871 i                    ---Execution time = 37s---

4°) Third Derivatives of a Function of N variables ( N < 5 )

"ZTD" employs the following formulae ( of order 10 , 11 , 11 respectively )

• h³ f'''_xxx ~ [ -e f(x-5h) - d f(x-4h) - c f(x-3h) - b f(x-2h) - a f(x-h) + a f(x+h) + b f(x+2h) + c f(x+3h) + d f(x+4h) + e f(x+5h) ]

with a = -1669/720 , b = 8738/5040 , c = -4869/10080 , d = 2522/30240 , e = -205/30240

• h³ f'''_xxy ~ SUM_k=1,2,...,5 A_k [ - 2 f_0k + 2 f_0-k + ( f_kk - f_-k-k + f_-kk - f_k-k ) ]

With A₁ = 5/6 , A₂ = -5/84 , A₃ = +5/756 , A₄ = -5/8064 , A₅ = +1/31500

• h³ f'''_xyz ~ SUM_k=1,2,...,5 A_k [ f_k00 - f_-k00 + f_0k0 - f_0-k0 + f_00k - f_00-k + f_kkk - f_-k-k-k - ( f_kk0 - f_-k-k0 + f_k0k - f_-k0-k + f_0kk - f_0-k-k ) ]

With A₁ = 5/6 , A₂ = -5/84 , A₃ = +5/756 , A₄ = -5/8064 , A₅ = +1/31500

Data Registers: • R00 = Function name ( Registers R00 thru R2n are to be initialized before executing "ZTD" )

• R01-R02 = x₁ , • R03-R04 = x₂ , .......... , • R2n-1-R2n = x_n n < 5

                                           R10 = h
                                           R11-R12 = ¶²f / ¶x_i¶x_j¶x_k       R13 ........ R21: temp
Flags: F09-F10
Subroutine:   A program which computes f(x₁, ..... ,x_n) assuming x₁ , ............ , x_n are in R01 thru R2n

-Lines 54-99-267 are three-byte GTOs

01 LBL "ZTD"
02 X>Y?
03 X<>Y
04 RDN
05 X>Y?
06 X<>Y
07 R^
08 X>Y?
09 X<>Y
10 CF 09
11 CF 10
12 X=Y?
13 SF 09
14 ST+ X
15 DSE X
16 STO 16
17 RDN
18 X=Y?
19 SF 10
20 ST+ X
21 DSE X
22 STO 17
23 RDN
24 ST+ X
25 DSE X
26 STO 18
27 X<>Y
28 ABS
29 STO 10
30 6
31 *
32 STO 15
33 FC? 09
34 GTO 00
35 RCL 17
36 X<> 18
37 STO 17
38 GTO 02
39 LBL 00
40 FC? 10
41 GTO 02
42 RCL 16
43 X<> 17
44 STO 16
45 LBL 02
46 RCL IND 16
47 STO 20

48 RCL IND 17
49 STO 13
50 RCL IND 18
51 STO 14
52 FS? 09
53 FC? 10
54 GTO 01
55 XEQ 05
56 205
57 ST* Z
58 *
59 STO 11
60 X<>Y
61 STO 12
62 XEQ 05
63 2522
64 ST* Z
65 *
66 ST- 11
67 X<>Y
68 ST- 12
69 XEQ 05
70 14607
71 ST* Z
72 *
73 ST+ 11
74 X<>Y
75 ST+ 12
76 XEQ 05
77 52428
78 ST* Z
79 *
80 ST- 11
81 X<>Y
82 ST- 12
83 XEQ 05
84 70098
85 ST* Z
86 *
87 X<>Y
88 RCL 12
89 +
90 X<>Y
91 RCL 11
92 +
93 30240
94 ST/ Z

95 /
96 RCL 20
97 STO IND 16
98 RDN
99 GTO 14
100 LBL 05
101 RCL 10
102 ST- 15
103 RCL 20
104 RCL 15
105 +
106 STO IND 16
107 XEQ IND 00
108 STO 13
109 X<>Y
110 STO 14
111 RCL 20
112 RCL 15
113 -
114 STO IND 16
115 XEQ IND 00
116 X<>Y
117 RCL 14
118 -
119 X<>Y
120 RCL 13
121 -
122 RTN
123 LBL 07
124 RCL 20
125 RCL 15
126 +
127 STO IND 16
128 RCL 13
129 LASTX
130 +
131 STO IND 17
132 XEQ IND 00
133 RCL 15
134 SIGN
135 ST* Z
136 *
137 ST+ 19
138 X<>Y
139 ST+ 21
140 RCL 13
141 RCL 15

142 -
143 STO IND 17
144 XEQ IND 00
145 RCL 15
146 SIGN
147 ST* Z
148 *
149 ST- 19
150 X<>Y
151 ST- 21
152 RCL 20
153 STO IND 16
154 XEQ IND 00
155 RCL 15
156 SIGN
157 ST+ X
158 ST* Z
159 *
160 ST+ 19
161 X<>Y
162 ST+ 21
163 RCL 13
164 STO IND 17
165 RCL 15
166 CHS
167 STO 15
168 X<0?
169 GTO 07
170 RCL 21
171 RCL 19
172 RTN
173 LBL 03
174 CLX
175 STO 19
176 STO 21
177 RCL 10
178 ST- 15
179 FC? 09
180 FS? 10
181 GTO 07
182 LBL 04
183 RCL 20
184 RCL 15
185 +
186 STO IND 16
187 XEQ IND 00
188 RCL 15

189 SIGN
190 ST* Z
191 *
192 ST+ 19
193 X<>Y
194 ST+ 21
195 RCL 13
196 RCL 15
197 +
198 STO IND 17
199 XEQ IND 00
200 RCL 15
201 SIGN
202 ST* Z
203 *
204 ST- 19
205 X<>Y
206 ST- 21
207 RCL 14
208 RCL 15
209 +
210 STO IND 18
211 XEQ IND 00
212 RCL 15
213 SIGN
214 ST* Z
215 *
216 ST+ 19
217 X<>Y
218 ST+ 21
219 RCL 13
220 STO IND 17
221 XEQ IND 00
222 RCL 15
223 SIGN
224 ST* Z
225 *
226 ST- 19
227 X<>Y
228 ST- 21
229 RCL 20
230 STO IND 16
231 XEQ IND 00
232 RCL 15
233 SIGN
234 ST* Z
235 *

236 ST+ 19
237 X<>Y
238 ST+ 21
239 RCL 13
240 RCL 15
241 +
242 STO IND 17
243 XEQ IND 00
244 RCL 15
245 SIGN
246 ST* Z
247 *
248 ST- 19
249 X<>Y
250 ST- 21
251 RCL 14
252 STO IND 18
253 XEQ IND 00
254 RCL 15
255 SIGN
256 ST* Z
257 *
258 ST+ 19
259 X<>Y
260 ST+ 21
261 RCL 13
262 STO IND 17
263 RCL 15
264 CHS
265 STO 15
266 X<0?
267 GTO 04
268 RCL 21
269 RCL 19
270 RTN
271 LBL 01
272 XEQ 03
273 157500
274 ST/ Z
275 /
276 STO 11
277 X<>Y
278 STO 12
279 XEQ 03
280 8064
281 ST/ Z
282 /

283 ST- 11
284 X<>Y
285 ST- 12
286 XEQ 03
287 756
288 ST/ Z
289 /
290 ST+ 11
291 X<>Y
292 ST+ 12
293 XEQ 03
294 84
295 ST/ Z
296 /
297 ST- 11
298 X<>Y
299 ST- 12
300 XEQ 03
301 6
302 ST/ Z
303 /
304 X<>Y
305 RCL 12
306 +
307 X<>Y
308 RCL 11
309 +
310 5
311 ST* Z
312 *
313 LBL 14
314 RCL 10
315 ENTER
316 X^2
317 *
318 ST/ Z
319 /
320 STO 11
321 X<>Y
322 STO 12
323 X<>Y
324 END

( 522 bytes / SIZE 022 )

STACK	INPUTS	OUTPUTS
T	h	/
Z	i	/
Y	j	Im ( ¶²f / ¶x_i¶x_j¶x_k )
X	k	Re ( ¶²f / ¶x_i¶x_j¶x_k )

Example: f(x,y,z,t) = exp(-x²) Ln( y² + z + t³ ) x = y = z = t = 2 + i

-Load this short routine:

01 LBL "T"
02 RCL 08
03 RCL 07
04 RCL 08
05 RCL 07
06 Z^2
07 Z*Z
08 RCL 04
09 RCL 03
10 Z^2
11 Z+Z
12 RCL 06
13 RCL 05
14 Z+Z
15 LNZ
16 RCL 02
17 RCL 01
18 Z^2
19 X<>Y
20 CHS
21 X<>Y
22 CHS
23 E^Z
24 Z*Z
25 RTN

"T"    ASTO 00
   2     STO 01   STO 03   STO 05   STO 07
   1     STO 02   STO 04   STO 06   STO 08

• f'''_xyz = ¶³f / ¶x¶y¶z = ?

-With h = 0.1

     0.1 ENTER^
      1    ENTER^
      2    ENTER^
      3    XEQ "ZTD"   >>>>    f'''_xyz = ¶³f / ¶x¶y¶z = 0.002045420 + 0.002544508 i                           ---Execution time = 3m49s---

• f'''_xtt = ¶³f / ¶x¶t² = ?

-With h = 0.1

     0.1 ENTER^
      1    ENTER^
      4    ENTER^    R/S   >>>>    f'''_xtt = ¶³f / ¶x¶t² = - 0.028361907 - 0.027466199 i                                ---Execution time = 1m43s---

• f'''_yyy = ¶³f / ¶y³ = ?

-With h = 0.1

0.1 ENTER^
2 ENTER^ ENTER^ R/S >>>> f'''_yyy = ¶³f / ¶y³ = - 0.002342127 - 0.001495556 i ---Execution time = 35s---

5°) Hessian Matrix

"ZHESS" calculates and sores the coefficients of the Hessian matrix H

H = [ h_i,j ] with h_i,j = ¶²f / ¶x_i¶x_j

Data Registers: • R00 = Function name ( Registers R00 thru R2n are to be initialized before executing "ZHESS" )

• R01-R02 = x₁ , • R03-R04 = x₂ , .......... , • R2n-1-R2n = x_n

R09 = n < 5 R10 = h R11 to R24 : temp R27 ...... = H

Flag: F10
Subroutine: a program which computes f(x₁, ..... ,x_n) assuming x₁ , ............ , x_n are in R01 thru R2n

01 LBL "ZHESS"
02 STO 09
03 STO 24
04 X^2
05 ST+ X
06 26
07 +
08 STO 26
09 STO 22
10 X<>Y
11 STO 10

12 LBL 01
13 RCL 22
14 RCL 09
15 ST+ X
16 E-5
17 *
18 +
19 STO 23
20 RCL 24
21 STO 25
22 LBL 02

23 RCL 10
24 RCL 24
25 RCL 25
26 XROM "ZSD"
27 X<>Y
28 STO IND 22
29 STO IND 23
30 RCL 23
31 1
32 ST- 22
33 -

34 R^
35 STO IND 22
36 STO IND Y
37 DSE 22
38 DSE 23
39 DSE 25
40 GTO 02
41 RCL 09
42 RCL 24
43 -
44 1

45 +
46 ST+ X
47 ST- 22
48 DSE 24
49 GTO 01
50 RCL 26
51 .1
52 %
53 27
54 +
55 END

( 98 bytes / SIZE 2n^2+27 )

STACK	INPUTS	OUTPUTS
Y	h	/
X	n < 5	37.eee

Example: f(x,y,z,t) = x² y³ z⁴ t⁵ x = y = z = t = 1 + 2 i

01 LBL "T"
02 RCL 08
03 RCL 07
04 5
05 XROM "Z^X"
06 RCL 06
07 RCL 05
08 Z^2
09 Z^2
10 Z*Z
11 RCL 04
12 RCL 03
13 Z*Z
14 RCL 04
15 RCL 03
16 Z^2
17 Z*Z
18 RCL 02
19 RCL 01
20 Z^2
21 Z*Z
22 RTN

"T" ASTO 00

1 STO 01 STO 03 STO 05 STO 07
2 STO 02 STO 04 STO 07 STO 08

if you choose h = 0.1

0.1 ENTER^
4 XEQ "ZHESS" >>>> 27.058 ---Execution time = 16m12s---

-And we find the Hessian matrix H in registers R27 to R58 - after rounding some values:

      23506 + 20592 i        70518 + 61776 i         94024 + 82368 i        117530 + 102960 i
      70518 + 61776 i        70518 + 61776 i       141036 + 123552 i      176295 + 154440 i
      94024 + 82368 i      141036 + 123552 i     141036 + 123552 i      235060 + 205920 i
    117530 + 102960 i    176295 + 154440 i     235060 + 205920 i      235060 + 205920 i

Note:

-Choosing h < 0 like h = -0.01 would give faster results.

6°) Laplacian of a Function of N variables ( N < 5 )

"ZLAPN" may be used to evaluate the Laplacian of a function of n variables, provided n < 5

Data Registers: • R00 = Function name ( Registers R00 thru R2n are to be initialized before executing "ZLAPN" )

• R01-R02 = x₁ , • R03-R04 = x₂ , .......... , • R2n-1-R2n = x_n

R09 = n < 5 R10 = h R22-R23 = Lap(f)

Flag: F10
Subroutine: a program which computes f(x₁, ..... ,x_n) assuming x₁ , ............ , x_n are in R01 thru R2n

STACK	INPUTS	OUTPUTS
Y	h	Im ( Lap(f) )
X	n	Re ( Lap(f) )

where n is the number of variables

Example: f(x,y,z,t) = exp(-x²) Ln( y² + z + t³ ) x = y = z = t = 2 + i

"T"    ASTO 00
   2     STO 01   STO 03   STO 05   STO 07
   1     STO 02   STO 04   STO 06   STO 08

-With h = 0.1

0.1 ENTER^
4 XEQ "ZLAPN" >>>> Lap(f) = ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z² + ¶²f / ¶t² = -2.479704858 - 1.481631606 i ---Execution time = 149s---

Note:

-Here again, choosing h < 0 would give faster results.

7°) Biharmonic Operator ( 3-Dim )

-The Biharmonic Operator ( also called Bilaplacian ) is defined as

D² f = ¶⁴f / ¶x⁴ + ¶⁴f / ¶y⁴ + ¶⁴f / ¶z⁴ + 2 ( ¶⁴f / ¶x²¶y² + ¶⁴f / ¶x²¶z² + ¶⁴f / ¶y²¶z² )

-"ZBHRM" uses the following approximation

   h⁴ D² f ~ A f₀₀₀
               + B₁ ( f₁₀₀ + f_-100 + f₀₁₀ + f_0-10 + f₀₀₁ + f_00-1 ) + C₁ ( f₁₁₀ + f_-1-10 + f_1-10 + f_-110 + f₁₀₁ + f_-10-1 + f_10-1 + f_-101 + f₀₁₁ + f_0-1-1 + f_01-1 + f_0-11 )
               + B₂ ( f₂₀₀ + f_-200 + f₀₂₀ + f_0-20 + f₀₀₂ + f_00-2 ) + C₂ ( f₂₂₀ + f_-2-20 + f_2-20 + f_-220 + f₂₀₂ + f_-20-2 + f_20-2 + f_-202 + f₀₂₂ + f_0-2-2 + f_02-2 + f_0-22 )
               + B₃ ( f₃₀₀ + f_-300 + f₀₃₀ + f_0-30 + f₀₀₃ + f_00-3 ) + C₃ ( f₃₃₀ + f_-3-30 + f_3-30 + f_-330 + f₃₀₃ + f_-30-3 + f_30-3 + f_-303 + f₀₃₃ + f_0-3-3 + f_03-3 + f_0-33 )
               + B₄ ( f₄₀₀ + f_-400 + f₀₄₀ + f_0-40 + f₀₀₄ + f_00-4 ) + C₄ ( f₄₄₀ + f_-4-40 + f_4-40 + f_-440+ f₄₀₄ + f_-40-4 + f_40-4 + f_-404 + f₀₄₄ + f_0-4-4 + f_04-4 + f_0-44 )
               + B₅ ( f₅₀₀ + f_-500 + f₀₅₀+ f_0-50 + f₀₀₅ + f_00-5 ) + C₅ ( f₅₅₀ + f_-5-50 + f_5-50 + f_-550 + f₅₀₅ + f_-50-5 + f_50-5 + f_-505 + f₀₅₅ + f_0-5-5 + f_05-5 + f_0-55 )

where f_ijk = f ( x+i.h , y+j.h , z+k.h ) and

     A = 41523361 / 540000 ,   B₁ = -4069 / 180             C₁ = +10 / 3
                                                 B₂ = +4969 / 1260          C₂ = -10 / 84
                                                 B₃ = -2201 / 3240           C₃ = +10 / 1134
                                                 B₄ = +371 / 4320            C₄ = -10 / 16128
                                                 B₅ = -5221 / 945000       C₅ = +2 / 78750

Data Registers: • R00 = Function name ( Registers R00 thru R06 are to be initialized before executing "ZBHRM" )

• R01-R02 = x • R03-R04 = y • R05-R06 = z

R07-R08 = D² f R10 = h R11 to R20: temp

Flags: /
Subroutine: A program which computes f(x,y,z) assuming x is in X-register, y is in Y-register and z is in Z-register upon entry.

01 LBL "ZBHRM"
02 ABS
03 STO 09
04 STO 10
05 RCL 01
06 STO 12
07 RCL 03
08 STO 13
09 RCL 05
10 STO 14
11 XEQ IND 00
12 6
13 ST* 09
14 ST* Z
15 *
16 STO 19
17 X<>Y
18 STO 20
19 XEQ 01
20 24
21 ST* Z
22 *
23 5221
24 ST* 15
25 ST* 16
26 CLX
27 RCL 16
28 RCL 15
29 Z-Z
30 125
31 ST/ Z
32 /

33 STO 07
34 X<>Y
35 STO 08
36 XEQ 01
37 75
38 ST* Z
39 *
40 10388
41 ST* 15
42 ST* 16
43 CLX
44 RCL 16
45 RCL 15
46 Z-Z
47 16
48 ST/ Z
49 /
50 ST- 07
51 X<>Y
52 ST- 08
53 XEQ 01
54 200
55 ST* Z
56 *
57 15407
58 ST* 15
59 ST* 16
60 CLX
61 RCL 16
62 RCL 15
63 Z-Z
64 3

65 ST/ Z
66 /
67 ST+ 07
68 X<>Y
69 ST+ 08
70 6
71 ST/ 07
72 ST/ 08
73 XEQ 01
74 150
75 ST* Z
76 *
77 4969
78 ST* 15
79 ST* 16
80 CLX
81 RCL 16
82 RCL 15
83 Z-Z
84 ST- 07
85 X<>Y
86 ST- 08
87 7
88 ST/ 07
89 ST/ 08
90 XEQ 01
91 600
92 ST* Z
93 *
94 4069
95 ST* 15
96 ST* 16

97 CLX
98 RCL 16
99 RCL 15
100 Z-Z
101 RCL 08
102 RCL 07
103 Z+Z
104 GTO 04
105 LBL 01
106 RCL 10
107 ST- 09
108 CLX
109 STO 15
110 STO 16
111 STO 17
112 STO 18
113 LBL 02
114 RCL 14
115 RCL 09
116 +
117 STO 05
118 XEQ IND 00
119 ST+ 15
120 X<>Y
121 ST+ 16
122 RCL 14
123 STO 05
124 RCL 13
125 RCL 09
126 +
127 STO 03
128 XEQ IND 00

129 ST+ 15
130 X<>Y
131 ST+ 16
132 RCL 13
133 STO 03
134 RCL 12
135 RCL 09
136 +
137 STO 01
138 XEQ IND 00
139 ST+ 15
140 X<>Y
141 ST+ 16
142 RCL 12
143 STO 01
144 RCL 09
145 CHS
146 STO 09
147 X<0?
148 GTO 02
149 STO 11
150 LBL 03
151 RCL 14
152 RCL 09
153 +
154 STO 05
155 RCL 13
156 RCL 11
157 +
158 STO 03
159 XEQ IND 00
160 ST+ 17

161 X<>Y
162 ST+ 18
163 RCL 13
164 STO 03
165 RCL 12
166 RCL 11
167 +
168 STO 01
169 XEQ IND 00
170 ST+ 17
171 X<>Y
172 ST+ 18
173 RCL 14
174 STO 05
175 RCL 13
176 RCL 09
177 +
178 STO 03
179 XEQ IND 00
180 ST+ 17
181 X<>Y
182 ST+ 18
183 RCL 13
184 STO 03
185 RCL 12
186 STO 01
187 RCL 11
188 CHS
189 STO 11
190 X<0?
191 GTO 03
192 RCL 09

193 CHS
194 STO 09
195 X<0?
196 GTO 03
197 RCL 18
198 RCL 20
199 ST- 16
200 ST+ X
201 -
202 RCL 17
203 RCL 19
204 ST- 15
205 ST+ X
206 -
207 RTN
208 LBL 04
209 180
210 RCL 10
211 X^2
212 X^2
213 *
214 ST/ Z
215 /
216 STO 07
217 X<>Y
218 STO 08
219 X<>Y
220 END

( 350 bytes / SIZE 021 )

STACK	INPUTS	OUTPUTS
Y	/	Im ( D² f )
X	h	Re ( D² f )

Example: f(x) = exp(-x²).Ln(y²+z) x = y = z = 2 + i

01 LBL "T"
02 RCL 04
03 RCL 03
04 Z^2
05 RCL 06
06 RCL 05
07 Z+Z
08 LNZ
09 RCL 02
10 RCL 01
11 Z^2
12 X<>Y
13 CHS
14 X<>Y
15 CHS
16 E^Z
17 Z*Z
18 RTN

"T"    ASTO 00
   2     STO 01   STO 03   STO 05
   1     STO 02   STO 04   STO 06

-If you choose h = 0.1

0.1 XEQ "ZBHRM" >>>> D² f = 13.75488667 - 29.72425444 i ---Execution time = 4m18s---

-The exact result is 13.75496303 - 29.72413228 i rounded to 8 decimals

Notes:

-91 evaluations of the function are performed.
-The formula is of order 10 but - due to cancellation of leading digits - there are seldom more than 4 or 5 significant digits.

8°) Curl, Divergence, Gradient & Laplacian

"ZCDGL" computes the Curl , divergence, gradient and Laplacian of a 3D-Vector Field F= ( f , g , h ) = ( X , Y , Z )

      • If the coordinates are rectangular, clear flags F01 & F02, the HP41 displays REC during the calculations
      • If the coordinates are cylindrical, set F01 and clear F02, the HP41 displays CYL during the calculations
      • If the coordinates are spherical, set F02, the HP41 displays SPH during the calculations

-Flag F01 is cleared if flag F02 is set.

Formulae:

-The formulae for the vector Laplacian in cylindrical & spherical coordinates are given by wolframalpha and may be found here
-In rectangular coordinates, the coordinates of the vector Laplacian are simply the Laplacians of the coordinates.

• Rectangular Coordinates x , y , z

Curl E = ( ¶h/¶y - ¶g/¶z , ¶f/¶z - ¶h/¶x , ¶g/¶x - ¶f/¶y )

Div E = ¶f/¶x + ¶g/¶y + ¶h/¶z

Grad f = ( ¶f/¶x , ¶f/¶y , ¶f/¶z ) and similar formulae for g & h

Lapl f = ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z² and similar formulae for g & h

• Cylindrical Coordinates r , f , z

Curl E = [ (1/r) ¶h/¶f - ¶g/¶z , ¶f/¶z - ¶h/¶r , (1/r) ( g + r ¶g/¶r - ¶f/¶f ) ]

Div E = ¶f/¶r + (1/r) f + (1/r) ¶g/¶f + ¶h/¶z

Grad f = ( ¶f/¶r , (1/r) ¶f/¶f , ¶f/¶z ) and similar formulae for g & h

Lapl f = ¶²f / ¶r² + (1/r) ¶f/¶r + (1/r²) ¶²f / ¶f² + ¶²f / ¶z² and similar formulae for g & h

• Spherical Coordinates r , q , f

Curl E = [ (1/(r.sin q) ) ( h cos q + sin q ¶h/¶q - ¶g/¶f ) , (1/(r.sin q) ) ¶f/¶f - ¶h/¶r - h / r , (1/r) ( g + r ¶g/¶r - ¶f/¶q ) ]

Div E = ¶f/¶r + (2/r) f + g (cos q)/(r sin q) + (1/r) ¶g/¶q + (1/(r sin q)) ¶h/¶f

Grad f = ( ¶f/¶r , (1/r) ¶f/¶q , (1/(r sin q)) ¶f/¶f ) and similar formulae for g & h

Lapl f = ¶²f / ¶r² + (2/r) ¶f/¶r + (1/r²) ¶²f / ¶q² + (1/(r²tan q)) ¶f / ¶q + (1/(r²sin²q)) ¶²f / ¶f² and similar formulae for g & h

Data Registers: R00 = temp R10 = h

                                        R01-R02 = x         R31 to R36 = Curl F     R39 to R44 = Grad f       R57 to R62 = Vect Lapl    R63-R64 = Lap f
                                        R03-R04 = y                                              R45 to R50 = Grad g                                                R65-R66 = Lap g
                                        R05-R06 = z         R37-R38 = Div F          R51 to R56 = Grad h                                                R67-R68 = Lap h

R69 to R74: temp ( if the coordinates are not spherical, R71-R72-R73-R74 are unused )

Flags: /
Subroutines: 3 programs named "X" "Y" "Z" which compute X(x,y,z) , Y(x,y,z) and Z(x,y,z)
assuming x , y , z are in registers R01-R02 , R03-R04 , R05-R06

-Lines 12-25-217-316-435 are three-byte GTOs

01 LBL "ZCDGL"
02 STO 10
03 FS? 02
04 CF 01
05 "REC"
06 FS? 01
07 "CYL"
08 FS? 02
09 "SPH"
10 AVIEW
11 FC? 02
12 GTO 00
13 RCL 04
14 RCL 03
15 XROM "SINZ"
16 STO 71
17 X<>Y
18 STO 72
19 RCL 04
20 RCL 03
21 XROM "TANZ"
22 STO 73
23 X<>Y
24 STO 74
25 GTO 00
26 LBL 04
27 RCL 28
28 RCL 27
29 RCL 02
30 RCL 01
31 Z/Z
32 RCL 20
33 FS? 02
34 ST+ X
35 RCL 19
36 FS? 02
37 ST+ X
38 Z+Z
39 RCL 02
40 RCL 01
41 Z/Z
42 RCL 26
43 RCL 25
44 Z+Z
45 FS? 02
46 GTO 04
47 RCL 30
48 RCL 29
49 Z+Z
50 STO 11
51 X<>Y
52 STO 12
53 X<>Y
54 RTN
55 LBL 04
56 STO 11
57 X<>Y
58 STO 12
59 RCL 30
60 RCL 29
61 RCL 72
62 RCL 71
63 Z^2
64 Z/Z
65 STO 13
66 X<>Y
67 STO 14
68 RCL 22
69 RCL 21
70 RCL 74
71 RCL 73
72 Z/Z
73 RCL 14
74 RCL 13
75 Z+Z
76 RCL 02
77 RCL 01

78 Z^2
79 Z/Z
80 ST+ 11
81 X<>Y
82 ST+ 12
83 RCL 12
84 RCL 11
85 RTN
86 LBL 05
87 RCL 30
88 RCL 29
89 RCL 28
90 RCL 27
91 Z+Z
92 RCL 26
93 RCL 25
94 Z+Z
95 RTN
96 LBL 12
97 RCL 10
98 1
99 ENTER
100 XROM "ZSD"
101 STO 25
102 X<>Y
103 STO 26
104 RCL 17
105 STO 69
106 RCL 18
107 STO 70
108 RCL 10
109 2
110 ST/ 69
111 ST/ 70
112 ENTER
113 XROM "ZSD"
114 STO 27
115 X<>Y
116 STO 28
117 RCL 10
118 3
119 ENTER
120 XROM "ZSD"
121 STO 29
122 X<>Y
123 STO 30
124 RCL 10
125 1
126 XROM "ZFD"
127 STO 19
128 X<>Y
129 STO 20
130 RCL 10
131 2
132 XROM "ZFD"
133 STO 21
134 X<>Y
135 STO 22
136 RCL 10
137 3
138 XROM "ZFD"
139 STO 23
140 X<>Y
141 STO 24
142 RTN
143 LBL 06
144 RCL 24
145 ST+ X
146 RCL 23
147 ST+ X
148 RCL 02
149 RCL 01
150 Z^2
151 Z/Z
152 RCL 72
153 RCL 71
154 Z/Z

155 RTN
156 LBL 07
157 RCL 70
158 RCL 69
159 RCL 02
160 RCL 01
161 Z^2
162 Z/Z
163 RTN
164 LBL 08
165 RCL 70
166 RCL 69
167 RCL 02
168 RCL 01
169 Z/Z
170 RTN
171 LBL 09
172 RCL 24
173 RCL 23
174 RCL 02
175 RCL 01
176 Z/Z
177 RCL 72
178 RCL 71
179 Z/Z
180 RTN
181 LBL 00
182 "X"
183 ASTO 00
184 XEQ 12
185 STO 34
186 STO 44
187 X<>Y
188 STO 33
189 STO 43
190 RCL 21
191 STO 41
192 CHS
193 STO 35
194 RCL 22
195 STO 42
196 CHS
197 STO 36
198 RCL 19
199 STO 39
200 STO 37
201 RCL 20
202 STO 40
203 STO 38
204 CLX
205 STO 59
206 STO 60
207 STO 61
208 STO 62
209 FC? 01
210 FS? 02
211 GTO 01
212 XEQ 05
213 STO 63
214 X<>Y
215 STO 64
216 X<>Y
217 GTO 03
218 LBL 01
219 RCL 22
220 RCL 21
221 RCL 02
222 RCL 01
223 Z/Z
224 STO 41
225 CHS
226 STO 35
227 X<>Y
228 STO 42
229 CHS
230 STO 36
231 FC? 02

232 GTO 02
233 XEQ 09
234 STO 33
235 STO 43
236 X<>Y
237 STO 34
238 STO 44
239 LBL 02
240 XEQ 08
241 FS? 02
242 ST+ X
243 ST+ 37
244 X<>Y
245 FS? 02
246 ST+ X
247 ST+ 38
248 XEQ 04
249 STO 63
250 X<>Y
251 STO 64
252 RCL 22
253 ST+ X
254 RCL 21
255 ST+ X
256 RCL 02
257 RCL 01
258 Z^2
259 Z/Z
260 STO 59
261 X<>Y
262 STO 60
263 FS? 01
264 GTO 01
265 XEQ 06
266 STO 61
267 X<>Y
268 STO 62
269 LBL 01
270 XEQ 07
271 RCL 12
272 RCL 11
273 R^
274 FS? 02
275 ST+ X
276 R^
277 FS? 02
278 ST+ X
279 Z-Z
280 LBL 03
281 STO 57
282 X<>Y
283 STO 58
284 "Y"
285 ASTO 00
286 XEQ 12
287 CHS
288 STO 32
289 X<>Y
290 CHS
291 STO 31
292 RCL 19
293 ST+ 35
294 STO 45
295 RCL 20
296 ST+ 36
297 STO 46
298 RCL 23
299 STO 49
300 RCL 24
301 STO 50
302 FC? 01
303 FS? 02
304 GTO 01
305 RCL 21
306 STO 47
307 ST+ 37
308 RCL 22

309 STO 48
310 ST+ 38
311 XEQ 05
312 STO 65
313 X<>Y
314 STO 66
315 X<>Y
316 GTO 03
317 LBL 01
318 XEQ 08
319 ST+ 35
320 X<>Y
321 ST+ 36
322 FS? 01
323 GTO 01
324 XEQ 09
325 STO 49
326 CHS
327 STO 31
328 X<>Y
329 STO 50
330 CHS
331 STO 32
332 RCL 70
333 RCL 69
334 RCL 74
335 RCL 73
336 Z/Z
337 LBL 01
338 RCL 22
339 RCL 21
340 FS? 02
341 Z+Z
342 RCL 02
343 RCL 01
344 Z/Z
345 ST+ 37
346 X<>Y
347 ST+ 38
348 RCL 22
349 RCL 21
350 RCL 02
351 RCL 01
352 Z/Z
353 STO 47
354 X<>Y
355 STO 48
356 XEQ 04
357 STO 65
358 X<>Y
359 STO 66
360 X<>Y
361 FS? 01
362 GTO 01
363 XEQ 06
364 RCL 74
365 RCL 73
366 Z/Z
367 ST+ 61
368 X<>Y
369 ST+ 62
370 RCL 70
371 RCL 69
372 RCL 74
373 RCL 73
374 Z/Z
375 LBL 01
376 RCL 22
377 RCL 21
378 FS? 02
379 Z+Z
380 RCL 02
381 RCL 01
382 Z^2
383 Z/Z
384 ST+ X
385 ST- 57

386 X<>Y
387 ST+ X
388 ST- 58
389 XEQ 07
390 FS? 01
391 GTO 01
392 RCL 72
393 RCL 71
394 Z^2
395 Z/Z
396 LBL 01
397 RCL 12
398 RCL 11
399 R^
400 R^
401 Z-Z
402 LBL 03
403 ST+ 59
404 X<>Y
405 ST+ 60
406 "Z"
407 ASTO 00
408 XEQ 12
409 FS? 01
410 GTO 01
411 FS? 02
412 GTO 02
413 STO 56
414 ST+ 38
415 X<>Y
416 STO 55
417 ST+ 37
418 RCL 21
419 STO 53
420 ST+ 31
421 RCL 22
422 STO 54
423 ST+ 32
424 RCL 19
425 STO 51
426 ST- 33
427 RCL 20
428 STO 52
429 ST- 34
430 XEQ 05
431 STO 67
432 X<>Y
433 STO 68
434 X<>Y
435 GTO 03
436 LBL 01
437 ST+ 38
438 X<>Y
439 ST+ 37
440 RCL 19
441 STO 51
442 ST- 33
443 RCL 20
444 STO 52
445 ST- 34
446 RCL 23
447 STO 55
448 RCL 24
449 STO 56
450 GTO 01
451 LBL 02
452 XEQ 08
453 STO 11
454 X<>Y
455 STO 12
456 X<>Y
457 RCL 20
458 RCL 19
459 Z+Z
460 ST- 33
461 X<>Y
462 ST- 34

463 RCL 12
464 RCL 11
465 RCL 74
466 RCL 73
467 Z/Z
468 ST+ 31
469 X<>Y
470 ST+ 32
471 RCL 24
472 RCL 23
473 RCL 02
474 RCL 01
475 Z/Z
476 RCL 72
477 RCL 71
478 Z/Z
479 STO 55
480 ST+ 37
481 X<>Y
482 STO 56
483 ST+ 38
484 RCL 19
485 STO 51
486 RCL 20
487 STO 52
488 LBL 01
489 RCL 22
490 RCL 21
491 RCL 02
492 RCL 01
493 Z/Z
494 STO 53
495 ST+ 31
496 X<>Y
497 STO 54
498 ST+ 32
499 XEQ 04
500 STO 67
501 X<>Y
502 STO 68
503 X<>Y
504 FS? 01
505 GTO 03
506 XEQ 06
507 ST- 57
508 X<>Y
509 ST- 58
510 X<>Y
511 RCL 74
512 RCL 73
513 Z/Z
514 ST- 59
515 X<>Y
516 ST- 60
517 XEQ 07
518 RCL 72
519 RCL 71
520 Z^2
521 Z/Z
522 RCL 12
523 RCL 11
524 R^
525 R^
526 Z-Z
527 LBL 03
528 ST+ 61
529 X<>Y
530 ST+ 62
531 57.068
532 39.056
533 37.038
534 CLD
535 31.036
536 END

( 963 bytes / SIZE 075 )

STACK	INPUTS	OUTPUTS
T	/	57.068
Z	/	39.056
Y	/	37.038
X	h	31.036

    31.036 = control number of the Curl ( 3 complex numbers )
    37.038 = --------------------- Divergence = 1 complex number
    39 056 = --------------------- Gradients = 3 x 3 complex numbers
    57.068 = --------------------- Laplacians = 1 complex number for the vector Laplacian & 3 x 1 complex numbers for the scalar Laplacians

Example1 - Rectangular Coordinates: CF 01 CF 02

F = ( f , g , h ) = [ exp(-x²) Ln(y²+z) , x²y²z² , exp(x) y²z ] With x = y = z = 1 + 2 i

-Load the routines:

01 LBL "X"
02 RCL 04
03 RCL 03
04 Z^2
05 RCL 06
06 RCL 05
07 Z+Z
08 LNZ
09 RCL 02
10 RCL 01
11 Z^2
12 X<>Y
13 CHS
14 X<>Y
15 CHS
16 E^Z
17 Z*Z
18 RTN
19 LBL "Y"
20 RCL 02
21 RCL 01
22 RCL 04
23 RCL 03
24 Z*Z
25 RCL 06
26 RCL 05
27 Z*Z
28 Z^2
29 RTN
30 LBL "Z"
31 RCL 02
32 RCL 01
33 E^Z
34 RCL 04
35 RCL 03
36 Z^2
37 Z*Z
38 RCL 06
39 RCL 05
40 Z*Z
41 RTN

CF 01 CF 02

1 STO 01 STO 03 STO 05
2 STO 02 STO 04 STO 06

-If you choose h = 0.1

0.1 XEQ "ZCDGL" >>>>        31.036 = control number of the Curl                                           ---Execution time = 6m43s---
                                   RDN         37.038 = --------------------- Divergence
                                   RDN         39 056 = --------------------- Gradients
                                   RDN         57.068 = --------------------- Laplacians

-So, Curl F = ( -94.98658723 + 52.12000491 i , -14.45014465 + 26.13586270 i , 80.96399935 - 90.16478376 i )

and in R37-E38 Div F = 194.2339651 + 117.6139838 i

-You also find in registers R39 to R56:

    Grad f = (118.7272587 + 205.5539813 i , 1.036000652 + 14.16478376 i , 2.936556771 + 1.209278241 i )
    Grad g = ( 82 - 76 i , 82 - 76 i , 82 - 76 i )
    Grad h = ( 17.38670142 - 24.92658446 i , -12.98658723 - 23.87999509 i , -6.493293575 - 11.93999755 i )

in registers R57 to R62

Lapl F = ( 688.9653703 - 896.0414166 i , -42 - 144 i , 5.237391160 - 24.50795524 i )

and in R63 thru R68:

    Lapl f = 688.9653703 - 896.0414166 i
    Lapl g = -42 - 144 i
    Lapl h = 5.237391160 - 24.50795524 i

-In cylindrical & spherical coordinates, the vector Laplacian is less easy to compute...

Example2 - Cylindrical Coordinates: SF 01 CF 02

F = ( f , g , h ) = ( r z² sin²f , r² z , r³ z cos f ) r = 2 + i , f = PI/5 + PI/3 i , z = 1 + 2 i

01 LBL "X"
02 RCL 04
03 RCL 03
04 XROM "SINZ"
05 RCL 06
06 RCL 05
07 Z*Z
08 Z^2
09 RCL 02
10 RCL 01
11 Z*Z
12 RTN
13 LBL "Y"
14 RCL 02
15 RCL 01
16 Z^2
17 RCL 06
18 RCL 05
19 Z*Z
20 RTN
21 LBL "Z"
22 RCL 04
23 RCL 03
24 XROM "COSZ"
25 RCL 06
26 RCL 05
27 Z*Z
28 RCL 02
29 RCL 01
30 Z*Z
31 RCL 02
32 RCL 01
33 Z^2
34 Z*Z
35 RTN

SF 01 CF 02

2 STO 01 PI 5 / STO 03 1 STO 05
1 STO 02 PI 3 / STO 04 2 STO 06

-If you choose h = 0.1

0.1 XEQ "ZCDGL" >>>>        31.036 = control number of the Curl                                           ---Execution time = 8m56s---
                                   RDN         37.038 = --------------------- Divergence
                                   RDN         39 056 = --------------------- Gradients
                                   RDN         57.068 = --------------------- Laplacians

-So, Curl F = ( 11.81071680 - 8.352454288 i , -21.62580411 - 51.22375785 i , 16.70295144 + 3.026594510 i )

and in R37-E38 Div F = -3.723500330 + 0.268718860 i

-You also find in registers R39 to R56:

    Grad f = ( -7.195386865 - 6.251906933 i , -16.70295144 + 11.97340549 i , -19.01490724 - 1.368587064 i )
    Grad g = ( 0 + 10 i , 0 , 3 + 4 i )
    Grad h = ( 2.610896869 + 49.85517079 i , 14.81071680 - 4.352454288 i , 10.66727337 + 12.77253276 i )

in registers R57 thru 62

Lapl F = ( 11.36372903 + 15.97898944 i , -5.572998956 + 22.25990497 i , 29.37437882 + 51.78637057 i )

and in R63 thru R68:

    Lapl f = 7.235192910 + 14.91730402 i
    Lapl g = 4 + 8 i
    Lapl h = 29.37437882 + 51.78637057 i

Note:

-In cylindrical coordinates, Lapl h = the 3rd coordinate of the vector Laplacian
-This not true in spherical coordinates.

Example2 - Spherical Coordinates: CF 01 SF 02

F = ( f , g , h ) = ( r sin² q cos² f , r²sin f , r³ cos q cos² f ) r = 2 + i , q = PI/4 + PI/7 i , f = PI/5 + PI/3 i

01 LBL "X"
02 RCL 04
03 RCL 03
04 XROM "SINZ"
05 STO 75
06 X<>Y
07 STO 76
08 RCL 06
09 RCL 05
10 XROM "COSZ"
11 RCL 76
12 RCL 75
13 Z*Z
14 Z^2
15 RCL 02
16 RCL 01
17 Z*Z
18 RTN
19 LBL "Y"
20 RCL 06
21 RCL 05
22 XROM "SINZ"
23 RCL 02
24 RCL 01
25 Z^2
26 Z*Z
27 RTN
28 LBL "Z"
29 RCL 06
30 RCL 05
31 XROM "COSZ"
32 Z^2
33 STO 75
34 X<>Y
35 STO 76
36 RCL 04
37 RCL 03
38 XROM "COSZ"
39 RCL 76
40 RCL 75
41 Z*Z
42 RCL 02
43 RCL 01
44 Z*Z
45 RCL 02
46 RCL 01
47 Z^2
48 Z*Z
49 RTN

CF 01 SF 02

2 STO 01 PI 4 / STO 03 PI 5 / STO 05
1 STO 02 PI 7 / STO 04 PI 3 / STO 06

-If you choose h = 0.1

0.1 XEQ "ZCDGL" >>>>        31.036 = control number of the Curl                                           ---Execution time = 15m46s---
                                   RDN         37.038 = --------------------- Divergence
                                   RDN         39 056 = --------------------- Gradients
                                   RDN         57.068 = --------------------- Laplacians

-So, Curl F = ( -10.00203205 - 10.02528911 i , -35.48241590 + 15.81684439 i , 0.985054409 + 11.60674974 i )

and in R37-E38 Div F = -11.27980803 - 8.335798791 i

-You also find in registers R39 to R56:

    Grad f = ( 1.541115319 - 0.369198220 i , 1.626418145 - 2.720319640 i , -2.650373943 - 2.249451987 i )
    Grad g = ( 1.740981702 + 5.924286734 i , 0 , 3.542190827 - 1.714238056 i )
    Grad h = ( 24.62403147 - 13.54972230 i , -8.967281608 - 2.712699366 i , -18.62992945 - 8.676217665 i )

in registers R57 to R62

Lapl F = ( 13.15816600 + 1.756716012 i , 14.47564115 - 10.35413557 i , 29.82975585 - 13.49497748 i )

and in R63 thru R68:

    Lapl f = -1.370425868 + 1.423609606 i
    Lapl g = 3.714085670 + 6.017232138 i
    Lapl h = 32.22058351 - 15.01929182 i

8-b) Curl, Divergence, Gradient & Laplacian#2

-This variant uses the spherical coordinates ( r L b ) where L = longitude & b = latitude.

[ instead of ( r , theta , phi ) with phi = L and theta = PI/2 - b = colatitude ]

-The formulae become, with a vector field ( f , g , h ) [ which would be ( f , -h , g ) with the version above ]

• Spherical Coordinates r , L , b

Curl E = [ 1/(r.cos b) ¶h/¶L + ( g Tan b ) / r - ( 1/r ) ¶g/¶b ) , ( 1/r ) ¶f/¶b - ¶h/¶r - h / r , g / r + ¶g/¶r - 1/(r.cos b) ¶f/¶L ) ]

Div E = ¶f/¶r + (2/r) f - h (tan b)/r + (1/r) ¶h/¶b + (1/(r cos b)) ¶g/¶L

Grad f = ( ¶f/¶r , (1/r.cos b) ¶f/¶L , (1/r) ¶f/¶b ) and similar formulae for g & h

Lapl f = ¶²f / ¶r² + (2/r) ¶f/¶r + (1/r²) ¶²f / ¶b² - (tan b/(r²)) ¶f / ¶b + (1/(r²cos²b)) ¶²f / ¶L² and similar formulae for g & h

Vectot Laplacian =

[ ¶²f / ¶r² + (2/r) ¶f/¶r + (1/r²) ¶²f / ¶b² - (tan b/(r²)) ¶f / ¶b + (1/(r²cos²b)) ¶²f / ¶L² - (2/r²) ¶h/¶b - (2/(r²cos b)) ¶g/¶L - 2 f / r² + ( 2 h tan b ) / r² ;
¶²g / ¶r² + (2/r) ¶g/¶r + (1/r²) ¶²g / ¶b² - (tan b/(r²)) ¶g / ¶b + (1/(r²cos²b)) ¶²g / ¶L² + (2/(r²cos b)) ¶f/¶L - (2 tan b )/(r²cos b)) ¶h/¶L - g / (r²cos b) ;
¶²h / ¶r² + (2/r) ¶h/¶r + (1/r²) ¶²h / ¶b² - (tan b/(r²)) ¶h / ¶b + (1/(r²cos²b)) ¶²h / ¶L² + (2 tan b )/(r²cos b)) ¶g/¶L + (2/r²) ¶f/¶b - h / (r²cos b) ]

"ZCDGL" computes the Curl , divergence, gradient and Laplacian of a 3D-Vector Field F= ( f , g , h ) = ( X , Y , Z )

-Flag F01 is cleared if flag F02 is set.

Data Registers: R00 = temp R10 = h

R69 to R76: temp ( if the coordinates are not spherical, R71 to R76 are unused )

Flags: /
Subroutines: 3 programs named "X" "Y" "Z" which compute X(x,y,z) , Y(x,y,z) and Z(x,y,z)
assuming x , y , z are in registers R01-R02 , R03-R04 , R05-R06

-Lines 14-36-242-340-421-433 are three-byte GTOs

01 LBL "ZCDGL"
02 STO 10
03 FS? 02
04 CF 01
05 "REC"
06 FS? 01
07 "CYL"
08 FS? 02
09 "SPH"
10 AVIEW
11 FC? 01
12 FS? 02
13 FS? 30
14 GTO 00
15 RCL 06
16 RCL 05
17 XROM "COSZ"
18 STO 71
19 X<>Y
20 STO 72
21 X<>Y
22 RCL 02
23 RCL 01
24 FS? 02
25 Z*Z
26 STO 75
27 X<>Y
28 STO 76
29 RCL 06
30 RCL 05
31 FS? 02
32 XROM "TANZ"
33 STO 73
34 X<>Y
35 STO 74
36 GTO 00
37 LBL 04
38 RCL 28
39 RCL 27
40 FS? 01
41 GTO 04
42 RCL 72
43 RCL 71
44 Z^2
45 Z/Z
46 LBL 04
47 RCL 02
48 RCL 01
49 Z/Z
50 RCL 20
51 FS? 02
52 ST+ X
53 RCL 19
54 FS? 02
55 ST+ X
56 Z+Z
57 RCL 02
58 RCL 01
59 Z/Z
60 RCL 26
61 RCL 25
62 Z+Z
63 FS? 02
64 GTO 04
65 RCL 30
66 RCL 29
67 Z+Z
68 STO 11
69 X<>Y
70 STO 12
71 X<>Y
72 RTN
73 LBL 04
74 STO 11

75 X<>Y
76 STO 12
77 RCL 24
78 RCL 23
79 RCL 74
80 RCL 73
81 Z*Z
82 RCL 30
83 RCL 29
84 Z-Z
85 RCL 02
86 RCL 01
87 Z^2
88 Z/Z
89 ST- 11
90 X<>Y
91 ST- 12
92 RCL 12
93 RCL 11
94 RTN
95 LBL 05
96 RCL 30
97 RCL 29
98 RCL 28
99 RCL 27
100 Z+Z
101 RCL 26
102 RCL 25
103 Z+Z
104 RTN
105 LBL 12
106 ASTO 00
107 RCL 10
108 1
109 ENTER
110 XROM "ZSD"
111 STO 25
112 X<>Y
113 STO 26
114 RCL 17
115 STO 69
116 RCL 18
117 STO 70
118 RCL 10
119 2
120 ST/ 69
121 ST/ 70
122 ENTER
123 XROM "ZSD"
124 STO 27
125 X<>Y
126 STO 28
127 RCL 10
128 3
129 ENTER
130 XROM "ZSD"
131 STO 29
132 X<>Y
133 STO 30
134 RCL 10
135 1
136 XROM "ZFD"
137 STO 19
138 X<>Y
139 STO 20
140 RCL 10
141 2
142 XROM "ZFD"
143 STO 21
144 X<>Y
145 STO 22
146 RCL 10
147 3
148 XROM "ZFD"

149 STO 23
150 X<>Y
151 STO 24
152 RTN
153 LBL 06
154 RCL 22
155 ST+ X
156 RCL 21
157 ST+ X
158 RCL 02
159 RCL 01
160 Z/Z
161 RCL 76
162 RCL 75
163 Z/Z
164 RTN
165 LBL 07
166 RCL 22
167 ST+ X
168 RCL 21
169 ST+ X
170 RCL 74
171 RCL 73
172 Z*Z
173 RCL 02
174 RCL 01
175 Z/Z
176 RCL 76
177 RCL 75
178 Z/Z
179 RTN
180 LBL 08
181 RCL 22
182 RCL 21
183 RCL 76
184 RCL 75
185 Z/Z
186 RTN
187 LBL 09
188 RCL 70
189 RCL 69
190 RCL 76
191 RCL 75
192 Z^2
193 Z/Z
194 ST- 11
195 X<>Y
196 ST- 12
197 RCL 12
198 RCL 11
199 RTN
200 LBL 10
201 RCL 70
202 RCL 69
203 RCL 02
204 RCL 01
205 Z/Z
206 RTN
207 LBL 00
208 "X"
209 XEQ 12
210 STO 34
211 STO 44
212 X<>Y
213 STO 33
214 STO 43
215 RCL 21
216 STO 41
217 CHS
218 STO 35
219 RCL 22
220 STO 42
221 CHS
222 STO 36

223 RCL 19
224 STO 39
225 STO 37
226 RCL 20
227 STO 40
228 STO 38
229 CLX
230 STO 59
231 STO 60
232 STO 61
233 STO 62
234 FC? 01
235 FS? 02
236 GTO 01
237 XEQ 05
238 STO 63
239 X<>Y
240 STO 64
241 X<>Y
242 GTO 03
243 LBL 01
244 XEQ 08
245 STO 41
246 CHS
247 STO 35
248 X<>Y
249 STO 42
250 CHS
251 STO 36
252 FS? 01
253 GTO 01
254 RCL 24
255 RCL 23
256 RCL 02
257 RCL 01
258 Z/Z
259 STO 33
260 STO 43
261 X<>Y
262 STO 34
263 STO 44
264 RCL 24
265 ST+ X
266 RCL 23
267 ST+ X
268 RCL 02
269 RCL 01
270 Z^2
271 Z/Z
272 STO 61
273 X<>Y
274 STO 62
275 LBL 01
276 XEQ 06
277 STO 59
278 X<>Y
279 STO 60
280 XEQ 10
281 FS? 02
282 ST+ X
283 ST+ 37
284 X<>Y
285 FS? 02
286 ST+ X
287 ST+ 38
288 CHS
289 X<>Y
290 CHS
291 RCL 02
292 RCL 01
293 Z/Z
294 STO 13
295 X<>Y
296 STO 14

297 XEQ 04
298 STO 63
299 X<>Y
300 STO 64
301 X<>Y
302 RCL 14
303 RCL 13
304 Z+Z
305 LBL 03
306 STO 57
307 X<>Y
308 STO 58
309 "Y"
310 XEQ 12
311 CHS
312 STO 32
313 X<>Y
314 CHS
315 STO 31
316 RCL 19
317 STO 45
318 ST+ 35
319 RCL 20
320 STO 46
321 ST+ 36
322 FC? 01
323 FS? 02
324 GTO 01
325 RCL 21
326 STO 47
327 ST+ 37
328 RCL 22
329 STO 48
330 ST+ 38
331 RCL 23
332 STO 49
333 RCL 24
334 STO 50
335 XEQ 05
336 STO 65
337 X<>Y
338 STO 66
339 X<>Y
340 GTO 03
341 LBL 01
342 XEQ 10
343 ST+ 35
344 X<>Y
345 ST+ 36
346 XEQ 06
347 ST- 57
348 X<>Y
349 ST- 58
350 FS? 02
351 GTO 02
352 RCL 23
353 STO 49
354 RCL 24
355 STO 50
356 RCL 22
357 RCL 21
358 GTO 01
359 LBL 02
360 RCL 70
361 RCL 69
362 RCL 74
363 RCL 73
364 Z*Z
365 ST+ 31
366 X<>Y
367 ST+ 32
368 XEQ 07
369 ST+ 61
370 X<>Y

371 ST+ 62
372 RCL 24
373 RCL 23
374 RCL 02
375 RCL 01
376 Z/Z
377 STO 49
378 X<>Y
379 STO 50
380 RCL 22
381 RCL 21
382 RCL 72
383 RCL 71
384 Z/Z
385 LBL 01
386 RCL 02
387 RCL 01
388 Z/Z
389 ST+ 37
390 X<>Y
391 ST+ 38
392 XEQ 08
393 STO 47
394 X<>Y
395 STO 48
396 XEQ 04
397 STO 65
398 X<>Y
399 STO 66
400 XEQ 09
401 LBL 03
402 ST+ 59
403 X<>Y
404 ST+ 60
405 "Z"
406 XEQ 12
407 FS? 02
408 GTO 02
409 STO 56
410 ST+ 38
411 X<>Y
412 STO 55
413 ST+ 37
414 RCL 19
415 STO 51
416 ST- 33
417 RCL 20
418 STO 52
419 ST- 34
420 FS? 01
421 GTO 01
422 RCL 21
423 STO 53
424 ST+ 31
425 RCL 22
426 STO 54
427 ST+ 32
428 XEQ 05
429 STO 67
430 X<>Y
431 STO 68
432 X<>Y
433 GTO 03
434 LBL 02
435 X<>Y
436 RCL 02
437 RCL 01
438 Z/Z
439 ST+ 37
440 STO 55
441 X<>Y
442 ST+ 38
443 STO 56
444 XEQ 10

445 RCL 20
446 RCL 19
447 Z+Z
448 ST- 33
449 X<>Y
450 ST- 34
451 RCL 32
452 RCL 31
453 RCL 02
454 RCL 01
455 Z/Z
456 STO 31
457 X<>Y
458 STO 32
459 RCL 70
460 RCL 69
461 RCL 74
462 RCL 73
463 Z*Z
464 STO 13
465 X<>Y
466 STO 14
467 X<>Y
468 RCL 02
469 RCL 01
470 Z/Z
471 ST- 37
472 X<>Y
473 ST- 38
474 RCL 19
475 STO 51
476 RCL 20
477 STO 52
478 RCL 14
479 RCL 13
480 RCL 24
481 RCL 23
482 Z-Z
483 RCL 02
484 RCL 01
485 Z^2
486 Z/Z
487 ST+ X
488 ST+ 57
489 X<>Y
490 ST+ X
491 ST+ 58
492 XEQ 07
493 ST- 59
494 X<>Y
495 ST- 60
496 LBL 01
497 XEQ 08
498 STO 53
499 ST+ 31
500 X<>Y
501 STO 54
502 ST+ 32
503 XEQ 04
504 STO 67
505 X<>Y
506 STO 68
507 X<>Y
508 FS? 02
509 XEQ 09
510 LBL 03
511 ST+ 61
512 X<>Y
513 ST+ 62
514 57.068
515 39.056
516 37.038
517 31.036
518 CLD
519 END

( 933 bytes / SIZE 077 )

STACK	INPUTS	OUTPUTS
T	/	57.068
Z	/	39.056
Y	/	37.038
X	h	31.036

Example - Spherical Coordinates: CF 01 SF 02

( The results are the same as above with rectangular & cylindrical coordinates ).

F = ( f , g , h ) = ( r² cos² L cos² b , r³ sin² L sin b , r⁴ sin L cos² b ) r = 1 + 2 i , L = PI/5 + PI/3 i , b = PI/4 - i PI/7

01 LBL "X"
02 RCL 06
03 RCL 05
04 XROM "COSZ"
05 STO 77
06 X<>Y
07 STO 78
08 RCL 04
09 RCL 03
10 XROM "COSZ"
11 RCL 78
12 RCL 77
13 Z*Z
14 RCL 02
15 RCL 01
16 Z*Z
17 Z^2
18 RTN
19 LBL "Y"
20 RCL 04
21 RCL 03
22 XROM "SINZ"
23 RCL 02
24 RCL 01
25 Z*Z
26 Z^2
27 RCL 02
28 RCL 01
29 Z*Z
30 STO 77
31 X<>Y
32 STO 78
33 RCL 06
34 RCL 05
35 XROM "SINZ"
36 RCL 78
37 RCL 77
38 Z*Z
39 RTN
40 LBL "Z"
41 RCL 06
42 RCL 05
43 XROM "COSZ"
44 RCL 02
45 RCL 01
46 Z^2
47 Z*Z
48 Z^2
49 STO 77
50 X<>Y
51 STO 78
52 RCL 04
53 RCL 03
54 XROM "SINZ"
55 RCL 78
56 RCL 77
57 Z*Z
58 RTN

CF 01 SF 02

1 STO 01 PI 5 / STO 03 PI 4 / STO 05
2 STO 02 PI 3 / STO 04 PI 7 / CHS STO 06

-If you choose h = 0.1

0.1 XEQ "ZCDGL" >>>>        31.036 = control number of the Curl                                           ---Execution time = 18m37s---
                                   RDN         37.038 = --------------------- Divergence
                                   RDN         39 056 = --------------------- Gradients
                                   RDN         57.068 = --------------------- Laplacians

-So, Curl F = ( -17.64085209 + 10.08005857 i , -19.50265611 + 53.26019404 i , -32.48966259 - 2.500308528 i )

and in R37-E38 Div F = 23.87095012 + 59.06252866 i

-You also find in registers R39 to R56:

    Grad f = ( 4.559023524 + 5.426064833 i , 1.848530060 - 7.550199864 i , -7.067057392 - 0.532516658 i )
    Grad g = ( -22.98084939 - 7.537881290 i , -3.112788790 + 20.31407327 i , -3.557656590 - 7.234384010 i )
    Grad h = ( 9.948478996 - 43.03416857 i , -14.04716098 - 0.876080263 i , 11.91046124 + 18.59755046 i )

in registers R57 to R62

Lapl F = ( -45.55171103 - 6.743439951 i , -29.18137714 - 16.67563547 i , -52.18565801 - 40.46155333 i )

and in R63 thru R68:

    Lapl f = 2.000000038 + 0.000000732 i
    Lapl g = -20.43061066 + 4.992444368 i
    Lapl h = -73.87863763 - 41.86112025 i

Note:

-This program is not in the ZDERIVE41 module

9°) Gaussian Curvature & Mean Curvature of a Surface

-"ZKGM" calculates the Gaussian curvature K_G and the mean curvature K_m of a surface defined by a function z = z(x,y)

-The Gaussian Curvature K_G is be computed by:

K_G = [ ( ¶²z / ¶x² ) ( ¶²z / ¶y² ) - ( ¶²z / ¶x¶y )² ] / [ 1 + ( ¶z / ¶x )² + ( ¶z / ¶y )² ]²

-And the mean curvature K_m is given by:

K_m = (1/2) [ t ( 1+p² ) + r ( 1+q² ) - 2 p q s ] / ( 1+p²+q² )^3/2

where p = ¶z / ¶x , q = ¶z / ¶y, r = ¶²z / ¶x², s = ¶²z / ¶x¶y, t = ¶²z / ¶y²

Data Registers: • R00 = f name ( Registers R00 thru R04 are to be initialized before executing "ZKGM" )

• R01-R02 = x • R03-R04 = y R10 = h R15 to R28: temp

R11-R12 = K_G R13-R14 = K_m
Flags: /
Subroutine: A program that takes x in R01-R02, y in R03-R04 and returns z(x,y) in X-register

01 LBL "ZKGM"
02 1
03 2
04 XROM "ZSD"
05 STO 27
06 X<>Y
07 STO 28
08 RCL 10
09 2
10 ENTER
11 XROM "ZSD"
12 STO 25
13 X<>Y
14 STO 26
15 RCL 10
16 1
17 ENTER
18 XROM "ZSD"
19 STO 23
20 X<>Y
21 STO 24
22 RCL 10

23 2
24 XROM "ZFD"
25 STO 21
26 X<>Y
27 STO 22
28 RCL 10
29 1
30 XROM "ZFD"
31 STO 19
32 X<>Y
33 STO 20
34 X<>Y
35 Z^2
36 STO 11
37 X<>Y
38 STO 12
39 X<>Y
40 1
41 +
42 RCL 26
43 RCL 25
44 Z*Z

45 STO 13
46 X<>Y
47 STO 14
48 RCL 22
49 RCL 21
50 Z^2
51 1
52 +
53 ST+ 11
54 X<>Y
55 ST+ 12
56 X<>Y
57 RCL 24
58 RCL 23
59 Z*Z
60 ST+ 13
61 X<>Y
62 ST+ 14
63 2
64 ST/ 13
65 ST/ 14
66 RCL 22

67 RCL 21
68 RCL 20
69 RCL 19
70 Z*Z
71 RCL 28
72 RCL 27
73 Z*Z
74 ST- 13
75 X<>Y
76 ST- 14
77 RCL 26
78 RCL 25
79 RCL 24
80 RCL 23
81 Z*Z
82 RCL 28
83 RCL 27
84 Z^2
85 Z-Z
86 RCL 12
87 RCL 11
88 Z^2

89 Z/Z
90 X<> 11
91 X<>Y
92 X<> 12
93 X<>Y
94 1.5
95 XROM "Z^X"
96 RCL 14
97 RCL 13
98 R^
99 R^
100 Z/Z
101 STO 13
102 X<>Y
103 STO 14
104 X<>Y
105 RCL 12
106 RCL 11
107 END

( 175 bytes / SIZE 029 )

STACK	INPUTS	OUTPUTS
T	/	Im K_m
Z	/	Re K_m
Y	/	Im K_G
X	h	Re K_G

Example: A surface is defined by z(x,y) = Ln ( 1 + x² y³ )

-Calculate the Gaussian & mean curvatures at the point (x,y) = (1+2i,1+2i)

01 LBL "T"
02 RCL 04
03 RCL 03
04 RCL 04
05 RCL 03
06 Z^2
07 Z*Z
08 RCL 02
09 RCL 01
10 Z^2
11 Z*Z
12 1
13 +
14 LNZ
15 RTN

"T" ASTO 00

1 STO 01 STO 03
2 STO 02 STO 04

and if you choose h = 0.1

     0.1 XEQ "ZKGM" >>>>   0.034972708 = R11                                                     ---Execution time = 2m32s---
                                     RDN -0.037342284 = R12
                                     RDN -0.116193627 = R13
                                     RDN   0.113513738 = R14

-So, K_G = 0.034972708 - 0.037342284 i
and K_m = -0.116193627 + 0.113513738 i

Note:

-With a sphere, you'll get K = 1 / R² where R = radius of the sphere.

10°) Curvature & Torsion of a Parametric 3D-Curve

-We assume that the curve is parameterized by 3 functions x(t) , y(t) , z(t)

-The curvature is calculated by khi = < r' x r'' , r' x r'' >^1/2 / < r' , r' >^3/2 where the vector r = [ x(t) , y(t) , z(t) ]
and the torsion by tau = < r' x r'' , r''' > / < r' x r'' , r' x r'' > x = cross-product and < , > = dot-product

Data Registers: R00: temp ( Registers R01-R02 are to be initialized before executing "ZKHT" )

                                    • R01-R02 = t           R10 = h                            R22-R23 = x'        R24-R25 = x"     R26-R27 = x'''
                                                                                                               R28-R29 = y'        R30-R31 = y"     R32-R33 = y'''
                                        R11-R12 = khi       R13-R14 = tau                  R34-R35 = z'        R36-R37 = z"     R38-R39 = z'''

Flags: /
Subroutines: 3 functions named "Z1" "Z2" "Z3" that compute x(t) , y(t) , z(t) assuming t is in registers R01-R02

01 LBL "ZKHT"
02 STO 10
03 3
04 STO 41
05 39
06 STO 40
07 FIX 00
08 CF 29
09 LBL 01
10 "Z"
11 ARCL 41
12 ASTO 00
13 RCL 10
14 1
15 ENTER
16 ENTER
17 XROM "ZTD"
18 X<>Y
19 STO IND 40
20 DSE 40
21 X<>Y
22 STO IND 40
23 RCL 10
24 1

25 ENTER
26 XROM "ZSD"
27 DSE 40
28 X<>Y
29 STO IND 40
30 DSE 40
31 X<>Y
32 STO IND 40
33 RCL 10
34 1
35 XROM "ZFD"
36 DSE 40
37 X<>Y
38 STO IND 40
39 DSE 40
40 X<>Y
41 STO IND 40
42 DSE 40
43 DSE 41
44 GTO 01
45 FIX 09
46 SF 29
47 RCL 29
48 RCL 28

49 RCL 37
50 RCL 36
51 Z*Z
52 STO 11
53 X<>Y
54 STO 12
55 RCL 31
56 RCL 30
57 RCL 35
58 RCL 34
59 Z*Z
60 ST- 11
61 X<>Y
62 ST- 12
63 RCL 25
64 RCL 24
65 RCL 35
66 RCL 34
67 Z*Z
68 STO 13
69 X<>Y
70 STO 14
71 RCL 23
72 RCL 22

73 RCL 37
74 RCL 36
75 Z*Z
76 ST- 13
77 X<>Y
78 ST- 14
79 RCL 23
80 RCL 22
81 RCL 31
82 RCL 30
83 Z*Z
84 STO 15
85 X<>Y
86 STO 16
87 RCL 25
88 RCL 24
89 RCL 29
90 RCL 28
91 Z*Z
92 ST- 15
93 X<>Y
94 ST- 16
95 RCL 23
96 RCL 22

97 Z^2
98 RCL 29
99 RCL 28
100 Z^2
101 Z+Z
102 RCL 35
103 RCL 34
104 Z^2
105 Z+Z
106 3
107 XROM "Z^X"
108 STO 19
109 X<>Y
110 STO 20
111 RCL 12
112 RCL 11
113 Z^2
114 RCL 14
115 RCL 13
116 Z^2
117 Z+Z
118 RCL 16
119 RCL 15
120 Z^2

121 Z+Z
122 STO 17
123 X<>Y
124 STO 18
125 X<>Y
126 RCL 20
127 RCL 19
128 Z/Z
129 SQRTZ
130 X<> 11
131 X<>Y
132 X<> 12
133 X<>Y
134 RCL 27
135 RCL 26
136 Z*Z
137 STO 40
138 X<>Y
139 STO 41
140 RCL 33
141 RCL 32
142 RCL 14
143 RCL 13
144 Z*Z

145 RCL 41
146 RCL 40
147 Z+Z
148 STO 13
149 X<>Y
150 STO 14
151 RCL 39
152 RCL 38
153 RCL 16
154 RCL 15
155 Z*Z
156 RCL 14
157 RCL 13
158 Z+Z
159 RCL 18
160 RCL 17
161 Z/Z
162 STO 13
163 X<>Y
164 STO 14
165 X<>Y
166 RCL 12
167 RCL 11
168 END

( 288 bytes / SIZE 042 )

STACK	INPUTS	OUTPUTS
T	/	Im tau
Z	/	Re tau
Y	/	Im khi
X	h	Re khi

Example: The curve defined by x(t) = e^t cos t , y(t) = e^t sin t , z(t) = Ln(t)

-Evaluate the curvature khi (t) and the torsion tau (t) for t = 1 + 2 i

-Load the following routines in main memory

01 LBL "Z1"
02 RCL 02
03 RCL 01
04 XROM "COSZ"
05 RCL 02
06 RCL 01
07 E^Z
08 Z*Z
09 RTN
10 LBL "Z2"
11 RCL 02
12 RCL 01
13 XROM "SINZ"
14 RCL 02
15 RCL 01
16 E^Z
17 Z*Z
18 RTN
19 LBL "Z3"
20 RCL 02
21 RCL 01
22 LNZ
23 RTN

1 STO 01
2 STO 02

and if you choose h = 0.1

    0.1   XEQ "ZKHT" >>>>    0.109325569 = R11                                                     ---Execution time = 4m38s---
                                     RDN   0.234754529 = R12
                                     RDN   0.054205527 = R13
                                     RDN   0.046420520 = R14

-So, khi = 0.109325569 + 0.234754529 i
and tau = 0.054205527 + 0.046420520 i

11°) Curvature & Torsion of a Parametric 4D-Curve

-The method follows the Gram-Schmidt process.
-The formulas are given in references [1] & [2] after replacing | V | by < V , V >^1/2

Data Registers: R00: temp ( Registers R01-R02 are to be initialized before executing "ZKHT4" )

                                    • R01-R02 = t           R10 = h                             R22-R23 = x'        R24-R25 = x"     R26-R27 = x'''      R46 to R73: temp
                                                                                                                R28-R29 = y'        R30-R31 = y"     R32-R33 = y'''
                                        R11-R12 = khi       R13-R14 = tau                  R34-R35 = z'        R36-R37 = z"     R38-R39 = z'''
                                                                                                                R40-R41 = t'        R42-R43 = t"      R44-R45 = t'''
Flags: /
Subroutines:   4 functions named "Z1" "Z2" "Z3" "Z4"   that compute z₁(t) , z₂(t) , z₃(t) , z₄(t) assuming t is in registers R01-R02

01 LBL "ZKHT4"
02 STO 10
03 4
04 STO 59
05 43
06 STO 58
07 FIX 00
08 CF 29
09 GTO 01
10 LBL 11
11 STO M
12 CLST
13 LBL 12
14 RCL IND M
15 DSE M
16 RCL IND M
17 Z^2
18 Z+Z
19 DSE M
20 GTO 12
21 SQRTZ
22 RTN
23 LBL 13
24 CLA
25 STO M
26 X<>Y
27 INT
28 STO N
29 CLST
30 LBL 14
31 RCL IND N
32 DSE N
33 RCL IND N
34 RCL IND M
35 DSE M
36 RCL IND M
37 Z*Z
38 ST+ O
39 X<>Y
40 ST+ P
41 DSE N
42 DSE M
43 GTO 14
44 RCL P
45 RCL O
46 RTN
47 LBL 01
48 "Z"
49 ARCL 59
50 ASTO 00
51 RCL 10

52 1
53 ENTER
54 ENTER
55 XROM "ZTD"
56 X<>Y
57 STO IND 58
58 DSE 58
59 X<>Y
60 STO IND 58
61 7
62 ST- 58
63 RCL 10
64 1
65 ENTER
66 XROM "ZSD"
67 X<>Y
68 STO IND 58
69 DSE 58
70 X<>Y
71 STO IND 58
72 7
73 ST- 58
74 RCL 10
75 1
76 XROM "ZFD"
77 X<>Y
78 STO IND 58
79 DSE 58
80 X<>Y
81 STO IND 58
82 15
83 ST+ 58
84 DSE 59
85 GTO 01
86 FIX 09
87 SF 29
88 27.019
89 STO 05
90 8.008
91 +
92 STO 06
93 LASTX
94 +
95 STO 07
96 LASTX
97 +
98 STO 08
99 LASTX
100 +
101 STO 09
102 RCL 05

103 XEQ 11
104 STO 03
105 X<>Y
106 STO 04
107 RCL 06
108 RCL 05
109 XEQ 13
110 STO 15
111 X<>Y
112 STO 16
113 RCL 06
114 INT
115 STO N
116 RCL 05
117 STO M
118 LBL 02
119 RCL IND M
120 DSE M
121 RCL IND M
122 RCL 16
123 RCL 15
124 Z*Z
125 RCL 04
126 RCL 03
127 Z^2
128 Z/Z
129 RCL IND N
130 DSE N
131 RCL IND N
132 R^
133 R^
134 Z-Z
135 RCL 04
136 RCL 03
137 Z/Z
138 X<>Y
139 STO IND 08
140 DSE 08
141 X<>Y
142 STO IND 08
143 DSE 08
144 DSE N
145 DSE M
146 GTO 02
147 8
148 ST+ 08
149 RCL 08
150 XEQ 11
151 STO 13
152 X<>Y
153 STO 14

154 RCL 06
155 XEQ 11
156 STO 11
157 X<>Y
158 STO 12
159 RCL 07
160 RCL 05
161 XEQ 13
162 CHS
163 STO 00
164 X<>Y
165 CHS
166 STO 17
167 LBL 03
168 RCL 16
169 RCL 15
170 RCL 04
171 RCL 03
172 Z/Z
173 Z^2
174 3
175 ST* Z
176 *
177 RCL 12
178 RCL 11
179 Z^2
180 Z-Z
181 RCL 17
182 RCL 00
183 Z+Z
184 RCL IND 05
185 DSE 05
186 RCL IND 05
187 Z*Z
188 STO M
189 X<>Y
190 STO N
191 RCL 16
192 ST+ X
193 RCL 15
194 ST+ X
195 RCL IND 06
196 CHS
197 DSE 06
198 RCL IND 06
199 CHS
200 Z*Z
201 RCL N
202 RCL M
203 Z+Z
204 RCL 04

205 RCL 03
206 Z^2
207 Z/Z
208 RCL IND 07
209 DSE 07
210 RCL IND 07
211 Z+Z
212 RCL 04
213 RCL 03
214 Z/Z
215 X<>Y
216 STO IND 09
217 DSE 09
218 X<>Y
219 STO IND 09
220 CLX
221 SIGN
222 ST- 09
223 ST- 07
224 ST- 06
225 DSE 05
226 GTO 03
227 8
228 ST+ 05
229 ST+ 07
230 ST+ 09
231 RCL 08
232 RCL 09
233 XEQ 13
234 CHS
235 STO 11
236 X<>Y
237 CHS
238 STO 12
239 LBL 04
240 RCL IND 08
241 DSE 08
242 RCL IND 08
243 RCL 12
244 RCL 11
245 Z*Z
246 RCL 14
247 RCL 13
248 Z^2
249 Z/Z
250 RCL IND 09
251 DSE 09
252 RCL IND 09
253 Z+Z
254 RCL 14
255 RCL 13

256 Z/Z
257 STO IND 09
258 CLX
259 SIGN
260 RCL 09
261 +
262 X<>Y
263 STO IND Y
264 DSE 09
265 CLX
266 DSE 08
267 GTO 04
268 8
269 ST+ 08
270 ST+ 09
271 RCL 07
272 RCL 08
273 XEQ 13
274 CHS
275 STO 11
276 X<>Y
277 CHS
278 STO 12
279 LBL 05
280 RCL IND 05
281 DSE 05
282 RCL IND 05
283 RCL 17
284 RCL 00
285 Z*Z
286 RCL 04
287 RCL 03
288 Z^2
289 Z/Z
290 RCL IND 07
291 DSE 07
292 RCL IND 07
293 Z+Z
294 STO M
295 X<>Y
296 STO N
297 RCL IND 08
298 DSE 08
299 RCL IND 08
300 RCL 12
301 RCL 11
302 Z*Z
303 RCL 14
304 RCL 13
305 Z^2
306 Z/Z

307 RCL N
308 RCL M
309 Z+Z
310 STO IND 08
311 CLX
312 SIGN
313 RCL 08
314 +
315 X<>Y
316 STO IND Y
317 CLX
318 SIGN
319 ST- 08
320 ST- 07
321 DSE 05
322 GTO 05
323 8
324 ST+ 08
325 RCL 08
326 XEQ 11
327 STO 00
328 X<>Y
329 STO 15
330 RCL 14
331 RCL 13
332 RCL 04
333 RCL 03
334 Z/Z
335 STO 11
336 X<>Y
337 STO 12
338 RCL 08
339 RCL 09
340 XEQ 13
341 RCL 04
342 RCL 03
343 Z/Z
344 RCL 15
345 RCL 00
346 Z/Z
347 STO 13
348 X<>Y
349 STO 14
350 X<>Y
351 RCL 12
352 RCL 11
353 CLA
354 END

( 561 bytes / SIZE 074 )

STACK	INPUTS	OUTPUTS
T	/	Im tau
Z	/	Re tau
Y	/	Im khi
X	h	Re khi

Example: The curve defined by

z₁(t) = t⁴ z₃(t) = t³ at the point t = 0.7 + 0.8 i
z₂(t) = t² z₄(t) = e^t

01 LBL "Z1"
02 RCL 02
03 RCL 01
04 Z^2
05 Z^2
06 RTN
07 LBL "Z2"
08 RCL 02
09 RCL 01
10 Z^2
11 RTN
12 LBL "Z3"
13 RCL 02
14 RCL 01
15 Z^2
16 RCL 02
17 RCL 01
18 Z*Z
19 RTN
20 LBL "Z4"
21 RCL 02
22 RCL 01
23 E^Z
24 RTN

0.7 STO 01
0.8 STO 02

and if you choose h = 0.1

    0.1   XEQ "ZKHT4" >>>> -0.124754341 = R11                                                     ---Execution time = 3m31s---
                                       RDN   0.370021139 = R12
                                       RDN   0.143785333 = R13
                                       RDN   0.048055204 = R14

-So, khi = -0.124754341 + 0.370021140 i
and tau = 0.143785368 + 0.048055232 i

12°) Complex Linear Systems

-We need a program to find the inverse of the matrix [ g_ij ] i-e [ g^ij ]
-But "ZLS" may also be used for another purpose

"ZLS" allows you to solve linear systems, including overdetermined and underdetermined systems.
You can also invert up to a 8x8 complex matrix.

The determinant of this matrix is also computed Re (det A) is stored in register R00.

X-register = Re (det A)
Y-register = Im (det A)

Gaussian elimination with partial pivoting is used.

You can choose the 1st data register Rbb - except R00 -
You store the first coefficient into Rbb for the real part, Rbb+1 for the imaginary part , ... , up to the last one into Ree ( column by column ) ( bb > 00 )

Then you key in bbb.eeerr XQ "ZLS" and the system will be solved.

where r is the number of rows of the matrix

>>>> If a pivot is smaller than about 5 E-7 , it will be considered to be zero.

Don't interrupt "ZLS" because registers P and Q are used ( there is no risk of crash, but their contents could be altered )
"ZLS" only employs the data registers where the coefficients are stored and R00

-Line 209 is a three-byte GTO

01 LBL "ZLS"
02 CLA
03 ENTER
04 FRC
05 ISG X
06 INT
07 1
08 %
09 +
10 E5
11 /
12 +
13 ENTER
14 INT
15 LASTX
16 FRC
17 ISG X
18 INT
19 RCL Y
20 +
21 DSE X
22 E3
23 /
24 +
25 STO O
26 SIGN
27 STO 00
28 X<>Y
29 STO N
30 LBL 01
31 RCL N

32 INT
33 RCL O
34 FRC
35 +
36 STO P
37 ENTER
38 CLX
39 ISG P
40 LBL 02
41 RCL IND Y
42 X^2
43 ISG Z
44 RCL IND Z
45 X^2
46 +
47 X<=Y?
48 GTO 02
49 X<>Y
50 X<> Z
51 STO P
52 X<>Y
53 ENTER
54 LBL 02
55 RDN
56 ISG Y
57 GTO 02
58 RCL N
59 ENTER
60 FRC
61 RCL P
62 INT

63 DSE X
64 +
65 X=Y?
66 GTO 04
67 LBL 03
68 RCL IND X
69 X<> IND Z
70 STO IND Y
71 CLX
72 SIGN
73 ST+ Z
74 +
75 RCL IND X
76 X<> IND Z
77 STO IND Y
78 CLX
79 SIGN
80 ST- Z
81 -
82 ISG Y
83 TEXT0
84 ISG X
85 GTO 03
86 SIGN
87 CHS
88 ST* 00
89 ST* M
90 LBL 04
91 CLX
92 SIGN
93 RCL N

94 +
95 RCL IND X
96 RCL IND N
97 RCL Y
98 RCL Y
99 R-P
100 X<>Y
101 CLX
102 PI
103 CHS
104 10^X
105 X^2
106 X>Y?
107 CLST
108 R^
109 R^
110 RCL M
111 RCL 00
112 Z*Z
113 STO 00
114 RDN
115 STO M
116 RDN
117 Z=0?
118 GTO 09
119 1/Z
120 RCL N
121 STO P
122 STO Q
123 SIGN
124 ST+ Q

125 RDN
126 LBL 05
127 RCL IND Q
128 RCL IND P
129 Z*Z
130 STO IND P
131 RDN
132 STO IND Q
133 ISG Q
134 RDN
135 ISG P
136 GTO 05
137 LBL 06
138 RCL N
139 ENTER
140 STO P
141 FRC
142 RCL O
143 INT
144 +
145 X=Y?
146 GTO 08
147 STO Q
148 ENTER
149 SIGN
150 +
151 RCL IND X
152 RCL IND Q
153 LBL 07
154 RCL P
155 ENTER

156 SIGN
157 +
158 RDN
159 RCL IND T
160 RCL IND P
161 Z*Z
162 ST- IND Q
163 CLX
164 RCL Q
165 INT
166 ISG X
167 CLX
168 RDN
169 ST- IND T
170 ISG P
171 RDN
172 ISG Q
173 GTO 07
174 LBL 08
175 ISG O
176 ISG O
177 GTO 06
178 RCL N
179 FRC
180 ISG X
181 INT
182 ST- O
183 LBL 09
184 RCL N
185 FRC
186 ISG X

187 INT
188 .1
189 %
190 +
191 ST+ O
192 2
193 RCL N
194 +
195 ENTER
196 ISG X
197 X=0?
198 GTO 00
199 E7
200 *
201 DSE X
202 STO N
203 LASTX
204 ST/ N
205 CLX
206 E2
207 MOD
208 X#0?
209 GTO 01
210 LBL 00
211 RCL M
212 RCL 00
213 CLA
214 END

( 321 bytes / SIZE var )

STACK	INPUTS	OUTPUTS
Y	/	Im (det A)
X	bbb.eeerr	Re (det A)

Example1: Find the solution of the system:

    ( 6 + 9 i ) x + ( 3 + 7 i ) y + ( 2 + 5 i ) z = -51 + 91 i
    ( 5 + 2 i ) x + ( 7 + 6 i ) y + ( 8 + 4 i ) z = 14 + 126 i
    ( 4 + i ) x + ( 4 + 9 i ) y + ( 1 + 2 i ) z = -29 + 68 i

If you choose to store the first element into R11, you have to store these 24 numbers:

    6   9        3 7      2 5     -51   91                  R11 R12   R17 R18   R23 R24    R29 R30
    5   2        7 6      8 4      14   126    into       R13 R14   R19 R20   R25 R26    R31 R32        respectively
    4   1        4 9      1 2     -29   68                  R15 R16   R21 R22   R27 R28    R33 R34

>>> The first register is R11, the last register is R34 and there are 3 rows, therefore the control number of the matrix is 11.03403

11.03403 XEQ "ZLS" >>>> 453 = R00 ---Execution time = 38s---
RDN -248

>>> So, Det A = 453 - 248 i

Registers R11 thru R28 now contain the unit matrix and the solution is in registers R29 thru R34

    x = 1 + 2 i
    y = 3 + 4 i            with very small roundoff errors in the last decimals
    z = 5 + 6 i

Example2: Invert the same matrix A

         [ [ 6 + 9 i    3 + 7 i    2 + 5 i ]
A = [ 5 + 2 i    7 + 6 i    8 + 4 i ]
           [ 4 + i      4 + 9 i    1 + 2 i ] ]

-If you choose again R11 for the 1st coefficient, store

    6   9        3 7      2 5     1 0     0 0    0 0             R11 R12   R17 R18   R23 R24    R29 R30    R35 R36   R41 R42
    5   2        7 6      8 4     0 0     1 0    0 0    into    R13 R14   R19 R20   R25 R26    R31 R32    R37 R38   R43 R44    respectively
    4   1        4 9      1 2     0 0     0 0    1 0             R15 R16   R21 R22   R27 R28    R33 R34    R39 R40   R45 R46

-The control number of this array is 11.04603

11.04603 XEQ "ZLS" >>>> 453 = R00 ---Execution time = 52s---
RDN -248

>>> So, Det A = 453 - 248 i ( of course the same result as above )

Registers R11 thru R28 now contain the unit matrix and the inverse matrix B = A^-1 is in registers R29 thru R46

           [ [ 0.061530559 - 0.116424771 i    -0.067405788 + 0.018285573 i    0.000854851 + 0.046825614 i ]
   B =    [ 0.034700221 + 0.045487097 i    -0.024546985 - 0.015646032 i    0.041917717 - 0.124954539 i   ]        rounded to 9D
             [ -0.054425544 + 0.018769239 i    0.160164671 - 0.042558855 i    -0.076010543 + 0.086422484 i ] ]

-You can check, after multiplying all the coefficients by det A = 453 - 248 i that

            [ [ -1 - 68 i    -26 + 25 i    12 + 21 i ]
    B =    [ 27 + 12 i    -15 - i       -12 - 67 i   ]   /   det A
              [ -20 + 22 i    62 - 59 i    -13 + 58 i ] ]

13°) Initialization - Metric Tensor & Christoffel Symbols

-Always execute "ZINIGC" before calculating the components of the tensors below

-You have to load in main memory the functions g_ij which must be called LBL "G11" LBL "G12" .... LBL "G22" and so on
-Since the metric tensor is symmetric, LBL "G21" ... are unuseful.

-They must take x¹ , x² , ........... , xⁿ in registers R01-R02 , R03-R04............. R2n-1R2n ( n < 5 ) and return g_ij in X & Y registers

-"ZINIGC" calculates and stores the n(n+1)/2 g_ij in registers R61-R62 ................ ( i <= j ) for a given point ( x¹, ........... , xⁿ )
-Then it calls "ZLS" to find the inverse of the metric matrix and stores g^ij in registers R61+n(n+1) .....

-The determinant g = det g_ij is stored in R59-R60

-The Christoffel symbols G^k_ij = (1/2) g^km ( ¶_i g_mj + ¶_j g_im - ¶_m g_ij ) are then computed and stored in the following registers.

-Only G^k_ij with i <= j because they are symmetric G^k_ji = G^k_ij

>>> After executing "ZINIGC" , do not disturb register R09 & R61 to R60+n(n+1)(n+2)

Data Registers: R00 = temp ( Registers R01 thru R2n are to be initialized before executing "ZINIGC" )

                                      • R01-R02 = x¹            R09 = n                              R61 to R60+n(n+1)                         = g_ij    with i <= j
                                      • R03-R04 = x²            R10 = h                              R61+n(n+1) to R60+2n(n+1)          = g^ij    with i <= j
                                         .............                                                                 R61+2n(n+1) to R60+n(n+1)(n+2) = G^k_ij with i <= j
                                      • R2n-1R2n = xⁿ          R59-60 = g = det g_ij

Flags /
Subroutines: The n(n+1)/2 functions g_ij with i <= j

-Lines 247-249-251 are three-byte GTOs

01 LBL "ZINIGC"
02 STO 09
03 STO 11
04 ENTER
05 X^2
06 +
07 STO 20
08 STO 58
09 X<>Y
10 STO 10
11 CLX
12 60
13 +
14 STO 13
15 RCL 09
16 X^2
17 ST+ X
18 +
19 STO 26
20 RCL 09
21 ST+ X
22 E5
23 /
24 +
25 STO 14
26 CLX
27 STO 15
28 LBL 04
29 RCL 11
30 STO 12
31 LBL 05
32 RCL 11
33 RCL 12
34 -ZDERIVE41
35 ASTO X
36 XEQ IND X
37 SIGN

38 X=0?
39 ENTER
40 X#0?
41 X<> L
42 X<>Y
43 STO IND 13
44 STO IND 14
45 STO IND 26
46 CLX
47 SIGN
48 CHS
49 RCL 14
50 +
51 X<>Y
52 DSE 13
53 DSE 26
54 STO IND 13
55 STO IND 26
56 STO IND Y
57 DSE 26
58 DSE 14
59 DSE 13
60 DSE 12
61 GTO 05
62 2
63 ST+ 15
64 RCL 15
65 ST- 26
66 RCL 14
67 FRC
68 RCL 26
69 +
70 STO 14
71 DSE 11
72 GTO 04
73 RCL 09
74 X^2

75 4
76 *
77 RCL 58
78 +
79 60
80 +
81 .1
82 %
83 +
84 RCL 09
85 X^2
86 ST+ X
87 DSE X
88 -
89 STO 12
90 CLRGX
91 RCL 09
92 1
93 +
94 ST+ X
95 E5
96 /
97 +
98 SIGN
99 LBL 06
100 STO IND L
101 ISG L
102 GTO 06
103 RCL 12
104 FRC
105 61
106 +
107 RCL 58
108 +
109 RCL 09
110 E5
111 /

112 +
113 STO 13
114 XROM "ZLS"
115 STO 59
116 X<>Y
117 STO 60
118 RCL 13
119 FRC
120 E3
121 *
122 INT
123 STO 12
124 RCL 58
125 ST+ X
126 60
127 +
128 STO 11
129 RCL 09
130 STO 13
131 ST+ X
132 STO 14
133 LBL 07
134 RCL IND 12
135 STO IND 11
136 DSE 11
137 DSE 12
138 DSE 14
139 GTO 07
140 DSE 13
141 FS? 30
142 GTO 00
143 RCL 13
144 ST+ X
145 STO 14
146 RCL 09
147 RCL 13
148 -

149 ST+ X
150 ST- 12
151 GTO 07
152 LBL 00
153 RCL 58
154 RCL 09
155 ST* 20
156 STO 23
157 2
158 ST/ 20
159 +
160 *
161 60
162 +
163 STO 25
164 STO 29
165 LBL 10
166 RCL 09
167 STO 21
168 LBL 01
169 RCL 21
170 STO 22
171 LBL 02
172 VIEW 20
173 CLX
174 STO 30
175 STO 31
176 RCL 09
177 STO 24
178 LBL 03
179 RCL 23
180 RCL 24
181 ENTER
182 CLX
183 ZCIJK
184 Z=0?
185 GTO 00

186 STO 28
187 X<>Y
188 STO 32
189 CLX
190 STO 26
191 STO 27
192 RCL 22
193 RCL 24
194 -ZDERIVE41
195 ASTO 00
196 RCL 10
197 RCL 21
198 XROM "ZFD"
199 ST+ 26
200 X<>Y
201 ST+ 27
202 RCL 21
203 RCL 24
204 -ZDERIVE41
205 ASTO 00
206 RCL 10
207 RCL 22
208 XROM "ZFD"
209 ST+ 26
210 X<>Y
211 ST+ 27
212 RCL 21
213 RCL 22
214 -ZDERIVE41
215 ASTO 00
216 RCL 10
217 RCL 24
218 XROM "ZFD"
219 ST- 26
220 X<>Y
221 ST- 27
222 RCL 27

223 RCL 26
224 RCL 32
225 RCL 28
226 Z*Z
227 ST+ 30
228 X<>Y
229 ST+ 31
230 LBL 00
231 DSE 24
232 GTO 03
233 RCL 31
234 RCL 30
235 2
236 ST/ Z
237 /
238 X<>Y
239 STO IND 25
240 DSE 25
241 X<>Y
242 STO IND 25
243 DSE 25
244 DSE 20
245 CLX
246 DSE 22
247 GTO 02
248 DSE 21
249 GTO 01
250 DSE 23
251 GTO 10
252 RCL 29
253 E3
254 /
255 61
256 +
257 CLD
258 END

( 410 bytes / SIZE 061+n(n+1)(n+2) )

STACK	INPUTS	OUTPUTS
Y	h	/
X	n < 5	bbb.eee(gij & G^k_ij)

where n is the number of variables

Example: A 4D-Riemannian manifold, with a metric defined by

    g₁₁ = 1 + x² y²           g₁₄ = x t     and the other    g_ij = 0
    g₂₂ = 1 + x² z²
    g₃₃ = 1 + y² z²
    g₄₄ = 1 + x² t²

-Load for instance these 10 routines in main memory:

01 LBL "G11"
02 RCL 02
03 RCL 01
04 RCL 04
05 RCL 03
06 Z*Z
07 Z^2
08 ENTER
09 CLX
10 SIGN
11 +
12 RTN

13 LBL "G22"
14 RCL 02
15 RCL 01
16 RCL 06
17 RCL 05
18 Z*Z
19 Z^2
20 ENTER
21 CLX
22 SIGN
23 +
24 RTN

25 LBL "G33"
26 RCL 06
27 RCL 05
28 RCL 04
29 RCL 03
30 Z*Z
31 Z^2
32 ENTER
33 CLX
34 SIGN
35 +
36 RTN

37 LBL "G44"
38 RCL 02
39 RCL 01
40 RCL 08
41 RCL 07
42 Z*Z
43 Z^2
44 ENTER
45 CLX
46 SIGN
47 +
48 RTN

49 LBL "G14"
50 RCL 02
51 RCL 01
52 RCL 08
53 RCL 07
54 Z*Z
55 RTN
56 LBL "G12"
57 ASTO X
58 RTN
59 LBL "G13"
60 ASTO X

61 RTN
62 LBL "G23"
63 ASTO X
64 RTN
65 LBL "G24"
66 ASTO X
67 RTN
68 LBL "G34"
69 ASTO X
70 RTN
71 END

>>> At the point x = y = z = t = 1 + 2 i

1 STO 01 STO 03 STO 05 STO 07
2 STO 02 STO 04 STO 06 STO 08

SIZE 181 at least

-We have n = 4 and if you choose h = - 0.01 ( h = +0.1 would be much slower )

- 0.01   ENTER^
      4      XEQ "ZINIGC"   >>>>   After calculating the inverse matrix [ g^ij ]
                                                      the HP41 displays a countdown   40 39 ........... 2   1

and finally returns the control number 61.180 ---Execution time = 12m08s---

-And you find in registers R61 thru R80

    R61-62 = g₁₁ = -6 - 24 i     R63-64 = g₁₂ = 0                R67-68 = g₁₃ = 0                 R73-74 = g₁₄ = -3 + 4 i
                                               R65-66 = g₂₂ = -6 - 24 i     R69-70 = g₂₃ = 0                 R75-76 = g₂₄ = 0
                                                                                          R71-72 = g₃₃ = -6 - 24 i     R77-78 = g₃₄ = 0
                                                                                                                                      R79-80 = g₄₄ = -6 - 24 i
-In registers R81 thru R100

R81-82 = g¹¹ = -0.011247 +0.038444 i rounded to 6D ............... R93-94 = g¹⁴ = -0.007464 + 0.003136 i

and so on .... until R99-R100 = g⁴⁴ = -0.011247 + 0.038444 i rounded to 6D

-And in registers R101 to R180, for example

R101-R102 = G¹₁₁ = 0.186872 - -0.412188 i

R161-R162 = G⁴₁₁ = 0.000238574 - 0.003612692 i .... and so on .... until R179-R180 = G⁴₄₄ = 0.098496984 - 0.392624655 i

-We have also R59-R60 = g = det g_ij = 197964 - 321984 i

Notes:

-Use ZCIJK below to recall g_ij , g^ij and G^k_ij

-Execute a SIZE at least

n	SIZE
2	85
3	121
4	181

-Registers R50 to R56 are unused by all these programs,
-So, you can use them to calculate the different Gij ( and synthetic registers M N O too )

-If n = 4, h = +0.1 and with relatively simple g_ij functions, all different from zero, V41 needs about 7 seconds to XEQ "ZINIGC"

14°) Recalling Gij & Cijk

-ZCIJK is an M-Code routine that recalls G^k_ij provided R09 = n
-This register has been initialized by "ZINIGC"

08B "K"
00A "J"
009   "I"
003   "C"
01A "Z"
378   C=c
03C RCR 3
106 A=C S&X
130 LDI S&X
009 009
146 A=A+C S&X
130 LDI S&X
200 200h
306 ?A<C S&X
381 GOTO
00A NEXIST
0A6 A<>C S&X
270 RAMSLCT
038 READATA
10E A=C ALL
046 C=0 S&X
270 RAMSLCT
0AE A<>C ALL
361   chk alpha
050   and SETDEC
10E A=C ALL
070 N=C ALL
04E C
35C =
050 1
01D C=
060 A+C
0B0 C=N ALL
13D C=
060 AB*C
070 N=C ALL
0F8 C=X
13D C=
060 AB*C
0B0 C=N ALL
025 C=
060 AB+C
070 N=C ALL
078 C=Z
10E A=C ALL
0B8 C=Y
30E ?A<C ALL
027 JC+04
068 Z=C
0AE A<>C ALL
0A8 Y=C
10E   A=C ALL
135 C=
060 A*C
0B8 C=Y
2BE C=-C
025 C=
061 AB+C
0B0 C=N ALL
025 C=
060 AB+C
078 C=Z
025 C=
060 AB+C
078 C=Z
025 C=
060 AB+C
068 Z=C
260 SETHEX
38D C
008 ->S&X
106 A=C S&X
130 LDI S&X
03C 03Ch=060d
146 A=A+C S&X
378 C=c
03C RCR 3
146 A=A+C S&X
130 LDI S&X
200 200h
306 ?A<C S&X
381 GOTO
00A NEXIST
0A6 A<>C S&X
270 RAMSLCT
106 A=C S&X
038 READATA
0EE B<>C ALL
0A6 A<>C S&X
266 C=C-1 S&X
270 RAMSLCT
038 READATA
10E A=C ALL
046 C=0 S&X
270 RAMSLCT
0AE A<>C ALL
0E8 X=C
0EE B<>C ALL
0A8 Y=C
3E0 RTN

STACK	INPUTS	OUTPUTS
Z	i	address-60
Y	j	Im G^k_ij
X	k	Re G^k_ij

Where "address" is the address of the data register that contains Im G^k_ij

Example1:

-With the metric above and after executing "ZINIGC"

      3 ENTER^
          ENTER^
      2 XEQ "ZCIJK" >>>> -0.186274510
                                  RDN    0.411764706

So, G²₃₃ = -0.186274510 + 0.411764706 i

-Register Z also contains 72 so the address of the data registers which contains G²₃₃ are R131 & R132

Example2: ZCIJK may also be used to recall g_ij if k = -1 and g^ij if k = 0

-For instance, 4 ENTER^ ENTER^ 0 ZCIJK gives g⁴⁴ = -0.011247060 + 0.038444497 i

and 2 ENTER^ ENTER^ 1 CHS ZCIJK yields g₂₂ = - 6 - 24 i

Notes:

-ZCIJK checks that the data register R09 does exist and that it does contain numbers
-However, it does not check X-Y-Z-registers for alpha data.

-Since n is a positive integer smaller than 10, ZCIJK only take into account the first digit of the mantissas

-The formula to get the address of the register Rmm that contains Im G^k_ij is

m = 60 + 2 i + ( j² - j ) + ( k + 1 ) n ( n + 1 )

15°) Curvature Tensor ( 0-4 )

-The curvature tensor R_ijkl ( 4 times covariant ) may be computed by the formula:

R_ijkl = (1/2) ( ¶_jk g_il + ¶_il g_jk - ¶_jl g_ik - ¶_ik g_jl ) + g_mp ( G^m_jkG^p_il - G^m_jlG^p_ik )

-This listing combines "ZRIJ" & "ZRIJKL"
-Lines 174-176 are three-byte GTOs

01 LBL "ZRIJ"
02 LBL 12
03 SF 09
04 STO 24
05 X<>Y
06 STO 22
07 RCL 09
08 STO 23
09 CLX
10 STO 36
11 STO 37
12 LBL 04
13 RCL 09
14 STO 25
15 LBL 05
16 RCL 23
17 RCL 25
18 ENTER
19 CLX
20 ZCIJK
21 Z=0?
22 GTO 00
23 STO 34
24 X<>Y
25 STO 35
26 GTO 13
27 LBL 11
28 RCL 35
29 RCL 34
30 Z*Z

31 ST+ 36
32 X<>Y
33 ST+ 37
34 LBL 00
35 DSE 25
36 GTO 05
37 DSE 23
38 GTO 04
39 RCL 37
40 RCL 36
41 RTN
42 GTO 12
43 LBL "ZRIJKL"
44 LBL 10
45 CF 09
46 STO 25
47 RDN
48 STO 24
49 RDN
50 STO 23
51 X<>Y
52 STO 22
53 LBL 13
54 RCL 09
55 STO 26
56 CLX
57 STO 30
58 STO 31
59 LBL 01
60 RCL 09

61 STO 27
62 LBL 02
63 RCL 26
64 RCL 27
65 ENTER
66 SIGN
67 CHS
68 ZCIJK
69 Z=0?
70 GTO 00
71 STO 32
72 X<>Y
73 STO 33
74 RCL 23
75 RCL 24
76 RCL 26
77 ZCIJK
78 STO 28
79 X<>Y
80 STO 29
81 RCL 22
82 RCL 25
83 RCL 27
84 ZCIJK
85 RCL 29
86 RCL 28
87 Z*Z
88 STO 28
89 X<>Y
90 STO 29

91 RCL 23
92 RCL 25
93 RCL 26
94 ZCIJK
95 STO M
96 X<>Y
97 STO N
98 RCL 22
99 RCL 24
100 RCL 27
101 ZCIJK
102 RCL N
103 RCL M
104 Z*Z
105 RCL 29
106 RCL 28
107 R^
108 R^
109 Z-Z
110 RCL 33
111 RCL 32
112 Z*Z
113 ST+ 30
114 X<>Y
115 ST+ 31
116 LBL 00
117 DSE 27
118 GTO 02
119 DSE 26
120 GTO 01

121 RCL 22
122 RCL 25
123 -ZDERIVE41
124 ASTO 00
125 RCL 10
126 RCL 23
127 RCL 24
128 XROM "ZSD"
129 STO 28
130 X<>Y
131 STO 29
132 RCL 23
133 RCL 24
134 -ZDERIVE41
135 ASTO 00
136 RCL 10
137 RCL 22
138 RCL 25
139 XROM "ZSD"
140 ST+ 28
141 X<>Y
142 ST+ 29
143 RCL 22
144 RCL 24
145 -ZDERIVE41
146 ASTO 00
147 RCL 10
148 RCL 23
149 RCL 25
150 XROM "ZSD"

151 ST- 28
152 X<>Y
153 ST- 29
154 RCL 23
155 RCL 25
156 -ZDERIVE41
157 ASTO 00
158 RCL 10
159 RCL 22
160 RCL 24
161 XROM "ZSD"
162 ST- 28
163 X<>Y
164 ST- 29
165 RCL 29
166 RCL 28
167 2
168 ST/ Z
169 /
170 RCL 31
171 RCL 30
172 Z+Z
173 FS? 09
174 GTO 11
175 RTN
176 GTO 10
177 END

( 325 bytes / SIZE 061+n(n+1)(n+2) )

STACK	INPUTS	OUTPUTS
T	i	/
Z	j	/
Y	k	Im R_ijkl
X	l	Re R_ijkl

Examples: The same metric and with x = y = z = t = 1 + 2 i

1 ENTER^
2 ENTER^
1 ENTER^
2 XEQ "ZRIJKL" >>>> R₁₂₁₂ = -2.759976 + 4.320320 i ---Execution time = 38s---

-Likewise

1 ENTER^
2 ENTER^
2 ENTER^
4 R/S >>>> R₁₁₁₄ = 1.011247 - 0.038444 i ---Execution time = 43s---

Notes:

-The curvature tensor has many components but it has also many symmetries
-So there are in fact N = n² ( n² - 1 ) / 12 independant components i-e

   N = 1 if n = 2
   N = 6 if n = 3
   N =20 if n = 4

16°) Curvature Tensor ( 1-3 )

-The 1-time contravariant 3-times covariant curvature tensor is given by:

R^l_ijk = g^lm R_mijk

-"ZRIJK^L" calculates these components

01 LBL "ZRIJK^L"
02 STO 34
03 RDN
04 STO 25
05 RDN
06 STO 24
07 X<>Y

08 STO 23
09 RCL 09
10 STO 22
11 CLX
12 STO 37
13 STO 38
14 LBL 03

15 RCL 22
16 RCL 34
17 ENTER
18 CLX
19 ZCIJK
20 Z=0?
21 GTO 00

22 STO 35
23 X<>Y
24 STO 36
25 RCL 22
26 RCL 23
27 RCL 24
28 RCL 25

29 XROM "ZRIJKL"
30 RCL 36
31 RCL 35
32 Z*Z
33 ST+ 37
34 X<>Y
35 ST+ 38

36 LBL 00
37 DSE 22
38 GTO 03
39 RCL 38
40 RCL 37
41 END

( 81 bytes / SIZE 061+n(n+1)(n+2) )

STACK	INPUTS	OUTPUTS
T	i	/
Z	j	/
Y	k	Im R^l_ijk
X	l	Re R^l_ijk

Examples: The same metric and with x = y = z = t = 1 + 2 i

2 ENTER^
1 ENTER^
2 ENTER^
1 XEQ "ZRIJK^L" >>>> R¹₂₁₂ = -0.127624 - 0.158155 i ---Execution time = 84s---

-Likewise

3 ENTER^
1 ENTER^
2 ENTER^
4 R/S >>>> R⁴₃₁₂ = 0.008407 - 0.040034 i ---Execution time = 76s---

17°) Ricci Tensor

-If we contract 2 indices of the curvature tensor ( the 2nd and the 4th ), we get the Ricci tensor which is symmetric

R_ij = g^km R_ikjm

-We can also contract other indices, but it would yield -R_ij

>>> The listing is given 2 paragraphs above ( with "ZRIJKL" )

STACK	INPUTS	OUTPUTS
Y	i	Im R_ij
X	j	Re R_ij

Examples: The same metric and with x = y = z = t = 1 + 2 i

1 ENTER^
1 XEQ "ZRIJ" >>>> R₁₁ = -0.135054 - 0.155065 i ---Execution time = 5m11s--

-Likewise

2 ENTER^
3 R/S >>>> R₂₃ = -0.127501 - 0.157186 i ---Execution time = 4m15s---

18°) Scalar Riemannian Curvature

-Contracting the Ricci tensor, we get the scalar curvature R = g^ij R_ij
-"ZR" calculates and stores R in registers R57-R58

-The listing combines "ZR" "ZEIJ" & "ZWIJKL"
-Line 217 is a three-byte GTO

01 LBL "ZR"
02 CLX
03 STO 57
04 STO 58
05 RCL 09
06 STO 22
07 LBL 01
08 RCL 22
09 STO 24
10 LBL 02
11 RCL 22
12 RCL 24
13 ENTER
14 CLX
15 ZCIJK
16 Z=0?
17 GTO 00
18 STO 38
19 X<>Y
20 STO 39
21 RCL 22
22 RCL 24
23 X=Y?
24 GTO 02
25 2
26 ST* 38
27 ST* 39
28 RDN
29 LBL 02
30 XROM "ZRIJ"
31 RCL 39
32 RCL 38

33 Z*Z
34 ST+ 57
35 X<>Y
36 ST+ 58
37 LBL 00
38 DSE 24
39 GTO 02
40 DSE 22
41 GTO 01
42 RCL 58
43 RCL 57
44 RTN
45 LBL "ZEIJ"
46 XEQ 12
47 LBL 10
48 STO 41
49 X<>Y
50 STO 40
51 X<>Y
52 1
53 CHS
54 ZCIJK
55 RCL 58
56 RCL 57
57 Z*Z
58 2
59 ST/ Z
60 /
61 STO 38
62 X<>Y
63 STO 39
64 RCL 40

65 RCL 41
66 XROM "ZRIJ"
67 RCL 39
68 RCL 38
69 Z-Z
70 CLD
71 RTN
72 GTO 10
73 LBL "ZWIJKL"
74 XEQ 12
75 LBL 13
76 STO 41
77 RDN
78 STO 40
79 RDN
80 STO 39
81 X<>Y
82 STO 38
83 CLX
84 STO 42
85 STO 43
86 RCL 39
87 RCL 40
88 1
89 CHS
90 ZCIJK
91 STO 44
92 X<>Y
93 STO 45
94 Z=0?
95 GTO 00
96 RCL 38

97 RCL 41
98 XROM "ZRIJ"
99 RCL 45
100 RCL 44
101 Z*Z
102 ST+ 42
103 X<>Y
104 ST+ 43
105 LBL 00
106 RCL 38
107 RCL 41
108 1
109 CHS
110 ZCIJK
111 RCL 45
112 RCL 44
113 Z*Z
114 STO 44
115 RDN
116 STO 45
117 RDN
118 Z=0?
119 GTO 00
120 STO 46
121 X<>Y
122 STO 47
123 RCL 39
124 RCL 40
125 XROM "ZRIJ"
126 RCL 47
127 RCL 46
128 Z*Z

129 ST+ 42
130 X<>Y
131 ST+ 43
132 LBL 00
133 RCL 39
134 RCL 41
135 1
136 CHS
137 ZCIJK
138 STO 48
139 X<>Y
140 STO 49
141 Z=0?
142 GTO 00
143 RCL 38
144 RCL 40
145 XROM "ZRIJ"
146 RCL 49
147 RCL 48
148 Z*Z
149 ST- 42
150 X<>Y
151 ST- 43
152 LBL 00
153 RCL 38
154 RCL 40
155 1
156 CHS
157 ZCIJK
158 RCL 49
159 RCL 48
160 Z*Z

161 STO 48
162 RDN
163 STO 49
164 RDN
165 Z=0?
166 GTO 00
167 STO 46
168 X<>Y
169 STO 47
170 RCL 39
171 RCL 41
172 XROM "ZRIJ"
173 RCL 47
174 RCL 46
175 Z*Z
176 ST- 42
177 X<>Y
178 ST- 43
179 LBL 00
180 RCL 49
181 RCL 48
182 RCL 45
183 RCL 44
184 Z-Z
185 RCL 58
186 RCL 57
187 Z*Z
188 RCL 09
189 1
190 -
191 ST/ Z
192 /

193 RCL 43
194 RCL 42
195 Z+Z
196 RCL 09
197 2
198 -
199 ST/ Z
200 /
201 STO 42
202 X<>Y
203 STO 43
204 RCL 38
205 RCL 39
206 RCL 40
207 RCL 41
208 XROM
       "ZRIJKL"
209 ST+ 42
210 X<>Y
211 ST+ 43
212 X<>Y
213 RCL 43
214 RCL 42
215 CLD
216 RTN
217 GTO 13
218 LBL 12
219 "XEQ ZR
        FIRST"
220 AVIEW
221 END

( 407 bytes / SIZE 061+n(n+1)(n+2) )

STACK	INPUT	OUTPUT
Y	/	Im (R)
X	/	Re (R)

Example: The same one

XEQ "ZR" >>>> R = 0.014648 - 0.008939 i ---Execution time = 24m18s---

Notes:

-With a 2D-surface, the scalar curvature is twice the Gaussian curvature given by KGM
-"ZR" uses 4 subroutine-levels, provided the "GIJ"s do not call any subroutine.

>>> Always execute "ZR" before computing Einstein tensor or Weyl tensor. Registers R57-R58 must contain R

-Execution time becomes less than 4 seconds with V41.

-If n = 4, h =+0.1 and with relatively simple g_ij functions, all different from zero, V41 needs about 30 seconds to XEQ "ZR"

19°) Einstein Tensor

-Einstein tensor is a symmetric tensor defined as

E_ij = R_ij - (1/2) R g_ij

-Only after executing "ZR" - which stores R into R57-R58 - execute ZEIJ to get the components of Einstein tensor

>>> The listing is given in paragraph above

STACK	INPUTS	OUTPUTS
Y	i	Im E_ij
X	j	Re E_ij

Examples: The same metric at the same point

1 ENTER^
1 XEQ "ZEIJ" >>>> "XEQ ZR FIRST" and finally E₁₁ = 0.016161 - 0.006103 i ---Execution time = 5m12s---

-Likewise you will find

2 ENTER^
3 R/S >>>> E₂₃ = -0.127501 - 0.157186 i

-And 4 ENTER^ R/S >>>> E₄₄ = 0.149732 + 0.147500 i

Notes:

-When you XEQ "ZEIJ" the HP41 displays "XEQ ZR FIRST" to remind that R57-R58 must contain the scalar curvature R before executing "ZEIJ"
-This message is not displayed again if you press R/S to get the other components.

-"ZEIJ" uses 4 subroutine-levels, provided the "GIJ" do not call any subroutine.

-The cosmolical constant L is omitted here: the complete tensor is E_ij= R_ij - (1/2) R g_ij + L g_ij

-Elie Cartan proved that this tensor is the unique tensor T of order 2 involving derivatives of order < 3
and linear with respect to the 2nd derivatives that satisfies the conservative equations

D_i Tⁱ_j = 0 for all j where D_i is the covariant differentiation.

-They are the equations of General Relativity !

20°) Weyl Tensor

-Weyl tensor is the covariant tensor of order 4 defined as

W_ijkl = R_ijkl + [ 1/(n-2) ] ( R_il g_jk - R_ik g_jl + R_jk g_il - R_jl g_ik ) + [ R/(n-1)/(n-2) ] ( g_ik g_jl - g_il g_jk )

-Like Einstein tensor, execute "ZR" first before executing "ZWIJKL" so that data registers R57-R58 = R

>>> The listing is given 2 paragraphs above.

STACK	INPUTS	OUTPUTS
T	i	Im R_ijkl
Z	j	Re R_ijkl
Y	k	Im W_ijkl
X	l	Re W_ijkl

Examples: The same metric at the same point

1 ENTER^
2 ENTER^
1 ENTER^
2 XEQ "ZWIJKL" >>>> "XEQ ZR FIRST" and finally W₁₂₁₂ = -0.867899 + 1.483204 i ---Execution time = 10m05s---

And you also get the curvature tensor in registers Z & T: R₁₂₁₂ = -2.759976 + 4.320320 i

-Likewise

2 ENTER^
3 ENTER^
1 ENTER^
3 R/S >>>> W₂₃₁₃ = 0.061810 + 0.096108 i

And in Z&T-registers: R₂₃₁₃ = -2.872549 + 4.156863 i

Notes:

-In fact, all the components of Weyl tensor = 0 if n < 4
-The number of independent components is n(n-3)(n+1)(n+2)/12

-When you XEQ "ZWIJKL" the HP41 displays "XEQ ZR FIRST" to remind that R57-58 must contain the scalar curvature R before executing "ZWIJKL"
-This message is not displayed again if you press R/S to get another component.

>>> When the program stops, registers Z&T contains R_ijkl

-"ZWIJKL" uses 4 subroutine-levels, provided the "GIJ" do not call any subroutine.

-If n = 4, h = +0.1 and with relatively simple g_ij functions, all different from zero, V41 needs about 12 seconds to XEQ "ZWIJKL"

21°) Covariant Derivative of a Vector Field

-The usual partial derivative of a vector or a tensor is not a tensor, but the covariant derivatives are.
-They involve the partial derivatives and the Christoffel symbols.
-The formulae to get a tensor are not the same if the vector is covariant or contravariant:

Covariant Vector:

D_k V_i = ¶_k V_i - G^m_ik V_m

Contravariant Vector:

D_k Vⁱ = ¶_k Vⁱ + Gⁱ_km V^m

-So, in both cases, we get a tensor of order 2 whose coefficients are calculated and stored by "ZDCOV" & "ZDCOV^"

Data Registers: R00: temp

• R01 thru Rnn coordinates of the point ( n < 5 ) R61 thru R60+n(n+1)(n+2) = g_ij , g^ij , G^k_ij

R09 = n R11 to R17 are used by "ZFD"
R10 = h

Flag: F01
Subroutines: ZCIJK & "ZFD"

-This listing combines "ZCURL" "ZDIV" "ZLAPG" "DCOV" "DCOV^" "V-V^" "V^-V"

01 LBL "ZCURL"
02 XEQ 14
03 X<>Y
04 STO 21
05 ST+ X
06 E5
07 /
08 +
09 STO 26
10 CLX
11 STO 25
12 LBL 06
13 RCL 21
14 STO 22
15 LBL 07
16 RCL 21
17 RCL 22
18 X=Y?
19 CLST
20 X=Y?
21 GTO 00
22 RCL 18
23 1
24 -
25 RCL 22
26 +
27 RCL IND X
28 STO 00
29 RCL 10
30 RCL 21
31 XROM "ZFD"
32 STO 23
33 X<>Y
34 STO 24
35 RCL 18
36 1
37 -
38 RCL 21
39 +
40 RCL IND X
41 STO 00
42 RCL 10
43 RCL 22
44 XROM "ZFD"
45 RCL 24
46 RCL 23
47 Z-Z
48 LBL 00
49 X<>Y
50 STO IND 19
51 CHS
52 STO IND 26
53 DSE 19
54 RCL 26
55 1
56 -
57 R^
58 STO IND 19
59 CHS
60 STO IND Y
61 DSE 19
62 DSE 26

63 DSE 22
64 GTO 07
65 2
66 ST+ 25
67 RCL 25
68 ST- 19
69 RCL 26
70 FRC
71 RCL 19
72 +
73 STO 26
74 DSE 21
75 GTO 06
76 RCL 20
77 E3
78 /
79 RCL 09
80 ENTER
81 X^2
82 +
83 RCL 09
84 2
85 +
86 *
87 61
88 +
89 +
90 RCL 18
91 SIGN
92 RDN
93 RTN
94 LBL "ZDIV"
95 STO 18
96 INT
97 RCL 09
98 STO 19
99 +
100 STO 25
101 CLX
102 STO 21
103 STO 22
104 DSE 25
105 LBL 04
106 RCL 09
107 STO 20
108 CLX
109 STO 23
110 STO 24
111 LBL 05
112 RCL 19
113 RCL 20
114 ENTER
115 ZCIJK
116 ST+ 23
117 X<>Y
118 ST+ 24
119 DSE 20
120 GTO 05
121 RCL IND 25
122 STO 00
123 XEQ IND X
124 RCL 24

125 RCL 23
126 Z*Z
127 STO 23
128 X<>Y
129 STO 24
130 RCL 10
131 RCL 19
132 XROM "ZFD"
133 RCL 24
134 RCL 23
135 Z+Z
136 ST+ 21
137 X<>Y
138 ST+ 22
139 DSE 25
140 DSE 19
141 GTO 04
142 RCL 18
143 SIGN
144 RCL 22
145 RCL 21
146 RTN
147 LBL "ZLAPG"
148 CLX
149 STO 26
150 STO 27
151 RCL 09
152 ENTER
153 STO 23
154 ENTER
155 STO 24
156 ST+ 24
157 3
158 +
159 *
160 2
161 +
162 *
163 61
164 +
165 STO 30
166 ST+ 24
167 DSE 24
168 LBL 08
169 RCL 10
170 RCL 23
171 XROM "ZFD"
172 X<>Y
173 STO IND 24
174 DSE 24
175 X<>Y
176 STO IND 24
177 DSE 24
178 DSE 23
179 GTO 08
180 RCL 09
181 STO 22
182 ST+ X
183 ST+ 24
184 LBL 01
185 RCL 22
186 STO 23

187 LBL 02
188 CLX
189 STO 28
190 STO 29
191 RCL 09
192 STO 25
193 LBL 03
194 RCL 22
195 RCL 23
196 RCL 25
197 ZCIJK
198 Z=0?
199 DSE 24
200 Z=0?
201 GTO 00
202 RCL IND 24
203 DSE 24
204 RCL IND 24
205 Z*Z
206 ST+ 28
207 X<>Y
208 ST+ 29
209 LBL 00
210 DSE 24
211 DSE 25
212 GTO 03
213 RCL 09
214 ST+ X
215 ST+ 24
216 RCL 10
217 RCL 22
218 RCL 23
219 XROM "ZSD"
220 ST- 28
221 X<>Y
222 ST- 29
223 RCL 22
224 RCL 23
225 X=Y?
226 GTO 00
227 2
228 ST* 28
229 ST* 29
230 RDN
231 LBL 00
232 ENTER
233 CLX
234 ZCIJK
235 RCL 29
236 RCL 28
237 Z*Z
238 ST- 26
239 X<>Y
240 ST- 27
241 DSE 23
242 GTO 02
243 DSE 22
244 GTO 01
245 RCL 24
246 E3
247 /
248 RCL 30

249 +
250 RCL 27
251 RCL 26
252 RTN
253 LBL "ZDCOV"
254 CF 01
255 GTO 00
256 LBL"ZDCOV^"
257 SF 01
258 LBL 00
259 XEQ 14
260 RCL 09
261 STO 22
262 LBL 11
263 RCL 09
264 STO 21
265 LBL 12
266 RCL 09
267 STO 23
268 CLX
269 STO 24
270 STO 25
271 LBL 13
272 RCL 21
273 RCL 22
274 RCL 23
275 FS? 01
276 X<>Y
277 ZCIJK
278 Z=0?
279 GTO 00
280 STO 26
281 X<>Y
282 STO 27
283 RCL 18
284 1
285 -
286 RCL 23
287 +
288 RCL IND X
289 XEQ IND X
290 RCL 27
291 FS? 01
292 CHS
293 RCL 26
294 FS? 01
295 CHS
296 Z*Z
297 ST- 24
298 X<>Y
299 ST- 25
300 LBL 00
301 DSE 23
302 GTO 13
303 RCL 18
304 1
305 -
306 RCL 22
307 +
308 RCL IND X
309 STO 00
310 RCL 10

311 RCL 21
312 XROM "ZFD"
313 RCL 25
314 RCL 24
315 Z+Z
316 X<>Y
317 STO IND 19
318 DSE 19
319 X<>Y
320 STO IND 19
321 DSE 19
322 DSE 21
323 GTO 12
324 DSE 22
325 GTO 11
326 RCL 20
327 E3
328 /
329 RCL 19
330 +
331 RCL 18
332 SIGN
333 ST+ Y
334 RDN
335 RTN
336 LBL "ZV-V^"
337 SF 01
338 GTO 00
339 LBL "ZV^-V"
340 CF 01
341 LBL 00
342 XEQ 14
343 RCL 09
344 X^2
345 LASTX
346 -
347 ST+ X
348 -
349 RCL 18
350 INT
351 RCL 09
352 STO 21
353 4
354 *
355 +
356 DSE X
357 X<Y?
358 X<>Y
359 STO 19
360 STO 20
361 LBL 09
362 CLX
363 STO 23
364 STO 24
365 RCL 09
366 STO 22
367 LBL 10
368 RCL 21
369 RCL 22
370 1
371 CHS
372 FS? 01

373 CLX
374 ZCIJK
375 RCL 18
376 INT
377 RCL 22
378 ST+ X
379 +
380 DSE X
381 STO T
382 RDN
383 RCL IND Z
384 DSE T
385 RCL IND T
386 Z*Z
387 ST+ 23
388 X<>Y
389 ST+ 24
390 DSE 22
391 GTO 10
392 RCL 24
393 STO IND 20
394 DSE 20
395 RCL 23
396 STO IND 20
397 DSE 20
398 DSE 21
399 GTO 09
400 RCL 19
401 E3
402 /
403 RCL 19
404 RCL 09
405 ST+ X
406 -
407 +
408 RCL 18
409 SIGN
410 ST+ Y
411 RDN
412 RTN
413 LBL 14
414 STO 18
415 RCL 09
416 ENTER
417 STO Z
418 5
419 +
420 *
421 2
422 +
423 *
424 60
425 +
426 STO 19
427 STO 20
428 RTN
429 END

( 744 bytes / SIZE var )

STACK	INPUT	OUTPUT
X	bbb.eee(V)	bbb.eee(DV)
L	/	bbb.eee(V)

Examples: With the same metric as above and still at the same point x = y = z = t = 1 + 2 i

Covariant Vector Field

V₁ = x² + y z t V₂ = x² y² t³ V₃ = x y z t V₄ = x² y² z² t²

-Load the corresponding routines in main memory ( the names are arbitrary, provided they have less than 7 characters - global LBLs )

01 LBL "V1"
02 RCL 08
03 RCL 07
04 RCL 06
05 RCL 05
06 Z*Z
07 RCL 04
08 RCL 03
09 Z*Z
10 RCL 02
11 RCL 01

12 Z^2
13 Z+Z
14 RTN
15 LBL "V2"
16 RCL 02
17 RCL 01
18 RCL 04
19 RCL 03
20 Z*Z
21 RCL 08
22 RCL 07

23 Z*Z
24 Z^2
25 RCL 08
26 RCL 07
27 Z*Z
28 RTN
29 LBL "V3"
30 RCL 02
31 RCL 01
32 RCL 04
33 RCL 03

34 Z*Z
35 RCL 06
36 RCL 05
37 Z*Z
38 RCL 08
39 RCL 07
40 Z*Z
41 RTN
42 LBL "V4"
43 RCL 02
44 RCL 01

45 RCL 04
46 RCL 03
47 Z*Z
48 RCL 06
49 RCL 05
50 Z*Z
51 RCL 08
52 RCL 07
53 Z*Z
54 Z^2
55 RTN

-Store the names of theses subroutines in 4 consecutive registers, without disturbing R61 to R180 , or R09 , R11 to R27 which are used by "DCOV"

-For instance V1 ASTO 31 ...... V4 ASTO 34
-Control number 31.034

SIZE 213 at least

-Still with h = -0.01 in R10

31.034 XEQ "ZDCOV" >>>> 181.212 ---Execution time = 2m49s---

-And we obtain the coefficients of the matrix [ D_i V_j ] , in column order, rounded to 6 decimals:

     122.576240 + 35.714716 i       156.135392 +   2.146502 i                    -11 - 2 i                   30.384990 + 277.137180 i
     -80.864608 - 81.853498 i        180.805784 + 132.422126 i   -119.686275 - 40.254902 i                 58 + 556 i
                  - 3 + 4 i                     -108.686275 - 38.254902 i      120.058824 + 39.431373 i                58 + 556 i
     -30.615010 - 274.862820 i                  351   + 132 i                             -11 - 2 i                  -24.025561 + 321.947514 i

-For example in R187-R188, D₃ V₁ = -3 + 4 i

Contravariant Vector Field

-If the same functions are the components of a contravariant vector field: V¹ = x² + y z t V² = x² y² t³ V³ = x y z t V⁴ = x² y² z² t²

31.034 XEQ "ZDCOV^" >>>> 181.212 ---Execution time = 2m49s---

-And we obtain the coefficients of the matrix [ D_i V^j ] , in column order, rounded to 6 decimals:

      77.335435 + 94.305170 i       355.656863 + 121.705882 i                 -11 - 2 i                   94.818411 + 858.464061 i
   -122.135203 - 34.150438 i       221.029412 + 92.549020 i    -142.058824 - 43.431373 i     48.800484 + 532.449814 i
                  - 3 + 4 i                    131.058824 + 41.431373 i       97.686275 + 36.254902 i                58 + 556 i
     -31.923111 - 271.977506 i                  351   + 132 i                             -11 - 2 i               136.006271 + 802.014928 i

-For example in R187-R188, D₃ V¹ = -3 + 4 i

22°) Curl of a Covariant Vector Field

-The partial derivatives of a vector is not a tensor and Curl_ij V is defined as D_i V_j - D_j V_i
-But if we apply the formulae given in the paragraph above, the Christoffel symbols vanish and finally it yields

Curl_ij V = ¶_i V_j - ¶_j V_i

-This tensor is obviously antisymmetric, so Curl_ii V = 0 and Curl_ji V = - Curl_ij V

-"ZCURL" calculates and stores all the coefficients of the matrix [ Curl_ij V ] , in column order

Data Registers: R00: temp

• R01 thru Rnn coordinates of the point ( n < 5 ) R61 thru R60+n(n+1)(n+2) = g_ij , g^ij , G^k_ij

                                  • R09 = n      R20 to R26 are used by "ZCURL"
                                  • R10 = h
Flags: /
Subroutine:         "ZFD"

>>> The listing is given in the paragraph "Covariant derivative of a vector field" above

STACK	INPUTS	OUTPUTS
X	bbb.eee	bbb.eee_CURL
L	/	bbb.eee

Example: Again the same metric at the same point and the same covariant vector field:

V₁ = x² + y z t V₂ = x² y² t³ V₃ = x y z t V₄ = x² y² z² t²

-If the function names are still in R31 R32 R33 R34 -> control number = 31.034

31.034 XEQ "ZCURL" >>>> 181.212 ---Execution time = 1m23s---

-And you get

           0               237 + 84 i       - 8 - 6 i           61 + 552 i
   -237 - 84 i               0             -11 - 2 i        -293 + 424 i
       8 + 6 i           11 + 2 i              0                69 + 558 i          For example,       Curl₁₃ V = ¶₁ V₃ - ¶₃ V₁ = - 8 - 6 i
    -61 - 552 i      293 - 424 i     -69 - 558 i              0

Notes

-This definition will not give the same results as "ZCDGL" in cylindrical and spherical coordinates.
because the usual formulae in 3D-Euclidean space involve what is called "Vraies Grandeurs" in French ( "True Magnitudes" in English ? )
which are neither the covariant nor the contravariant components, except with an orthonormal basis
-In tensor calculus, this "vraie grandeur" is the product of the contravariant component V^k by the modulus of the corresponding vector e_k of the basis
or the quotient of the covariant component V_k by the same quantity:

V^k e_k = V_k / e_k = (g_kk)^1/2 V^k = (g_kk) ^-1/2 V_k ( No summation here )

23°) Divergence of a Contravariant Vector Field

-The divergence of a vector field V is given by

Div V = D_i Vⁱ = ¶_i Vⁱ + Gⁱ_ik V^k

-It is a scalar

Data Registers: R00: temp

• R01 thru R2n coordinates of the point ( n < 5 ) R61 thru R60+n(n+1)(n+2) = g_ij , g^ij , G^k_ij

• R09 = n R19 to R25 are used by "ZDIV"
• R10 = h

Flags: /
Subroutines: ZCIJK "ZFD"

>>> The listing is given in the paragraph "Covariant derivative of a vector field" above

STACK	INPUTS	OUTPUTS
Y	/	Im Div V
X	bbb.eee	Re Div V
L	/	bbb.eee

Example: Again the same metric at the same point and the same contravariant vector field:

V¹ = x² + y z t V² = x² y² t³ V³ = x y z t V⁴ = x² y² z² t²

-> control number = 31.034

31.034 XEQ "ZDIV" >>>> 532.057392 ---Execution time = 37s---
X<>Y 1025.124020

-Thus, Div V = 532.057392 + 1025.124020 i

Note:

-For the same reason mentioned above, the divergence here will be generally different from the usual definition in Euclidean calculus.

24°) Laplacian and Gradient of a Scalar Field

-The gradient of a scalar field f is defined as Grad_k f = ¶_k f i-e simply the usual partial derivatives: it is a covariant vector

-The Laplacian of the same scalar field f is a scalar defined by

Lapl(f) = g^jk ( ¶_jk f - Gⁱ_jk ¶_i f )

Data Registers: • R00 = f name ( Register R00 is to be initialized before executing "ZLAPG" )

• R01 thru R2n coordinates of the point ( n < 5 ) R61 thru R60+n(n+1)(n+2) = g_ij , g^ij , G^k_ij

• R09 = n R22 to R30 are used by "ZLAPG"
• R10 = h R11 thru R21 are used by "ZSD"

Flags: /
Subroutines: ZCIJK "ZFD" "ZSD" and a function that computes f assuming x¹ is in R01-R02 .................. xⁿ is in R2n-1R2n

>>> The listing is given in the paragraph "Covariant derivative of a vector field" above

STACK	INPUTS	OUTPUTS
Z	/	bbb.eee_GRAD
Y	/	Im Lapl (f)
X	/	Re Lapl (f)

With bbb = 61+n(n+1)(n+2)

Example: Again the same metric at the same point and the following scalar field

f(x,y,z,t ) = x² + y z t i-e V¹ in the paragraph above

V1 ASTO 00

    XEQ "ZLAPG" >>>>   -0.127458                                                                        ---Execution time = 3m12s---
                             RDN      0.179583
                             RDN    181.188

-So, Lap(f) = -0.127458 + 0.179583 i

And you get the covariant components of the gradient in R181 to R188:

Grad(f) = [ 2+4i , -3+4i , -3+4i , -3+4i ]

Notes:

-For the same reason mentioned above, the gradient here will be generally different from the usual definition in Euclidean calculus.
-Contrariwise, the 2 definitions of the Laplacian match because f & Lapl(f) are scalars.

-If you prefer the contravariant components of the gradient, use the folowing routine "ZV-V^"

25°) Covariant Vector <> Contravariant Vector

-"V-V^" transforms the covariant components of a vector into its contravariant components
-"V^-V" performs the reverse operation.

Data Registers: R00 R12-R14: temp

• R01 thru Rnn coordinates of the point ( n < 8 )

• R09 = n R20 to R26 are used by "V<>V^" R61 thru R60+n(n+1)(n+2) = g_ij , g^ij , G^k_ij
• R10 = h

Flag: F01                     SF01 = V>V^
                                      CF01 = V^>V
Subroutine:    ZCIJK

>>> The listing is given in the paragraph "Covariant derivative of a vector field" above

STACK	INPUTS	OUTPUTS
X	bbb.eee(V)	bbb.eee(V')
L	/	bbb.eee

Example: Again the same metric at the same point.

-Suppose the covariant components of a vector is [ 2+3i , 3+4i , 4+5i , 6+7i ]
-Let's calculate the contravariant components.

-If you store these 4 complex numbers into R41 thru R48

41.048 XEQ "ZV-V^" >>>> 181.188 ---Execution time = 20s---

-And you get the contravariant components in R181 to R188

[ -0.204560 +0.009713 i , -0.186275 + 0.078431 i , -0.235294 + 0.107843 i , -0.360928 + 0.135817 i ] ( rounded to 6D )

-To recalculate the covariant components

181.188 XEQ "ZV^-V^" >>>> 189.196 and the covariant components [ 2+3i , 3+4i , 4+5i , 6+7i ] are in R189 to 196

Note:

-The HP41 stores the results in a block of consecutive registers starting at RBB
where BB = Max ( bb+2n , 61+n(n+1)(n+2) )

26°) Eigenvalues & Eigenvector of a Complex Matrix

-"ZEGVL" & "ZEGVH" employ a variant of the Jacobi's algorithm:

*** The eigenvalues of A are the diagonal-elements of an upper triangular matrix T equal to the infinite product ...U_k^-1.U_k-1^-1....U₂^-1.U₁^-1.A.U₁.U₂....U_k-1.U_k....
where the U_k are unitary matrices. The eigenvectors are the columns of U = U₁.U₂....U_k-1.U_k.... if A is Hermitian ( i-e if A equals its transconjugate )
Actually, T is diagonal if A is Hermitian.

-First, the HP41 finds the greatest element a_{i j} below the main diagonal
-Then, U₁ is determined so that it places a zero in position ( i , j ) in U₁^-1.A.U₁

U₁ has the same elements as the Identity matrix, except that u_{i i} = u_{j j} = x and u_{i j} = y + i.z , u_{j i} = -y + i.z

with y + i.z = C.x , x = ( 1 + | C |² )^-1/2 and C = 2.a_{i j} / [ a_{i i} - a_{i j} + ( ( a_{i i} - a_{i j} )² + 4 a_{i j} a_{j i} )^1/2 ]

-The process is repeated until the greatest rounded sub-diagonal element is zero.

-The successive greatest | a_{i j} | ( i > j ) are displayed when the routine is running.

*** If the matrix is non-Hermitian, and if all the eigenvalues are distinct, we define an upper triangular matrix W = [ w_{i j} ] by

   w_{i j} = 0    if i > j
   w_{i i} = 1
   w_{i j} = Sum_{k = i+1, i+2 , .... , j}      t_{i k} w_{k j} / ( t_{j j} - t_{i i} )      if i < j

-Then, the n columns of U.W are the eigenvectors corresponding to the n eigenvalues t₁₁ = k₁ , ............... , t_nn = k_n

Data Registers: • R00 = n ( All the • Registers are to be initialized before executing "ZEGVL" or "ZEGVH" )

                            •   R01,R02 = a₁₁     .................      • R2n²-2n+1,R2n²-2n+2 = a_1n
                              .................................................................................................                        ( the n² elements of A in column order )
   INPUTS
                            •   R2n-1,R2n = a_n1 .................      • R2n²-1,R2n² = a_nn

--------- odd registers = real part of the coefficients , even registers = imaginary part of the coefficients

   DURING            R01 thru R2n² = the "infinite" product   ...U_k^-1.U_k-1^-1....U₂^-1.U₁^-1.A.U₁.U₂....U_k-1.U_k....
       THE                R2n²+1 thru R4n² = U₁.U₂........U_k.....
ITERATIONS      R4n²+1 thru R4n²+2n: temp ( for non-Hermitian matrices only )

---------

                             R00 = n
                             R01,R02 = k₁          R2n+1,R2n+2 = V₁⁽¹⁾      ....................    R2n²+1,R2n²+2 = V₁⁽ⁿ⁾
                             .......................            ........................                ....................     ..................................
OUTPUTS
                             R2n-1,R2n = k_n        R4n-1,R4n =    V_n⁽¹⁾      .....................   R2n²+2n-1,R2n²+2n = V_n⁽ⁿ⁾

( n eigenvalues ; 1st eigenvector ; .................................. ; n^th eigenvector )

Flag: F02 SF02 for a Hermitian matrix
CF02 for a non-Hermitian matrix
Subroutines: /

01 LBL "ZEGVL"
02 CF 02
03 GTO 14
04 LBL "ZEGVH"
05 SF 02
06 LBL 14
07 RCL 00
08 X^2
09 ST+ X
10 .1
11 %
12 ST+ X
13 +
14 ISG X
15 CLRGX
16 RCL 00
17 1
18 +
19 ST+ X
20 E5
21 /
22 +
23 SIGN
24 LBL 00
25 STO IND L
26 ISG L
27 GTO 00
28 LBL 01
29 RCL 00
30 ST+ X
31 E3
32 /
33 ISG X
34 STO M
35 2
36 +
37 ENTER
38 ENTER
39 CLX
40 LBL 02
41 RCL IND Z
42 X^2
43 SIGN
44 ST+ T
45 CLX
46 RCL IND T
47 ST* X
48 ST+ L
49 X<> L
50 X>Y?
51 STO Z
52 X>Y?
53 +
54 RDN
55 ISG Z
56 GTO 02
57 2
58 ST+ M
59 CLX
60 RCL M
61 ST+ T
62 RDN
63 ISG Z
64 GTO 02
65 SQRT
66 VIEW X
67 ENTER
68 RND
69 X=0?
70 GTO 06
71 X<> Z

72 INT
73 STO M
74 2
75 ST- Y
76 /
77 STO Y
78 RCL 00
79 ST/ Z
80 MOD
81 ST+ X
82 RCL X
83 LASTX
84 *
85 +
86 X<>Y
87 INT
88 ST+ X
89 RCL X
90 RCL 00
91 *
92 +
93 2
94 ST+ Z
95 +
96 STO N
97 X<>Y
98 STO O
99 +
100 RCL M
101 -
102 RCL IND X
103 DSE Y
104 RCL IND Y
105 4
106 ST* Z
107 *
108 RCL IND M
109 RCL M
110 DSE X
111 RDN
112 RCL IND T
113 Z*Z
114 RCL IND N
115 RCL IND O
116 -
117 DSE N
118 DSE O
119 RCL IND N
120 STO N
121 CLX
122 RCL IND O
123 ST- N
124 RDN
125 STO O
126 RCL N
127 Z^2
128 Z+Z
129 SQRTZ
130 RCL O
131 ST+ Z
132 X<> N
133 +
134 RCL IND M
135 ST+ X
136 DSE M
137 RCL IND M
138 ST+ X
139 R^
140 R^
141 Z=0?
142 FS? 30

143 Z/Z
144 RCL Y
145 X^2
146 RCL Y
147 X^2
148 +
149 ENTER
150 SIGN
151 ST- M
152 +
153 SQRT
154 1/X
155 STO O
156 ST* Z
157 *
158 STO N
159 X<>Y
160 X<> M
161 2
162 /
163 STO Y
164 RCL 00
165 ST/ Z
166 MOD
167 X<>Y
168 INT
169 RCL 00
170 ST+ X
171 ST* Y
172 ST* Z
173 +
174 .1
175 %
176 +
177 RCL 00
178 ST+ X
179 -
180 1
181 ST+ Z
182 +
183 RCL 00
184 X^2
185 ST+ X
186 ST+ Z
187 .1
188 %
189 +
190 +
191 XEQ 03
192 RCL 00
193 ST+ X
194 ST- Z
195 -
196 RCL 00
197 X^2
198 ST+ X
199 ST- Z
200 .1
201 %
202 +
203 -
204 XEQ 03
205 RCL 00
206 ST/ Z
207 /
208 INT
209 X<>Y
210 INT
211 RCL 00
212 X^2
213 ST+ X

214 E3
215 /
216 +
217 RCL 00
218 ST+ X
219 1
220 CHS
221 ST* M
222 ST* N
223 ST+ Z
224 ST+ T
225 +
226 E5
227 /
228 ST+ Z
229 +
230 XEQ 03
231 GTO 01
232 LBL 03
233 STO P
234 LBL 04
235 CLX
236 RCL IND P
237 RCL O
238 *
239 RCL IND Y
240 RCL N
241 *
242 +
243 PI
244 SIGN
245 ST+ T
246 CLX
247 RCL IND T
248 RCL M
249 *
250 -
251 STO Q
252 CLX
253 RCL IND Y
254 RCL O
255 *
256 RCL IND P
257 RCL N
258 *
259 -
260 PI
261 SIGN
262 ST+ P
263 CLX
264 RCL IND P
265 RCL M
266 *
267 -
268 X<> IND Y
269 RCL M
270 *
271 RCL IND P
272 RCL O
273 *
274 +
275 PI
276 SIGN
277 ST+ Z
278 CLX
279 RCL IND Z
280 RCL N
281 *
282 +
283 X<> IND P
284 RCL N

285 *
286 RCL IND Y
287 RCL O
288 *
289 X<>Y
290 -
291 RCL P
292 SIGN
293 ST- L
294 X<> Q
295 X<> IND L
296 RCL M
297 *
298 +
299 STO IND Y
300 ISG Y
301 CLX
302 ISG P
303 GTO 04
304 X<> P
305 RTN
306 LBL 05
307 RCL P
308 SIGN
309 ST- L
310 RCL IND P
311 RCL IND L
312 RCL IND 04
313 DSE 04
314 RCL IND 04
315 Z*Z
316 RTN
317 LBL 06
318 CLD
319 FS? 02
320 GTO 11
321 RCL 00
322 X^2
323 ST+ X
324 LASTX
325 1
326 +
327 ST+ X
328 E5
329 /
330 +
331 STO O
332 RCL 00
333 DSE X
334 E3
335 /
336 ISG X
337 STO M
338 LBL 07
339 RCL 00
340 ST+ X
341 ENTER
342 X^2
343 +
344 STO Y
345 RCL M
346 INT
347 ST+ X
348 STO 03
349 -
350 E3
351 /
352 +
353 STO 04
354 RCL O
355 RCL 03

356 -
357 2 E-5
358 -
359 STO P
360 CLX
361 STO 03
362 STO N
363 STO IND Y
364 SIGN
365 ST- Y
366 STO IND Y
367 LBL 08
368 XEQ 05
369 ST+ 03
370 X<>Y
371 ST+ N
372 DSE P
373 DSE 04
374 GTO 08
375 RCL IND O
376 RCL IND P
377 -
378 PI
379 SIGN
380 ST- 04
381 ST- O
382 ST- P
383 CLX
384 RCL IND O
385 RCL IND P
386 -
387 RCL N
388 RCL 03
389 R^
390 R^
391 Z/Z
392 STO IND 04
393 CLX
394 SIGN
395 ST+ 04
396 ST+ O
397 X<>Y
398 STO IND 04
399 ISG M
400 GTO 07
401 RCL 00
402 STO 03
403 DSE M
404 LBL 09
405 RCL 00
406 X^2
407 ENTER
408 STO 04
409 ST+ 04
410 LASTX
411 ST+ 04
412 RCL M
413 INT
414 ST- 04
415 *
416 +
417 STO N
418 E3
419 /
420 +
421 RCL 03
422 ST+ N
423 +
424 RCL 00
425 E5
426 /
427 +

428 2
429 ST* 04
430 ST* N
431 *
432 STO P
433 LBL 10
434 XEQ 05
435 RCL N
436 DSE X
437 RDN
438 ST+ IND T
439 X<>Y
440 ST+ IND N
441 PI
442 INT
443 ST+ 04
444 ISG P
445 GTO 10
446 DSE 03
447 GTO 09
448 RCL M
449 FRC
450 E-3
451 -
452 STO M
453 DSE O
454 ISG M
455 GTO 07
456 LBL 11
457 RCL 00
458 X^2
459 ST+ X
460 STO M
461 LASTX
462 ST+ X
463 STO Z
464 1
465 +
466 E5
467 /
468 +
469 LBL 12
470 RCL IND X
471 STO IND Z
472 CLX
473 SIGN
474 ST- Z
475 -
476 RCL IND X
477 STO IND Z
478 DSE Y
479 RDN
480 DSE Y
481 GTO 12
482 CLX
483 RCL M
484 R^
485 1
486 ST+ Z
487 +
488 E3
489 /
490 +
491 RCL M
492 E6
493 /
494 +
495 REGMOVE
496 RCL 02
497 RCL 01
498 CLA
499 END

( 780 bytes / SIZE 4n²+1 if A is Hermitian / SIZE 4n²+2n+1 if A is non-Hermitian )

STACK	INPUTS	OUTPUTS
Y	/	y₁
X	/	x₁

where k₁ = x₁ + y₁ = the first eigenvalue.

Example1:

             1      4 -7.i   3-4.i
A = 4+7.i      6      1-5.i             A is equal to its transconjugate so A is Hermitian
          3+4.i   1+5.i      7

              1 0   4 -7 3 -4           R01 R02    R07 R08    R13 R14
Store   4 7   6 0   1 -5   into   R03 R04    R09 R10    R15 R16    respectively     and n = 3 STO 00
             3 4   1 5   7 0            R05 R06    R11 R12    R17 R18

-The matrix is Hermitian, so for instance FIX 8 XEQ "ZEGVH" yields ---Execution time = 4mn21s---

           k₁ = 15.61385274 + 0.i = ( X , Y ) = ( R01 , R02 )
           k₂ = 5.230678473 + 0.i = ( R03 , R04 )
           k₃ = -6.844531166 + 0.i = ( R05 , R06 )

and the 3 corresponding eigenvectors:

    V₁( 0.481585120 + 0.108201123 i , 0.393148652 + 0.527305194 i , -0.142958847 + 0.550739893 i )    = ( R07 R08 , R09 R10 , R11 R12 )
    V₂( -0.436763638 + 0.075843695 i , 0.210271354 - 0.463555613 i , -0.516799072 + 0.526598642 i )   = ( R13 R14 , R15 R16 , R17 R18 )
    V₃( -0.706095475 + 0.247553486 i , 0.326530385 + 0.449069516 i , 0.363126603 - 0. i )                      = ( R19 R20 , R21 R22 , R23 R24 )

-Due to roundoff errors, registers R02 R04 R06 contain very small numbers instead of 0
but actually, the eigenvalues of a Hermitian matrix are always real.

Example2:

           1+2.i   2+5.i   4+7.i
   A = 4+7.i   3+6.i   3+4.i        A is non-Hermitian
           3+4.i   1+7.i   2+4.i

              1 2   2 5   4 7           R01 R02    R07 R08    R13 R14
   Store   4 7   3 6   3 4   into   R03 R04    R09 R10    R15 R16    respectively     and n = 3 STO 00
              3 4   1 7   2 4            R05 R06    R11 R12    R17 R18

-The matrix is non-Hermitian, so FIX 8 XEQ "ZEGVL" produces ---Execution time = 5mn20s---

           k₁ = 7.656606017 + 15.61073839 i = ( X , Y ) = ( R01 , R02 )
           k₂ = 1.661248139 - 1.507335318 i   = ( R03 , R04 )
           k₃ = -3.317854157 - 2.103403077 i = ( R05 , R06 )   and the 3 corresponding eigenvectors

    V₁( 0.523491429 + 0.008555748 i , 0.650242685 - 0.046353000 i , 0.546288478 + 0.049882590 i )      = ( R07 R08 , R09 R10 , R11 R12 )
    V₂( -0.4777850326 - 0.196368365 i , 0.660547554 - 0.265888146 i , -0.356842635 + 0.318267519 i ) = ( R13 R14 , R15 R16 , R17 R18 )
    V₃( -0.567221101 - 0.574734664 i , 0.141603480 + 0.549379028 i , 0.480533314 - 0.090337846 i )     = ( R19 R20 , R21 R22 , R23 R24 )

Notes:

-The precision is controlled by the display format
-If the matrix is non-hermitian, "ZEGVL" will fail if 2 or more eigenvalues are equal.

27°) Complex Polar <> Rectangular conversion

-When x y and r u are complexes, the formulae:

x = r cos u
y = r sin u

may be solved by r = sqrt ( x² + y² ) and u = - i Ln [ (x + i y ) / r ]

-This assumes that r # 0 unless x = y = 0
-If r = 0 and x # 0 or y # 0, the HP41 will display NONEXISTENT

Data Registers: R00: unused ( Registers R01 thru R04 are to be initialized before executing "ZP-R" or "ZR-P" )

INPUT • R01-R02 = r • R03-R04 = u or • R01-R02 = x • R03-R04 = y

OUTPUT R01-R02 = x , R03-R04 = y or R01-R02 = r , R03-R04 = u

Flags: /
Subroutines: "SINZ" "COSZ" Z+Z LNZ Z^2

01 LBL "ZP-R"
02 RCL 04
03 RCL 03
04 XROM "SINZ"
05 RCL 02
06 RCL 01
07 Z*Z
08 X<> 03
09 X<>Y

10 X<> 04
11 X<>Y
12 XROM "COSZ"
13 RCL 02
14 RCL 01
15 Z*Z
16 STO 01
17 X<>Y
18 STO 02

19 X<>Y
20 RCL 04
21 RCL 03
22 R^
23 R^
24 RTN
25 LBL "ZR-P"
26 RCL 04
27 RCL 03

28 Z^2
29 RCL 02
30 RCL 01
31 Z^2
32 Z+Z
33 SQRTZ
34 X<> 01
35 X<>Y
36 X<> 02

37 X<>Y
38 RCL 03
39 RCL 04
40 CHS
41 Z+Z
42 RCL 02
43 RCL 01
44 Z=0?
45 FS? 30

46 GTO 00
47 RCL 03
48 RCL 04
49 Z=0?
50 GTO 01
51 SF 41
52 LBL 00
53 Z/Z
54 LNZ

55 CHS
56 LBL 01
57 STO 04
58 X<>Y
59 STO 03
60 RCL 02
61 RCL 01
62 END

( 99 bytes / SIZE 005 )

STACK	INPUT1	OUTPUT1	INPUT2	OUTPUT2
T	/	Im y	/	Im u
Z	/	Re y	/	Re u
Y	/	Im x	/	Im r
X	/	Re x	/	Re r

Example1: r = 1 + 2 i , u = 3 + 4 i

1 STO 01 3 STO 03
2 STO 02 4 STO 04

    XEQ "ZP-R" >>>>   -19.33263894 = R01
                          RDN    -57.92104456 = R02
                          RDN     57.88736456 = R03
                          RDN    -19.30933718 = R04

-Thus,

x = -19.33263894 - 57.92104456 i
y = 57.88736456 - 19.30933718 i

Example2: With the results obtained above:

    XEQ "ZR-P" or simply R/S >>>>   1 = R01        with small roundoff-errors in the last decimals
                                                   RDN   2 = R02
                                                   RDN   3 = R03
                                                   RDN   4 = R04

r = 1 + 2 i , u = 3 + 4 i

28°) Complex Spherical <> Rectangular conversion

"ZS-R" & "ZR-S" actuallly deal with ( complex ) longitudes & latitudes

       x = r cos b cos L
       y = r cos b sin L
       z = r sin b

where x , y , z = rectangular coordinates, r = ( x² + y² + z² )^1/2 , L = longitude , b = latitude

Data Registers: R00: unused ( Registers R01 thru R04 are to be initialized before executing "ZP-R" or "ZR-P" )

INPUT • R01-R02 = r • R03-R04 = L • R05-R06 = b or • R01-R02 = x • R03-R04 = y • R05-R06 = z

OUTPUT R01-R02 = x R03-R04 = y R05-R06 = z or R01-R02 = r R03-R04 = L R05-R06 = b

Flags: /
Subroutines: "ZP-R" "ZR-P"

01 LBL "ZS-R"
02 XROM "ZP-R"
03 RCL 03
04 X<> 05
05 STO 03
06 RCL 04
07 X<> 06
08 STO 04
09 XROM "ZP-R"
10 RTN
11 LBL "ZR-S"
12 XROM "ZR-P"
13 RCL 03
14 X<> 05
15 STO 03
16 RCL 04
17 X<> 06
18 STO 04
19 XROM "ZR-P"
20 END

( 44 bytes / SIZE 007 )

STACK	INPUT1	OUTPUT1	INPUT2	OUTPUT2
T	/	Im y	/	Im L
Z	/	Re y	/	Re L
Y	/	Im x	/	Im r
X	/	Re x	/	Re r

Example1: r = 5 + 7 i , L = 1 - 0.04 i , b = 1.2 - 0.07 i

5 STO 01 1 STO 03 1.2 STO 05
7 STO 02 -0.04 STO 04 -0.07 STO 06

    XEQ "ZS-R" >>>>   0.638343618 = R01
                          RDN    1.597236347 = R02
                          RDN    2.406270999 = R03
                          RDN    3.631178539 = R04

and R05 = 4.362046249 R06 = 5.786920482

So,

    x = 0.638343618 + 1.597236347 i
    y = 2.406270999 + 3.631178539 i
    z = 4.362046249 + 5.786920482 i

Example2: With the results obtained above:

    XEQ "ZR-S" or simply R/S >>>>     5 =     R01        with small roundoff-errors in the last decimals
                                                   RDN     7 =     R02
                                                   RDN     1 =     R03
                                                   RDN -0.04 = R04   and R05 = 1.2    R06 = -0.07

And we find again r = 5 + 7 i , L = 1 - 0.04 i , b = 1.2 - 0.07 i

Note:

-If you want the spherical coordinates ( r , u , v ) that are used in "ZCDGL" with SF02

u = PI/2 - b and v = L ( u = co-latitude )

29°) Complex Dot-Product & Hermitian Product

-The complex dot product of 2 vectors V( x₁ ,..........., x_n ) V'( y₁ ,........, y_n ) is defined as usual by

< V , V' > = x₁ y₁ + ................. + x_n y_n

-Whereas the Hermitian product is < V , V' >_H = x₁ y₁bar + ................. + x_n y_nbar where ybar = conjugate of y

( if y = 2 + 3 i ybar = 2 - 3 i )

"ZDOT" calculates < V , V' >
"ZDOTH" calculates < V , V' >_H

01 LBL "ZDOT"
02 CF 04
03 GTO 00
04 LBL "ZDOTH"
05 SF 04
06 LBL 00
07 CLA
08 STO N
09 X<>Y
10 STO M
11 LBL 01
12 RCL IND M
13 ISG M
14 RCL IND M
15 X<>Y
16 RCL IND N
17 ISG N
18 RCL IND N
19 FS? 04
20 CHS
21 X<>Y
22 Z*Z
23 ST+ O
24 X<>Y
25 ST+ P
26 ISG M
27 FS? 30
28 GTO 00
29 ISG N
30 GTO 01
31 LBL 00
32 RCL P
33 RCL O
34 CLA
35 END

( 74 bytes )

STACK	INPUTS	OUTPUTS
Y	bbb.eee(V)	y
X	bbb.eee(V')	x

With x + i y = < V , V' > or < V , V' >_H

Example: V(1+2i,3+4i,5+6i,7+9i) V'(2+5i,3+7i,8+6i,1+4i)

-Store for instance

1 2 3 4 5 6 7 9 into R11 R12 R13 R14 R15 R16 R17 R18
2 5 3 7 8 6 1 4 into R21 R22 R23 R24 R25 R26 R27 R28

11.018 ENTER^
21.028 XEQ "ZDOT" >>>> -52
X<>Y 157 < V , V' > = -52 + 157 i

-Likewise,

11.018 ENTER^
21.028 XEQ "ZDOTH" >>>> 168
X<>Y -11 < V , V' >_H = 168 - 11 i

Notes:

-"ZDOT" clears flag F04 whereas "ZDOTH" sets F04
-Synthetic register P is used, so don't interrupt these routines.
-The alpha register is cleared.
-These programs only use the registers where the complex coordinates are stored
-They stop when the shorter vector has been used:
in the example above, 11.018 ENTER^ 21.049 or 11.100 ENTER^ 21.028 would return the same result.

30°) E^Z-1

exp(x+i.y) - 1 is calculated by ( exp(x) - 1 ) cos y + i ( exp(x) - 1 ) sin y - 2 sin² y/2 + i sin y to get a petter precision for small arguments.

Data Registers: /
Flags: /
Subroutines: /

01 LBL "E^Z-1"
02 X<>Y
03 STO P
04 X<>Y
05 STO O
06 E^X
07 RAD
08 P-R
09 X<>Y
10 STO N
11 RDN
12 STO M
13 ABS
14 ENTER
15 SIGN
16 ST- M
17 E^X-1
18 X<Y?
19 GTO 00
20 CLX
21 RCL P
22 COS
23 X<>Y
24 X<> O
25 E^X-1
26 *
27 STO M
28 X<> P
29 2
30 /
31 SIN
32 X^2
33 ST+ X
34 ST- M
35 ENTER
36 LBL 00
37 RDN
38 X<> N
39 RCL M
40 DEG
41 CLA
42 END

( 66 bytes / SIZE 000 )

STACK	INPUTS	OUTPUTS
Y	y	Im ( e^z - 1 )
X	x	Re ( e^z - 1 )

With z = x + i y

Examples:

    3 ENTER^
    2 XEQ "E^Z-1"   >>>>    -8.315110095
                                X<>Y    1.042743656

-Thus, e^2+3.i - 1 = -8.315110095 + 1.042743656 i

   0.002 ENTER^
   0.001      R/S    >>>>    0.0009984981667
                           X<>Y    0.002001999666

-So, e^{0.001 + 0.002.i} - 1 = 0.0009984981667 + 0.002001999666 i

PI 2 / 61 CHS R/S >>>> -1
X<>Y 3.221340286 E-27

-And e^-61+PI/2 = -1 + 3.221340286 E-27 i

31°) Complex Einstein's Functions

4 functions are often called "Einstein Functions"

E₁(x) = x² exp(x) / ( exp(x) - 1 )² E₃(x) = Ln ( 1 - exp(-x) )
E₂(x) = x / ( exp(x) - 1 ) E₄(x) = E₂(x) - E₃(x)

Data Registers: /
Flags: /
Subroutines: "E^Z-1" E^Z Z*Z Z^2 Z/Z LNZ Z-Z

01 LBL "EMZ"
02 GTO IND Z
03 LBL 01
04 RCL Y
05 RCL Y
06 E^Z
07 Z*Z
08 Z*Z

09 R^
10 R^
11 XROM "E^Z-1"
12 Z^2
13 Z/Z
14 RTN
15 LBL 02
16 RCL Y

17 RCL Y
18 XROM "E^Z-1"
19 Z/Z
20 RTN
21 LBL 03
22 X<>Y
23 CHS
24 X<>Y

25 CHS
26 XROM "E^Z-1"
27 X<>Y
28 CHS
29 X<>Y
30 CHS
31 LNZ
32 RTN

33 LBL 04
34 XEQ 02
35 R^
36 R^
37 XEQ 03
38 Z-Z
39 END

( 67 bytes / SIZE 000 )

STACK	INPUTS	OUTPUTS
Z	m	/
Y	y	Im ( e^z - 1 )
X	x	Re ( e^z - 1 )

With z = x + i y and m = 1 , 2 , 3 , 4

Examples: z = 2 + 3 i

1 ENTER^ 3 ENTER^ 2 XEQ "EMZ" >>>>   0.658995987   X<>Y -1.198572640
2 ENTER^ 3 ENTER^ 2 XEQ "EMZ" >>>> -0.192258331   X<>Y -0.384898831
3 ENTER^ 3 ENTER^ 2 XEQ "EMZ" >>>>   0.125876182   X<>Y    0.016840416
4 ENTER^ 3 ENTER^ 2 XEQ "EMZ" >>>> -0.318134512   X<>Y -0.401739247

-Thus,

E₁(2+3.i) = 0.658995987 - 1.198572640 i
E₂(2+3.i) = -0.192258331 - 0.384898831 i
E₃(2+3.i) = 0.125876182 + 0.016840416 i
E₄(2+3.i) = -0.318134512 - 0.401739247 i

32°) Complex Planck Radiation Function

-Planck radiation function is defined as Plk(x) = 1 / [ x⁵ ( exp(1/x) - 1 ) ]

Data Registers: /
Flags: /
Subroutines: "E^Z-1" Z^2 Z*Z 1/Z

01 LBL "ZPLCK"
02 Z=0?
03 RTN
04 RCL Y
05 RCL Y
06 Z^2
07 Z^2
08 Z*Z
09 R^
10 R^
11 1/Z
12 XROM "E^Z-1"
13 Z*Z
14 1/Z
15 END

( 35 bytes / SIZE 000 )

STACK	INPUTS	OUTPUTS
Y	y	Im Plk (z)
X	x	Re Plk(z)

With z = x + i y

Example:

   0.2 ENTER^
   0.1 XEQ "ZPLCK"   >>>>   -13.00641622
                                    X<>Y   -221.0591332

-So Plk ( 0.1 + 0.2 i ) = -13.00641622 - 221.0591332 i

References:

[1] https://en.wikipedia.org/w/index.php?title=Differential_geometry_of_curves&oldid=740874421
[2] https://en.wikipedia.org/w/index.php?title=Frenet%E2%80%93Serret_formulas&oldid=731811089
[3] Elie Cartan - "Leçons sur la géométrie des espaces de Riemann" - Gauthier-Villars, Paris ( in French )
[4] Denis-Papin & Kaufmann - "Cours de calcul Tensoriel Appliqué" - Albin Michel ( in French )
[5] List of Formulas in Riemannian Geometry
[6] Nikodem J. Poplawski - "Spacetime and fields"- Department of Physics, Indiana University, Bloomington, IN 47405, USA
[7] http://www.wolframalpha.com/entities/mathworld/planck%E2%80%99s_radiation_function/mc/ph/55/
[8] Gerald Kaiser - "Quantum Physics, Relativity, and Complex Spacetime: Towards a New Synthesis"
[9] Giampiero Esposito - "From Spinor Geometry to Complex General Relativity"
[10] Jean E. Charon - "Complex General Relativity"

-The last reference is controversial, but there are many interesting ideas in Charon's theory...

hp41programs

Complex Derivatives & Riemannian Manifolds for the HP-41