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,

            ak bk = a1 b1 + a2 b2 + a3 b3

         gij ai aj = g11 a1 a1 + g12 a1 a2 + g13 a1 a3
                    + g21 a2 a1 + g22 a2 a2 + g23 a2 a3
                    + g31 a3 a1 + g32 a3 a2 + g33 a3 a3

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

        gij ai aj = g11 a1 a1 + g22 a2 a2 + g33 a3 a3  + 2 ( g12 a1 a2 + g13 a1 a3 + g23 a2 a3 )
 

-The partial derivative     m gij   means   gij / xm
-Likewise,          km gij   =  gij / xkxm
 
 
 
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  = ( ez - e-z )/2                           ArcSinh z  = Ln [ z + ( z2 + 1 )1/2 ]
    Cosh z = ( ez + e-z )/2                          ArcCosh z = Ln [ z + ( z + 1 )1/2 ( z - 1 )1/2 ]
    Tanh z = ( e2z - 1 )/( e2z + 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

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

Example:       z = 2+3.i

-Calculate  z2 ;  z1/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)2   =   -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  ;  zr = 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' ; zz'
 

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    zz' = 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.( f1 - f-1 ) - 600.( f2 - f-2 ) + 150.( f3 - f-3 ) - 25.( f4 - f-4 ) + 2.( f5 - f-5 ) ]   with  fk =  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 = x1 ,    R03-R04 = x2 , .......... ,    R2n-1-R2n = xn      n < 5

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

Flags:  /
Subroutine:   A program which computes  f(x1, ..... ,xn)  assuming  x1 , ......... , xn are in R01 thru R2n
 
 
      STACK        INPUTS      OUTPUTS
           Y            h     Im( f / xi )
           X             i     Re( f / xi )

 
Example:     f(x) = exp(-x2).Ln(y2+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 / x1 =  0.469272437 -  0.006070241 i

-Likewise,

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

  0.1   ENTER^
    3       R/S      >>>>   f 'x3 = f / x3 = 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 / x1 = 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

         d2f/dx2 = (1/25200.h2).[ -73766 f0 + 42000.( f1 + f-1 ) - 6000.( f2 + f-2 ) + 1000.( f3 + f-3 ) - 125.( f4 + f-4 ) + 8.( f5 + f-5 ) ]   with  fk =  f(x+k.h)

 and for the mixed derivative:

          f "xy = 2f / xy = (1/50400.h2).[ 73766 f00 + 42000.( f11 + f -1-1 - f10 - f -10 - f01 - f0-1 ) + 6000.( - f 22 - f -2-2 + f20 + f -20 + f02 + f0-2 )
                                                         + 1000.( f33 + f -3-3 - f30 - f -30 - f03 - f0-3 ) + 125.( - f 44 - f -4-4 + f40 + f -40 + f04 + f0-4 )
                                                               + 8.( f55 + f -5-5 - f50 - f -50 - f05 - f0-5 ) ]

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

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

          d2f/dx2 = (1/12h2).[ - f(x-2h) + 16.f(x-h) - 30.f(x) +16.f(x+h) - f(x+2h) ]

 and, for the crossed derivative

           f "xy = 2f / xy = (1/24.h2).[ 30 f00 + 16.( f11 + f -1-1 - f10 - f -10 - f01 - f0-1 ) + ( - f 22 - f -2-2 + f20 + f -20 + f02 + f0-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 = x1 ,    R03-R04 = x2 , .......... ,    R2n-1-R2n = xn      n < 5

                                           R10 = h
                                           R11-R12 =  2f / xixj        R13  ........  R21: temp
Flag:  F10
Subroutine:   A program which computes  f(x1, ..... ,xn)  assuming  x1 , ............ , xn are in R01 thru R2n
 
 
      STACK        INPUTS      OUTPUTS
           Z             h             /
           Y             j  Im ( 2f / xixj  )
           X             i  Re ( 2f / xixj  )

 
Example:    f(x) = exp(-x2).Ln(y2+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 = 2f / x2 =  -1.702735593 - 1.010546610 i

 Likewise,

 0.1   ENTER^
   1    ENTER^
   2       R/S     >>>>   f "xy = 2f / xy = 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 = 2f / x2 = -1.702735267 - 1.010547033 i                   ---Execution time = 14s---

   -0.01   ENTER^
       1      ENTER^
       2         R/S     >>>>   f "xy = 2f / xy = 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 )

       h3 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

       h3 f'''xxy ~  SUMk=1,2,...,5  Ak [ - 2 f0k + 2 f0-k  + ( fkk - f-k-k + f-kk - fk-k ) ]

 With    A1 = 5/6 ,  A2 = -5/84  ,  A3 = +5/756  ,  A4 = -5/8064  ,  A5 = +1/31500

       h3 f'''xyz ~  SUMk=1,2,...,5  Ak [ fk00 - f-k00 + f0k0 - f0-k0 + f00k - f00-k + fkkk - f-k-k-k - ( fkk0 - f-k-k0 + fk0k - f-k0-k + f0kk - f0-k-k ) ]

 With    A1 = 5/6 ,  A2 = -5/84  ,  A3 = +5/756  ,  A4 = -5/8064  ,  A5 = +1/31500
 
 

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

                                R01-R02 = x1 ,    R03-R04 = x2 , .......... ,    R2n-1-R2n = xn      n < 5

                                           R10 = h
                                           R11-R12 =  2f / xixjxk       R13  ........  R21: temp
Flags:  F09-F10
Subroutine:   A program which computes  f(x1, ..... ,xn)  assuming  x1 , ............ , xn are in R01 thru R2n
 
 
 
      STACK        INPUTS      OUTPUTS
           T             h             /
           Z             i             /
           Y             j  Im ( 2f / xixjxk )
           X             k  Re ( 2f / xixjxk )

 
Example:             f(x,y,z,t) = exp(-x2) Ln( y2 + z + t3 )     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 = 3f / xyz  = ?

-With   h = 0.1

     0.1  ENTER^
      1    ENTER^
      2    ENTER^
      3    XEQ "ZTD"   >>>>    f'''xyz = 3f / xyz  = 0.002045420 + 0.002544508 i                           ---Execution time = 3m49s---

    f'''xtt = 3f / xt2  = ?

-With   h = 0.1

     0.1  ENTER^
      1    ENTER^
      4    ENTER^    R/S   >>>>    f'''xtt = 3f / xt2  = - 0.028361907 - 0.027466199 i                                ---Execution time = 1m43s---

    f'''yyy = 3f / y3  = ?

-With   h = 0.1

     0.1  ENTER^
      2    ENTER^   ENTER^    R/S   >>>>    f'''yyy = 3f / y3  = - 0.002342127 - 0.001495556 i                       ---Execution time = 35s---
 

Hessian Matrix
 

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

       H = [ hi,j ]   with  hi,j = 2f / xixj
 

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

                                R01-R02 = x1 ,    R03-R04 = x2 , .......... ,    R2n-1-R2n = xn

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

Flag:  F10
Subroutine:   a program which computes  f(x1, ..... ,xn)  assuming  x1 , ............ , xn are in R01 thru R2n
 
 
      STACK        INPUTS      OUTPUTS
           Y             h             /
           X          n < 5        37.eee

 
Example:    f(x,y,z,t) = x2 y3 z4 t5     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 = x1 ,    R03-R04 = x2 , .......... ,    R2n-1-R2n = xn

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

Flag:  F10
Subroutine:   a program which computes  f(x1, ..... ,xn)  assuming  x1 , ............ , xn 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(-x2) Ln( y2 + z + t3 )     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) =  2f / x2 + 2f / y2 + 2f / z2 + 2f / t2 = -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

       D2 f = 4f / x4 + 4f / y4 + 4f / z4 + 2 ( 4f / x2y2 + 4f / x2z2 + 4f / y2z2 )
 

-"ZBHRM" uses the following approximation
 

   h4 D2 f  ~  A f000
               + B1 ( f100 + f-100 + f010 + f0-10 + f001 + f00-1 )  + C1 ( f110 + f-1-10 + f1-10 + f-110 + f101 + f-10-1 + f10-1 + f-101 + f011 + f0-1-1 + f01-1 + f0-11 )
               + B2 ( f200 + f-200 + f020 + f0-20 + f002 + f00-2 )  + C2 ( f220 + f-2-20 + f2-20 + f-220 + f202 + f-20-2 + f20-2 + f-202 + f022 + f0-2-2 + f02-2 + f0-22 )
               + B3 ( f300 + f-300 + f030 + f0-30 + f003 + f00-3 )  + C3 ( f330 + f-3-30 + f3-30 + f-330 + f303 + f-30-3 + f30-3 + f-303 + f033 + f0-3-3 + f03-3 + f0-33 )
               + B4 ( f400 + f-400 + f040 + f0-40 + f004 + f00-4 )  + C4 ( f440 + f-4-40 + f4-40 + f-440+ f404 + f-40-4 + f40-4 + f-404 + f044 + f0-4-4 + f04-4 + f0-44 )
               + B5 ( f500 + f-500 + f050+ f0-50 + f005 + f00-5 )  + C5 ( f550 + f-5-50 + f5-50 + f-550 + f505 + f-50-5 + f50-5 + f-505 + f055 + f0-5-5 + f05-5 + f0-55 )

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

     A = 41523361 / 540000  ,   B1 =  -4069 / 180             C1 =  +10 / 3
                                                 B2 =  +4969 / 1260          C2 =  -10 / 84
                                                 B3 =  -2201 / 3240           C3 =  +10 / 1134
                                                 B4 =  +371 / 4320            C4 =  -10 / 16128
                                                 B5 =  -5221 / 945000       C5 =  +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 = D2 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 ( D2 f )
           X             h      Re ( D2 f )

 
Example:      f(x) = exp(-x2).Ln(y2+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"  >>>>  D2 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  = 2f / x2 + 2f / y2 + 2f / z2      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  = 2f / r2 + (1/r) f/r + (1/r2) 2f / ¶f2 + 2f / z2      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  = 2f / r2 + (2/r) f/r + (1/r2) 2f / ¶q2 + (1/(r2tan q)f / ¶q + (1/(r2sin2q)) 2f / ¶f2      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(-x2) Ln(y2+z) , x2y2z2 , exp(x) y2z ]     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 z2 sin2f , r2 z , r3 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 sin2 q cos2 f , r2sin f  , r3 cos q cos2 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 KG and the mean curvature Km of a surface defined by a function  z = z(x,y)

-The Gaussian Curvature KG is be computed by:

   KG = [ ( 2z / x2 ) ( 2z / y2 ) - ( 2z / xy )2 ] / [ 1 + ( z / x )2 + ( z / y )2 ]2

-And the mean curvature Km is given by:

    Km = (1/2) [ t ( 1+p2 ) + r ( 1+q2 ) - 2 p q s ] / ( 1+p2+q2 )3/2

   where  p = z / x , q = z / y,  r = 2z / x2,  s = 2z / xy,  t = 2z / y2
 
 

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 = KG     R13-R14 = Km
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 Km
           Z             /       Re Km
           Y             /       Im KG
           X             h       Re KG

 
Example:   A surface is defined by  z(x,y) = Ln ( 1 + x2 y3 )    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,   KG = 0.034972708 - 0.037342284 i
 and   Km = -0.116193627 + 0.113513738 i

Note:

-With a sphere, you'll get   K =  1 / R2  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) = et cos t  ,  y(t) = et 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 z1(t) , z2(t) , z3(t) , z4(t)  assuming t is in registers R01-R02
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             /             /
           Y             /          tau(t)
           X             h          khi(t)

 
Example:   The curve defined by

     z1(t) = t4       z3(t) = t3      at the point  t = 0.7 + 0.8 i
     z2(t) = t2       z4(t) = et
 
 
 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 [ gij ]  i-e  [ gij ]
-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 gij 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  x1 , x2 , ........... , xn  in registers R01-R02 ,  R03-R04.............  R2n-1R2n   ( n < 5 ) and return  gij  in X & Y registers

-"ZINIGC"  calculates and stores the n(n+1)/2  gij  in registers  R61-R62 ................  ( i <= j )   for a given point  (  x1, ........... , xn )
-Then it calls "ZLS" to find the inverse of the metric matrix and stores gij  in registers  R61+n(n+1)  .....

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

-The Christoffel symbols  Gkij = (1/2) gkmi gmj + j gimm gij )   are then computed and stored in the following registers.

-Only  Gkij with  i <= j  because they are symmetric  Gkji = Gkij
 

>>> 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 = x1            R09 = n                              R61 to R60+n(n+1)                         =  gij    with  i <= j
                                        R03-R04 = x2            R10 = h                               R61+n(n+1) to R60+2n(n+1)          =  gij    with i <= j
                                         .............                                                                 R61+2n(n+1) to R60+n(n+1)(n+2) =  Gkij  with  i <= j
                                        R2n-1R2n = xn          R59-60 = g = det gij

Flags /
Subroutines:   The n(n+1)/2 functions    gij     with  i <= j
 
 
      STACK        INPUTS      OUTPUTS
           Y             h             /
           X          n < 5  bbb.eee(gij & Gkij)

   where n is the number of variables

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

    g11 = 1 + x2 y2           g14 =  x t     and the other    gij = 0
    g22 = 1 + x2 z2
    g33 = 1 + y2 z2
    g44 = 1 + x2 t2

-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 [ gij ]
                                                      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 = g11 = -6 - 24 i     R63-64 = g12 = 0                R67-68 = g13 = 0                 R73-74 = g14 = -3 + 4 i
                                               R65-66 = g22 = -6 - 24 i     R69-70 = g23 = 0                 R75-76 = g24 = 0
                                                                                          R71-72 = g33 =  -6 - 24 i     R77-78 = g34 = 0
                                                                                                                                      R79-80 = g44 =  -6 - 24 i
-In registers R81 thru R100

    R81-82 = g11 = -0.011247 +0.038444 i   rounded to 6D      ...............   R93-94 = g14 = -0.007464 + 0.003136 i

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

-And in registers R101 to R180, for example

   R101-R102 = G111 = 0.186872 - -0.412188 i

   R161-R162 = G411 = 0.000238574 - 0.003612692 i      .... and so on ....     until   R179-R180 = G444 = 0.098496984 - 0.392624655 i
 

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

Notes:

-Use  ZCIJK below to recall gij , gij  and  Gkij

-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 gij functions, all different from zero, V41 needs about 7 seconds to XEQ "ZINIGC"
 

Recalling Gij & Cijk
 

-ZCIJK is an M-Code routine that recalls  Gkij provided  R09 = n
-This register has been initialized by "ZINIGC"
 
 
      STACK        INPUTS      OUTPUTS
           Z             i     address-60
           Y             j       Im  Gkij
           X             k       Re  Gkij

   Where  "address" is the address of the data register that contains  Im  Gkij

Example1:

-With the metric above and after executing "ZINIGC"

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

 So,    G233 =  -0.186274510 + 0.411764706 i

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

Example2:     ZCIJK  may also be used to recall  gij  if  k = -1  and  gij  if  k = 0

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

 and   2  ENTER^   ENTER^   1  CHS   ZCIJK    yields   g22 = - 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 Gkij  is

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

Curvature Tensor ( 0-4 )
 

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

     Rijkl = (1/2) ( jk gil + il gjk - jl gik - ik gjl ) + gmp ( GmjkGpil - GmjlGpik )
 
 
 
      STACK        INPUTS      OUTPUTS
           T             i            /
           Z             j            /
           Y             k       Im  Rijkl
           X             l       Re  Rijkl

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

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

-Likewise

  1  ENTER^
  2  ENTER^
  2  ENTER^
  4     R/S      >>>>  R1114 = 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 =  n2 ( n2 - 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:

      Rlijk = glm Rmijk

-"ZRIJK^L" calculates these components
 
 
      STACK        INPUTS      OUTPUTS
           T             i            /
           Z             j            /
           Y             k       Im  Rlijk
           X             l       Re  Rlijk

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

  2  ENTER^
  1  ENTER^
  2  ENTER^
  1  XEQ "ZRIJK^L"  >>>>  R1212 = -0.127624 - 0.158155 i                                     ---Execution time = 84s---

-Likewise

  3  ENTER^
  1  ENTER^
  2  ENTER^
  4     R/S      >>>>  R4312 = 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

    Rij = gkm  Rikjm

-We can also contract other indices, but it would yield  -Rij
 
 
      STACK        INPUTS      OUTPUTS
           Y             i       Im  Rij
           X             j       Re  Rij

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

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

-Likewise

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

Scalar Riemannian Curvature
 

-Contracting the Ricci tensor, we get the scalar curvature    R = gij  Rij
-"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 gij functions, all different from zero, V41 needs about 30 seconds to XEQ "ZR"
 

Einstein Tensor
 

-Einstein tensor is a symmetric tensor defined as

   Eij = Rij - (1/2) R gij

-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  Eij
           X             j         Re  Eij

 
Examples:    The same metric at the same point

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

-Likewise you will find

  2  ENTER^
  3     R/S      >>>>  E23 = -0.127501 - 0.157186 i

-And  4  ENTER^  R/S  >>>>    E44 = 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     Eij= Rij - (1/2) R gij + L gij

-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

    Di Tij = 0  for all j    where    Di is the covariant differentiation.

-They are the equations of General Relativity !
 

Weyl Tensor
 

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

   Wijkl = Rijkl + [ 1/(n-2) ] ( Ril gjk - Rik gjl + Rjk gil - Rjl gik ) + [ R/(n-1)/(n-2) ] ( gik gjl - gil gjk )

-Like Einstein tensor, execute "ZR" first before executing "ZWIJKL" so that data registers R57-R58 = R
 
 
      STACK        INPUTS      OUTPUTS
           T             i       Im Rijkl
           Z             j       Re Rijkl
           Y             k       Im Wijkl
           X             l       Re Wijkl

 
Examples:    The same metric at the same point

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

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

-Likewise

  2  ENTER^
  3  ENTER^
  1  ENTER^
  3     R/S      >>>>  W2313 = 0.061810 + 0.096108 i

 And in Z&T-registers:     R2313 = -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  Rijkl

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

-If n = 4,  h = +0.1  and with relatively simple gij 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:

    Dk Vi =    k Vi - Gmik Vm

Contravariant Vector:

    Dk Vi =    k Vi + Gikm Vm
 

-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)  =  gij ,  gijGkij

                                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

    V1 = x2 + y z t      V2 = x2 y2 t3     V3 =  x y z t    V4 =  x2 y2 z2 t2
 

-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 [ Di Vj ] , 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,  D3 V1 = -3 + 4 i

    Contravariant Vector Field

-If the same functions are the components of a contravariant vector field:    V1 = x2 + y z t      V2 = x2 y2 t3     V3 =  x y z t    V4 =  x2 y2 z2 t2
 

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

-And we obtain the coefficients of the matrix [ Di Vj ] , 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,  D3 V1 = -3 + 4 i
 

Curl of a Covariant Vector Field
 

-The partial derivatives of a vector is not a tensor   and  Curlij V is defined as   Di Vj - Dj Vi
-But if we apply the formulae given in the paragraph above, the Christoffel symbols vanish and finally it yields

     Curlij V = i Vj - j Vi

-This tensor is obviously antisymmetric, so  Curlii V = 0  and  Curlji V = - Curlij V

-"ZCURL"  calculates and stores all the coefficients of the matrix  [ Curlij 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)  =  gij ,  gijGkij

                                    R09 = n      R20 to R26 are used by "ZCURL"
                                    R10 = h
Flags: /
Subroutine:         "ZFD"
 
 
 
      STACK        INPUTS      OUTPUTS
           X        bbb.eee     bbb.eeeCURL
           L             /        bbb.eee

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

    V1 = x2 + y z t      V2 = x2 y2 t3     V3 =  x y z t    V4 =  x2 y2 z2 t2

-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,       Curl13 V = 1 V3 - 3 V1 = - 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 Vk by the modulus of the corresponding vector ek of the basis
  or the quotient of the covariant component Vk by the same quantity:

    Vk ek = Vk / ek = (gkk)1/2 Vk = (gkk) -1/2 Vk         ( No summation here )
 

Divergence of a Contravariant Vector Field
 

-The divergence of a vector field V is given by

    Div V =     Di Vi =     i Vi + Giik Vk

-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)  =  gij ,  gijGkij

                                    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:

  V1 = x2 + y z t      V2 = x2 y2 t3     V3 =  x y z t    V4 =  x2 y2 z2 t2

-> 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  Gradk 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) = gjk ( jk f  -  Gijk 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)  =  gij ,  gijGkij

                                    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 x1 is in R01-R02 ..................  xn is in R2n-1R2n
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             /     bbb.eeeGRAD
           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 )  = x2 + y z t    i-e    V1  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)  =  gij ,  gijGkij
                                    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   ...Uk-1.Uk-1-1....U2-1.U1-1.A.U1.U2....Uk-1.Uk....
 where the Uk are unitary matrices. The eigenvectors are the columns of  U = U1.U2....Uk-1.Uk.... 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 ai j below the main diagonal
-Then,  U1  is determined so that it places a zero in position ( i , j )  in U1-1.A.U1

  U1  has the same elements as the Identity matrix, except that  ui i = uj j = x  and  ui j = y + i.z ,  uj i = -y + i.z

     with  y + i.z  = C.x  ,  x = ( 1 + | C |2 ) -1/2  and   C = 2.ai j / [ ai i - ai j + ( ( ai i - ai j )2 + 4 ai j aj i )1/2 ]

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

-The successive greatest | ai 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 = [ wi j ]  by

   wi j = 0    if  i > j
   wi i = 1
   wi j = Sumk = i+1, i+2 , .... , j      ti k wk j / ( tj j - ti i )      if  i < j

-Then, the n columns of  U.W  are the eigenvectors corresponding to the n eigenvalues  t11 = k1 , ............... , tnn = kn
 

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

                               R01,R02 = a11     .................        R2n2-2n+1,R2n2-2n+2 = a1n
                              .................................................................................................                        ( the n2 elements of A in column order )
   INPUTS
                               R2n-1,R2n = an1  .................        R2n2-1,R2n2 = ann

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

   DURING            R01 thru R2n2 = the "infinite" product   ...Uk-1.Uk-1-1....U2-1.U1-1.A.U1.U2....Uk-1.Uk....
       THE                R2n2+1 thru R4n2 =  U1.U2........Uk.....
 ITERATIONS      R4n2+1 thru R4n2+2n: temp  ( for non-Hermitian matrices only )

  ---------

                             R00 = n
                             R01,R02 = k1          R2n+1,R2n+2 = V1(1)      ....................    R2n2+1,R2n2+2 = V1(n)
                             .......................            ........................                ....................     ..................................
  OUTPUTS
                             R2n-1,R2n = kn        R4n-1,R4n =    Vn(1)      .....................   R2n2+2n-1,R2n2+2n = Vn(n)
 

                            ( n eigenvalues     ;        1st eigenvector      ;  .................................. ;    nth eigenvector )
 

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

  ( SIZE 4n2+1 if A is Hermitian / SIZE 4n2+2n+1 if A is non-Hermitian )
 
 
      STACK        INPUTS      OUTPUTS
           Y             /             y1
           X             /            x1

  where  k1 = x1 + y1  = 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---

           k1 =  15.61385274 + 0.i  = ( X , Y ) = ( R01 , R02 )
           k2 =  5.230678473 + 0.i  = ( R03 , R04 )
           k3 = -6.844531166 + 0.i  = ( R05 , R06 )

  and the 3 corresponding eigenvectors:

    V1( 0.481585120 + 0.108201123 i , 0.393148652 + 0.527305194 i , -0.142958847 + 0.550739893 i )    =  ( R07 R08 , R09 R10 , R11 R12 )
    V2( -0.436763638 + 0.075843695 i , 0.210271354 - 0.463555613 i , -0.516799072 + 0.526598642 i )   =  ( R13 R14 , R15 R16 , R17 R18 )
    V3( -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---

           k1 = 7.656606017 + 15.61073839 i  = ( X , Y ) = ( R01 , R02 )
           k2 = 1.661248139 - 1.507335318 i   = ( R03 , R04 )
           k3 = -3.317854157 - 2.103403077 i  = ( R05 , R06 )   and the 3 corresponding eigenvectors

    V1( 0.523491429 + 0.008555748 i , 0.650242685 - 0.046353000 i , 0.546288478 + 0.049882590 i )      =  ( R07 R08 , R09 R10 , R11 R12 )
    V2( -0.4777850326 - 0.196368365 i , 0.660547554 - 0.265888146 i , -0.356842635 + 0.318267519 i )  =  ( R13 R14 , R15 R16 , R17 R18 )
    V3( -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 ( x2 + y2 )  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  = ( x2 + y2 + z2 )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( x1 ,..........., xn )  V'(  y1 ,........, yn )     is defined as usual by

          < V , V' > = x1 y1 + ................. + xn yn

-Whereas the Hermitian product is     < V , V' >H = x1 y1bar + ................. + xn ynbar      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...