ZDERIVE41 Module

Overview

-This 8K-module extends most of the programs in DERIVE41 to complex holomorphic functions and complex Riemannian manifolds:

   Several differential operators
   Usual problems in Euclidean differential geometry
   Tensors in Riemaniann geometry ( Complex General Relativity )
   Complex linear systems and eigenvalues & eigenvectors of complex matrices ( hermitian or not )

There are also programs to compute elemetary complex functions ( exponential, logarithm, circular & hyperbolic functions )

Last Update:

-I've added the vector Laplacian in the "ZCDGL" program
-The formulas are relatively complicated in cylindrical and even more in spherical coordinates,
so I had to remove "ZPLCK" "ZIN" and "ZOUT" to free enough space.

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.

About Flags:

F09-F10

-F09 is used by "ZSD" and F09-F10 by "ZTD" to compute 2nd & 3rd derivatives

F07

-This flag is used by several circular & hyperbolic functions

F01-F02

-With"ZCDGL" ,

Set F01 for cylindrical coordinates
Set F02 for spherical coordinates

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

XROM Function Desciption

21,00
21,01
21,02
21,03
21,04
21,05
21,06
21,07
21,08
21,09
21,10
21,11
21,12
21,13
21,14
21,15
21,16
21,17
21,18
21,19
21,20
21,21
21,22
21,23
21,24
21,25
21,26
21,27
21,28
21,29
21,30
21,31
21,32
21,33
21,34
21,35
21,36
21,37
21,38
21,39
21,40
21,41
21,42
21,43
21,44
21,45
21,46
21,47
21,48
21,49
21,50
21,51
21,52
21,53
21,54
21,55
21,56
21,57
21,58
21,59
21,60
21,61
21,62
21,63 -ZDERIVE41
Z+Z
Z-Z
Z*Z
Z/Z
1/Z
Z^2
SQRTZ
E^Z
LNZ
"Z^X"
"Z^Z"
"SHZ"
"CHZ"
"THZ"
"ASHZ"
"ACHZ"
"ATHZ"
"SINZ"
"COSZ"
"TANZ"
"ASINZ"
"ACOSZ"
"ATANZ"
Z=0?
"ZFD"
"ZSD"
"ZTD"
"ZHESS"
"ZLAPN"
"ZBHRM"
"ZCDGL"
"ZKHT"
"ZKGM"
"ZLS"
-ZRIEMANN
"ZINIGC"
ZCIJK
"ZRIJKL"
"ZRIJK^L"
"ZRIJ"
"ZR"
"ZEIJ"
"ZWIJKL"
"ZCURL"
"ZDIV"
"ZLAPG"
"ZDCOV"
"ZDCOV^"
"ZV-V^"
"ZV^-V"
-MISCELLAN
"ZEGVL"
"ZEGVH"
"ZKHT4"
"ZP-R"
"ZR-P"
"ZS-R"
"ZR-S"
"ZDOT"
"ZDOTH"
"1"
"2"
"3" Section Header
Sum of 2 complexes
Difference between 2 complexes
Product
Quotient
Inverse
Square
Square Root
Exponential
Logarithm
Raising a Complex to a Real exponent
Raising a Complex to a Complex exponent
Hyperbolic Sine
Hyperbolic Cosine
Hyperbolic Tangent
Inverse hyperbolic Sine
Inverse hyperbolic Cosine
Inverse hyperbolic Tangent
Sine
Cosine
Tangent
Inverse Sine
Inverse Cosine
Inverse Tangent
Skips the next line in a program if z # 0
First Derivatives of a Function of N variables
Second Derivatives of a Function of N variables
Third Derivatives of a Function of N variables
Hessian Matrix
Laplacian of a Function of N variables
Biharmonic Operator ( 3 Dim. )
Curl, Divergence, Gradient & Laplacian ( 3D )
Curvature & Torsion of a Curve ( 3 Dim )
Gaussian Curvature & Mean Curv. of a Surface.
Complex Linear Systems
Section Header
Initialization ( Metric Tensor & Christoffel Symbols )
Recalling Gij & Cijk
Curvature Tensor ( 0-4 )
Curvature Tensor ( 1-3 )
Ricci Tensor
Scalar Riemannian Curvature
Einstein Tensor
Weyl Tensor
Curl of a Covariant Vector Field
Divergence of a Contravariant Vector Field
Laplacian and Gradient of a Scalar Field
Covariant Derivative of a covariant vector field
Covariant Derivative of a contravariant vector field
Covariant Vector -> Contravariant Vector
Contravariant Vector -> Covariant Vector
Section Header
Eigenvalues & eigenvector of a complex matrix
Eigenvalues & eigenvector of a hermitian matrix
Curvature & Torsion of a 4D-curve
Complex Polar -> Rectangular conversion
Complex Rectangular -> Polar conversion
Complex Spherical -> Rectangular Conversion
Complex Rectangular -> Spherical conversion
Complex Dot-Product
Complex Hermitian Product
3 Subroutines
called by
"ZCDGL"

GIJ name in Alpha register

-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

Complex Operations & Elementary Functions

Z+Z Z-Z Z*Z Z/Z 1/Z Z^2 SQRTZ 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.

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

1°) First case: Function of one complex number

      STACK        INPUTS      OUTPUTS

           Y             y             v

           X             x             u

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

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

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

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

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

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 )

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

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

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

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

      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.

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

"LAPN" 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

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

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

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.

      STACK        INPUTS     OUTPUTS

           Y             /      Im ( D² f )

           X             h      Re ( D² f )

"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.

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

      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

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

      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.

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

      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

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

      STACK        INPUTS      OUTPUTS

           Z             /             /

           Y             /          tau(t)

           X             h          khi(t)

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 = 3m36s---
                                       RDN   0.370021140 = R12
                                       RDN   0.143785368 = R13
                                       RDN   0.048055232 = R14

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

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

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

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

      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"

Recalling Gij & Cijk

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

      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 )

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 )

      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

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

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

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

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

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

      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"

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

      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 !

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

      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"

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"

     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

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"

      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 )

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"

      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.

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

      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^"

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

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

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: /

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

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

      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

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"

      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 )

Complex Dot-Product&Complex 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

      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.

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

XROM	Function	Desciption
21,00 21,01 21,02 21,03 21,04 21,05 21,06 21,07 21,08 21,09 21,10 21,11 21,12 21,13 21,14 21,15 21,16 21,17 21,18 21,19 21,20 21,21 21,22 21,23 21,24 21,25 21,26 21,27 21,28 21,29 21,30 21,31 21,32 21,33 21,34 21,35 21,36 21,37 21,38 21,39 21,40 21,41 21,42 21,43 21,44 21,45 21,46 21,47 21,48 21,49 21,50 21,51 21,52 21,53 21,54 21,55 21,56 21,57 21,58 21,59 21,60 21,61 21,62 21,63	-ZDERIVE41 Z+Z Z-Z Z*Z Z/Z 1/Z Z^2 SQRTZ E^Z LNZ "Z^X" "Z^Z" "SHZ" "CHZ" "THZ" "ASHZ" "ACHZ" "ATHZ" "SINZ" "COSZ" "TANZ" "ASINZ" "ACOSZ" "ATANZ" Z=0? "ZFD" "ZSD" "ZTD" "ZHESS" "ZLAPN" "ZBHRM" "ZCDGL" "ZKHT" "ZKGM" "ZLS" -ZRIEMANN "ZINIGC" ZCIJK "ZRIJKL" "ZRIJK^L" "ZRIJ" "ZR" "ZEIJ" "ZWIJKL" "ZCURL" "ZDIV" "ZLAPG" "ZDCOV" "ZDCOV^" "ZV-V^" "ZV^-V" -MISCELLAN "ZEGVL" "ZEGVH" "ZKHT4" "ZP-R" "ZR-P" "ZS-R" "ZR-S" "ZDOT" "ZDOTH" "1" "2" "3"	Section Header Sum of 2 complexes Difference between 2 complexes Product Quotient Inverse Square Square Root Exponential Logarithm Raising a Complex to a Real exponent Raising a Complex to a Complex exponent Hyperbolic Sine Hyperbolic Cosine Hyperbolic Tangent Inverse hyperbolic Sine Inverse hyperbolic Cosine Inverse hyperbolic Tangent Sine Cosine Tangent Inverse Sine Inverse Cosine Inverse Tangent Skips the next line in a program if z # 0 First Derivatives of a Function of N variables Second Derivatives of a Function of N variables Third Derivatives of a Function of N variables Hessian Matrix Laplacian of a Function of N variables Biharmonic Operator ( 3 Dim. ) Curl, Divergence, Gradient & Laplacian ( 3D ) Curvature & Torsion of a Curve ( 3 Dim ) Gaussian Curvature & Mean Curv. of a Surface. Complex Linear Systems Section Header Initialization ( Metric Tensor & Christoffel Symbols ) Recalling Gij & Cijk Curvature Tensor ( 0-4 ) Curvature Tensor ( 1-3 ) Ricci Tensor Scalar Riemannian Curvature Einstein Tensor Weyl Tensor Curl of a Covariant Vector Field Divergence of a Contravariant Vector Field Laplacian and Gradient of a Scalar Field Covariant Derivative of a covariant vector field Covariant Derivative of a contravariant vector field Covariant Vector -> Contravariant Vector Contravariant Vector -> Covariant Vector Section Header Eigenvalues & eigenvector of a complex matrix Eigenvalues & eigenvector of a hermitian matrix Curvature & Torsion of a 4D-curve Complex Polar -> Rectangular conversion Complex Rectangular -> Polar conversion Complex Spherical -> Rectangular Conversion Complex Rectangular -> Spherical conversion Complex Dot-Product Complex Hermitian Product 3 Subroutines called by "ZCDGL"

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