DIFF GEOM

# Differential Geometry(2) for the HP41

Overview

0°)  2 M-Code Routines
1°)  Curvature & Torsion of a Parametric Curve
2°)  Gaussian Curvature of a Surface
3°)  Biharmonic Operator
4°)  Curl, Divergence, Gradient & Laplacian
5°)  Laplacian of a function of n variables ( n < 10 )
6°)  Linear Systems
7°)  Riemannian Geometry

a) Initialization
b) Recalling Gij & Cij^k
c) Curvature Tensor
d) Ricci Tensor
e) Scalar Curvature
f) Einstein Tensor
g) Weyl Tensor
h) Covariant Derivative of a Vector or a Tensor of order 2
i) Curl, Divergence, Laplacian & Gradient
j) Covariant Vector <> Contravariant Vector

8°)  An Example of a Non-Diagonal Metric
9°)  Schwarzchild Metric
10°) Non-Riemannian Geometry

a) Curvature Tensor(s)
b) Ricci Tensor & Segmental Curvature
c) Scalar Curvature(s)
d) Torsion Tensor
e) Torsion Vector

11°)  Differential Geometry QRG

-Almost all these programs require the evaluation of derivatives
-Formulas of order 10 are usually used, but the execution time may be prohibitive to compute all tensors..
-So, you can also choose 4th-order methods, to get faster - but less accurate - results.

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

-Several programs come from "Numerical Differentiation for the HP41", sometime with minor modifications...
-All are included in a module published here.

>>> Last update:

-I've slightly modified several programs.
-Instead of using F00 for diagonal metrics, I've used a trick to get the derivative faster when the function = 0
-Simply store an alpha string in X-register and "FD" ( if CF 03 ) or "SD" will return 0 without evaluating the function 10 times.

-Several routines have been added to calculate a few tensors in non-Riemannian manifolds.

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

0°) GIJ name & CIJK name in Alpha register

-The M-code routine Fµ places "GIJ" in the alpha register where  I = Inf ( i , j )  &  J = Sup ( i , j )
if X-register = i  and Y-register = j ,
-So, the program lines ( with FIX 0  CF 29 ):

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

are replaced by

 Fµ

-A similar function ( called "DIFF MANIF" - a section header in the module )
places CIJK in alpha provided Z = i  Y = j  X = k

086   "F"
009   "I"
00E  "N"
001  "A"
00D  "M"
020   " "
006   "F"
006   "F"
009   "I"
004   "D"
02D  "-"
288   SETF 7
023   JNC +04
0CC  "µ"
006   "F"
284   CLRF 7
04E  C=0
168  M=C
1A8  N=C
1E8  O=C
228  P=C
28C  ?FSET 7
033  JNC +06
29C  PT=7
110  LD@PT- 4
0D0  LD@PT- 3
0D0  LD@PT- 3
023  JNC+04
09C  PT=5
110  LD@PT- 4
1D0  LD@PT- 7
01C  PT=3
0D0  LD@PT- 3
010  LD@PT- 0
0D0  LD@PT- 3
28C  ?FSET 7
07B  JNC+15d
10E  A=C ALL
078  C=Z
13C  RCR 8
05C  PT=4
102  A=C @PT
0B8  C=Y
0FC  RCR 10
21C  PT=2
102  A=C @PT
0F8  C=X
37C  RCR 12
39C  PT=0
102  A=C @PT
083  JNC+16d
0EE B<>C ALL
0B8  C=Y
10E  A=C ALL
0F8  C=X
30E  ?A<C ALL
017  JC+02
0AE  A<>C ALL
06E  A<>B ALL
37C  RCR 12
39C  PT=0
102  A=C @PT
0EE  B<>C ALL
21C  PT=2
0FC  RCR 10
102  A=C @PT
0AE  A<>C ALL
168  M=C
3E0  RTN

( 69 words )

1°) Curvature & Torsion of a Parametric Curve

-A curve  may be defined by 3 functions X(t) Y(t) Z(t)  in a 3-dimensional Euclidean space.
-"KHT" returns the curvature c ( khi ) and the torsion t ( tau ) of the curve at a point defined by t

-The curvature is calculated by   c = | r' x r'' | / | r' |3                                      where  the vector  r = [ X(t) , Y(t) , Z(t) ]
and the torsion by   t = ( r' x r'' ) .r''' / | r' x r'' |2                        x = cross-product  and  . = dot-product

Data Registers:              R00 =  | r' |

R04 = z'        R07 = y'        R10 = | r' x r'' |2       R13 = x'        R16 = x
R05 = z''       R08 = y''        R11 = t                    R14 = x''
R06 = z'''      R09 = y'''       R12 = h                    R15 = x'''
Flags: /
Subroutines:    "dF"  and  3 functions named  "X"  "Y"  "Z"   that compute x(t) , y(t) , z(t)  assuming t is in X-register upon entry

 01 LBL "KHT"  02 LBL 13  03 "Z"  04 ASTO 00  05 XROM "dF"  06 STO 12   07 RDN  08 STO 13  09 X<>Y  10 STO 14  11 "Y"  12 ASTO 00  13 RCL 02  14 RCL 01  15 XROM "dF"  16 STO 09  17 RDN  18 STO 10         19 X<>Y  20 STO 11  21 "X"  22 ASTO 00  23 RCL 02   24 RCL 01  25 XROM "dF"  26 RCL 12  27 RCL 09  28 R-P  29 RCL 03  30 R-P  31 STO 00  32 RCL 13  33 RCL 09  34 *  35 RCL 12  36 RCL 10  37 *  38 -  39 X^2  40 STO 06  41 LASTX  42 RCL 05  43 STO 08  44 *  45 RCL 12  46 RCL 04  47 STO 07  48 *  49 RCL 13  50 RCL 03  51 *  52 -  53 ENTER  54 X^2  55 ST+ 06  56 CLX  57 RCL 11  58 * 59 +  60 RCL 10  61 RCL 03  62 *  63 RCL 09  64 RCL 07   65 *  66 -  67 ENTER  68 X^2  69 ST+ 06  70 CLX  71 RCL 14  72 *  73 +  74 RCL 06  75 ST/ Y  76 SQRT  77 RCL 00   78 3  79 Y^X  80 /  81 X<>Y  82 STO 04         83 X<>Y  84 STO 03  85 RTN  86 GTO 13  87 LBL "dF"  88 LBL 10  89 STO 01  90 X<>Y  91 STO 02  92 X<>Y  93 XEQ IND 00  94 STO 06  95 CLX  96 STO 08  97 XEQ 01  98 STO 04  99 RCL 07 100 STO 03 101 STO 05 102 42 103 ST* 04 104 ST+ X 105 ST* 03 106 70098 107 CHS 108 ST* 05 109 XEQ 01 110 6 111 * 112 ST- 04 113 RCL 07 114 24 115 * 116 ST- 03 117 2184.5 118 * 119 ST+ 05 120 XEQ 01 121 ST+ 04 122 RCL 07  123 6 124 * 125 ST+ 03 126 2434.5 127 * 128 ST- 05 129 XEQ 01 130 8 131 / 132 ST- 04 133 RCL 07  134 ST- 03 135 2522 136 * 137 ST+ 05 138 25 139 ST* 03 140  E3 141 ST* 04 142 XEQ 01 143 8 144 * 145 ST+ 04 146 RCL 07        147 ST+ X 148 ST+ 03 149 RCL 07 150 205 151 * 152 ST- 05 153 RCL 02 154 ST/ 04 155 2520 156 * 157 ST/ 03 158 10 159 * 160 ST/ 04 161 1.2 162 * 163 RCL 02 164 X^2 165 * 166 ST/ 05 167 RCL 05 168 RCL 04 169 RCL 03 170 RTN 171 GTO 10 172 LBL 01 173 RCL 02 174 ST+ 08 175 RCL 01 176 RCL 08 177 + 178 XEQ IND 00 179 STO 07 180 RCL 01 181 RCL 08 182 - 183 XEQ IND 00 184 RCL 07  185 X<>Y 186 ST- 07 187 + 188 RCL 06 189 ST+ X 190 - 191 RTN 192 LBL 02 193 RCL 14 194 RCL 15 195 X<> Z 196 GTO IND 16 197 LBL 03 198 RCL 15  199 X<>Y 200 RCL 13        201 GTO IND 16 202 LBL 04 203 RCL 14 204 RCL 13 205 GTO IND 16 206 LBL 07 207 RCL 11 208 X<>Y 209 GTO IND 12 210 LBL 08 211 RCL 10 212 GTO IND 12 213 LBL "dF2" 214 LBL 11 215 XROM "MDV" 216 STO 09 217 RCL 01 218 STO 10 219 RCL 02 220 STO 11 221 8 222 X<> 00 223 STO 12 224 RCL 03 225 RCL 11 226 XROM "dF" 227 STO 13 228 RDN 229 STO 14 230 X<>Y 231 STO 15 232 DSE 00 233 RCL 02 234 RCL 10 235 XROM "dF" 236 RCL 12 237 STO 00 238 RCL 15 239 STO 08 240 RCL 09 241 SIGN 242 RCL 14  243 STO 07 244 RCL 04 245 RCL 13 246 STO 06 247 RCL 03 248 RTN 249 GTO 11 250 LBL "dF3" 251 LBL 05 252 STO 13 253 RDN 254 STO 14        255 RDN 256 STO 15 257 X<>Y 258 STO 02  259 4 260 X<> 00 261 STO 16 262 RCL 02 263 RCL 15 264 XROM "dF" 265 STO 09 266 RDN 267 STO 10 268 X<>Y 269 STO 11 270 DSE 00 271 RCL 02 272 RCL 14 273 XROM "dF" 274 STO 17 275 RDN 276 STO 18 277 X<>Y 278 STO 19 279 DSE 00 280 RCL 02 281 RCL 13 282 XROM "dF" 283 RCL 16 284 STO 00 285 RCL 19 286 STO 08 287 RCL 06 288 STO 12 289 RCL 04 290 RCL 18 291 STO 07 292 + 293 RCL 10 294 + 295 RCL 09 296 RCL 17 297 STO 06 298 RCL 03  299 RTN 300 GTO 05 301 LBL "KGS" 302 LBL 06 303 XROM "dF2" 304 X^2 305 X<>Y 306 X^2 307 + 308 1 309 + 310 X^2 311 X<> Z 312 * 313 RCL 09        314 X^2 315 - 316 X<>Y 317 / 318 RTN 319 GTO 06 320 LBL "MDV" 321 LBL 12 322 STO 01 323 RDN 324 STO 02 325 X<>Y 326 STO 03 327 X<> T 328 XEQ IND 00 329 STO 04 330 CLX 331 STO 06 332 XEQ 01 333 42 334 * 335 STO 05 336 XEQ 01 337 6 338 * 339 ST- 05 340 XEQ 01 341 ST+ 05 342 XEQ 01 343 8 344 / 345 ST- 05 346  E3 347 ST* 05 348 XEQ 01 349 8 350 * 351 RCL 05 352 + 353 50400 354 RCL 03 355 X^2 356 * 357 / 358 STO 05  359 RTN 360 GTO 12 361 LBL 01 362 RCL 03 363 ST+ 06 364 RCL 02 365 RCL 01 366 RCL 06 367 + 368 XEQ IND 00 369 STO 07 370 RCL 02        371 RCL 01 372 RCL 06 373 - 374 XEQ IND 00 375 ST+ 07 376 RCL 02 377 RCL 06 378 + 379 RCL 01 380 XEQ IND 00 381 ST+ 07 382 RCL 02 383 RCL 06 384 - 385 RCL 01 386 XEQ IND 00 387 ST+ 07 388 RCL 02 389 RCL 01 390 RCL 06 391 ST+ Z 392 + 393 XEQ IND 00 394 ST- 07 395 RCL 02 396 RCL 01 397 RCL 06 398 ST- Z 399 - 400 XEQ IND 00 401 RCL 07 402 - 403 RCL 04 404 ST+ X 405 + 406 RTN 407 GTO 12 408 END

( 613 bytes / SIZE var. )

 STACK INPUTS OUTPUTS Y h t (t) X t c (t)

Example:   The curve defined by   x(t) = et cos t  ,  y(t) = et sin t  ,  z(t) = Ln(t)        ( t expressed in radians )

-Evaluate the curvature c (t) and the torsion t (t) for t = 1.3

-Load the following routines in main memory

 01  LBL "X"  02  E^X  03  LASTX  04  COS  05  *  06  RTN  07  LBL "Y"  08  E^X  09  LASTX  10  SIN  11  *  12  RTN  13  LBL "Z"  14  LN  15  RTN

XEQ "RAD"   and if you choose  h = 0.1

0.1   ENTER^
1.3    XEQ "KHT"  >>>> c (1.3)  = 0.194807823                       ---Execution time = 45s---
X<>Y    t (1.3)  = 0.123665258

-The exact values are   c (1.3) = 0.194807823  &  t (1.3) = 0.123665529

Notes:

-If you want to estimate the curvature and torsion with other parameters, simply press   h ENTER^  t   R/S

-With an helicoidal curve:    x(t) = r Cos t , y(t) = r Sin t ,  z(t) = a t
where r & a are positive constants,

c = r / ( a2 + r2 )
= a / ( a2 + r2 )

-This long program also contains "KGS" to calculate the Gaussian Curvature of a surface ( 3D)
and several subroutines for numerical differentiation.

2°) Gaussian Curvature of a Surface

-We assume that the surface is defined by a function of 2 variables  z = z(x,y)

-The Gaussian Curvature K may be computed by:

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

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

R01 to R15: temp
Flags: /
Subroutines:  "dF" "dF2" "MDV" and a program that takes x in X-register and y in Y-register and returns f(x,y) in X-register

 STACK INPUTS OUTPUTS Z h / Y y / X x Gauss Curv.

Example:   Calculate the Gaussian Curvature K of the surface defined by z = Exp(-x.y) Ln( 2 + x2 + y2 )  at the point  x = 1 , y = 2

-Load the short routine below in main memory

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

Z  ASTO 00

If you choose  h = 0.1

0.1  ENTER^
2    ENTER^
1    XEQ "KGS"   >>>>    K = 0.057571609                                              ---Execution time = 78s---

-With  x = y = 1 , you'll get   K = -0.106545128

-The exact results are

K = + 0.0575715829 ...
K = - 0.1065451428 ....

Notes:

-If you want to estimate the Gaussian curvature, simply press   h  ENTER^  y  ENTER^  x   R/S

-If  K > 0  we have an elliptic point
-If  K < 0  we have a hyperbolic point
-If  K = 0  we have a parabolic point

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

Derivatives of a Function of 1 variable

-"dF"  calculates the 1st,  2nd & 3rd derivatives of a function of 1 variable f(x)

Formulae:

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 ) ] + O(h10)
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 ) ] + O(h10)

(  f(x+k.h) is denoted fk to simplify these expressions )

d3f/dx3 = [ -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) ] / h3 + O(h10)

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

-They are exact for any polynomial of degree < 11

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

R01 = x     R03 = df/dx    R05 = d3f/dx3    R06 to R08: temp
R02 = h     R04 = d2f/dx2
Flags: /
Subroutine:      A program which computes f(x) assuming x is in X-register upon entry

 STACK INPUTS OUTPUTS Z / f '''(x) Y h f ''(x) X x f '(x)

Example:      f(x) = exp(x) + ln(x)     x = 2

 01  LBL "T"   02  E^X   03  LASTX   04  LN   05  +   06  RTN

"T"  ASTO 00

-With  h = 0.1

0.1   ENTER^
2     XEQ "dF"  >>>>   f '(2) = 7.889056107 = R03                 ---Execution time = 13 seconds---
RDN   f "(2) = 7.139056635 = R04
RDN  f '''(2) = 7.639051753 = R05

Derivatives of a Function of 2 variables

-"dF2" calls "dF" and "MDV" to compute the 5 partial derivatives:

f 'x = f / x ; f 'y = f / y ;  f "xx = 2f / x2 ; f "yy = 2f / y2  and   f "xy  = 2f / xy

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

R01 = x                     R03 = f 'x = f / x            R06 =  f 'y = f / y           R09 =  f "xy  = 2f / xy
R02 = y                     R04 = f "xx = 2f / x2       R07 = f "yy = 2f / y2
R05 = f "'xxx = 3f / x3     R08 = f "'yyy = 3f / y3

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

 STACK INPUTS OUTPUTS T / f "yy = ¶2f / ¶y2 Z h f "xx = ¶2f / ¶x2 Y y f 'y = ¶f / ¶y X x f 'x = ¶f / ¶x

Example:     f(x) = exp(-x2).Ln(y)    x = 1 , y = 2

 01  LBL "T"  02  X^2  03  CHS  04  E^X  05  X<>Y  06  LN  07  *  08  RTN

T  ASTO 00

-If we choose h = 0.1

0.1   ENTER^
2     ENTER^
1     XEQ "dF2"  >>>>  f 'x  = f / x  =   -0.509989198     the exact result is    -0.509989195                ---Execution time = 57s---
RDN   f 'y  = f / y  =    0.183939720      the exact result is     0.183939721
RDN  f "xx  = 2f / x2 =  0.509989188      the exact result is     0.509989195
RDN  f "yy = 2f / y2 = -0.091969868      the exact result is    -0.091969860
LASTX  f "xy  = 2f / xy = -0.367879468    the exact result is    -0.367879441

-You have also:

R05 = f "'xxx = 3f / x3   = 1.019980483  exact result = 1.019978390
R08 = f "'yyy = 3f / y3 =  0.091970052   exact result = 0.091969860

Derivatives of a Function of 3 variables

-"dF3" calls  "dF" to compute the 6 partial derivatives  f 'x = f / x ; f 'y = f / y ; f 'z = f / z ;  f "xx = 2f / x2 ;  f "yy = 2f / y2 ;  f "zz = 2f / z2

and  f "'xxx = 3f / x3      f "'yyy = 3f / y3      f "'zzz = 3f / z3

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

R03 = f 'x = f / x            R06 =  f 'y = f / y           R09 = f 'z = f / z           R12 = f(x,y,z)
R04 = f "xx = 2f / x2       R07 = f "yy = 2f / y2       R10 = f "zz = 2f / z2      R13 = x  R14 = y  R15 = z
R05 = f "'xxx = 3f / x3     R08 = f "'yyy = 3f / y3     R11 = f "'zzz = 3f / z3

Flags: /
Subroutines:   "dF" and 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 T h Df Z z f 'z = ¶f / ¶z Y y f 'y = ¶f / ¶y X x f 'x = ¶f / ¶x

Df = Laplacian ( f ) =   2f / x2 + 2f / y2 + 2f / z2

Example:     f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 2 , z = 3

 01  LBL "T"  02  X^2  03  CHS  04  E^X  05  RDN  06  X^2  07  +  08  LN  09  R^  10  *  11  RTN

T  ASTO 00

-If we choose h = 0.1

0.1   ENTER^
3     ENTER^
2     ENTER^
1     XEQ "dF3"  >>>>   f 'x  = f / x  =   -1.413720683              R04 = f "xx = 2f / x2 =  1.413720682             ---Execution time = 52s---
RDN   f 'y  = f / y  =    0.210216822              R07 = f "yy = 2f / y2 = -0.015015424
RDN   f 'z  = f / z   =    0.052554204              R10 = f "zz2f / z2 = -0.007507569
RDN    2f / x2 + 2f / y2 + 2f / z2 =     1.409197689      and

R05 = f "'xxx  = 3f / x3  =  2.863447421
R08 = f "'yyy = 3f / y3  = -0.042900947
R11 = f "'zzz  = 3f / z3 =   0.002145569

Mixed Derivative of a Function of 2 variables

"MDV" calculates  2f / xy   with the formula:

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 ) ] + O(h10)

-This formula is exact for any polynomial of degree < 11               (   fkm =  f ( x + k.h , y + m.h )  )

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

R01 = x        R04 = f(x,y)
R02 = y        R05 = f "xy = 2f / xy     R06: temp
R03 = h
Flags: /
Subroutine:         A program which computes f(x,y) assuming x is in X-register and y is in Y-register upon entry.

 STACK INPUTS OUTPUTS Z h / Y y / X x f "xy = ¶2f / ¶x¶y

Example:     f(x) = exp(-x2).Ln(y)    x = 1 , y = 2

 01  LBL "T"   02  X^2   03  CHS   04  E^X   05  X<>Y   06  LN   07  *   08  RTN

T  ASTO 00

-If we choose h = 0.1

0.1   ENTER^
2     ENTER^
1     XEQ "MDV"  >>>>    f "xy  = 2f / xy = -0.367879468              ---Execution time = 28s---

-The exact result is    -0.367879441.....

3°) Biharmonic Operator

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

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

-The 4th derivative  d4f / dx4  may be approximated by the following formula, of order 10

( 192654 F - 140196 G + 52428 H - 9738 I +1261 J - 82 K ) / ( 15120 h4 )

where      G = f(x+h) + f(x-h)  ,  H = f(x+2h) + f(x-2h)  ,  I = f(x+3h) + f(x-3h)  ,  J = f(x+4h) + f(x-4h)  , K = f(x+5h) + f(x-5h)

and   F = f(x)

-Likewise,  4f / x2y2   may be computed by

4f / x2 y2 ~  [ 20881861 f00 /  3240000 - ( 5/3 ) ( 2 g1 - h1 ) +  ( 5/84 ) ( 2 g2 - h2 ) - ( 5/1134 ) ( 2 g3 - h3 )

+ ( 5/16128 ) ( 2 g4 - h4 ) - ( 1/78750 ) ( 2 g5 - h5 ) ] / h4

where  fij = f ( x+i.h , y+j.h )  ;  gk =  fk0 + f-k0 + f0k + f0-k  and   hk =  fkk + f-k-k + fk-k + f-kk

-Combining these formulas, we get 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            ( Register R00 is to be initialized before executing "BHRM" )

R01 = x        R04 = h
R02 = y        R05 = D2 f
R03 = z        R06-R07-R08-R09:  temp
Flags: /
Subroutine:         A program which computes f(x,y,z) assuming x is in X-register, y is in Y-register and z is in Z-register upon entry.

 01 LBL "BHRM"  02 STO 01  03 RDN  04 STO 02  05 RDN  06 STO 03  07 X<>Y  08 STO 04  09 STO 06  10 X<>Y  11 R^  12 R^  13 XEQ IND 00  14 6  15 ST* 06  16 *  17 STO 09    18 XEQ 01          19 24  20 *  21 RCL 07  22 5221  23 *  24 -  25 125  26 / 27 STO 05  28 XEQ 01  29 75  30 *  31 RCL 07  32 10388  33 *  34 -  35 16  36 /  37 ST- 05  38 XEQ 01  39 200  40 *  41 RCL 07          42 15407  43 *  44 -  45 3  46 /  47 ST+ 05  48 6  49 ST/ 05  50 XEQ 01  51 150  52 * 53 RCL 07  54 4969  55 *  56 -  57 ST- 05  58 7  59 ST/ 05  60 XEQ 01  61 600  62 *  63 RCL 07  64 4069  65 *  66 -  67 GTO 04  68 LBL 01  69 RCL 04          70 ST- 06  71 CLX  72 STO 07  73 STO 08  74 LBL 02  75 RCL 03  76 RCL 06  77 +  78 RCL 02 79 RCL 01  80 XEQ IND 00  81 ST+ 07  82 RCL 03  83 RCL 02  84 RCL 06  85 +  86 RCL 01  87 XEQ IND 00  88 ST+ 07  89 RCL 03  90 RCL 02  91 RCL 01  92 RCL 06  93 +  94 XEQ IND 00  95 ST+ 07  96 RCL 06  97 CHS  98 STO 06  99 X<0? 100 GTO 02 101 STO 10 102 LBL 03 103 RCL 03 104 RCL 06 105 + 106 RCL 02 107 RCL 10 108 + 109 RCL 01 110 XEQ IND 00 111 ST+ 08 112 RCL 03 113 RCL 06 114 + 115 RCL 02 116 RCL 01 117 RCL 10 118 + 119 XEQ IND 00 120 ST+ 08 121 RCL 03 122 RCL 02 123 RCL 06 124 + 125 RCL 01 126 RCL 10 127 + 128 XEQ IND 00 129 ST+ 08 130 RCL 10 131 CHS 132 STO 10 133 X<0? 134 GTO 03 135 RCL 06 136 CHS 137 STO 06 138 X<0? 139 GTO 03 140 RCL 08         141 RCL 09  142 ST- 07 143 ST+ X 144 - 145 RTN 146 LBL 04 147 RCL 05  148 + 149 180 150 RCL 04 151 X^2 152 X^2 153 * 154 / 155 STO 05 156 END

( 231 bytes / SIZE 010 )

 STACK INPUTS OUTPUTS T h / Z z / Y y / X x D2 f

Example:     f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 2 , z = 3

 01  LBL "T"  02  X^2  03  CHS  04  E^X  05  RDN  06  X^2  07  +  08  LN  09  R^  10  *  11  RTN

T  ASTO 00

-If you choose h = 0.1

0.1   ENTER^
3     ENTER^
2     ENTER^
1     XEQ "BHRM"  >>>> D2 f = -14.34282056 = R05             ---Execution time = 1m41s---

-The exact result is  -14.34264116  ( relative error about 2 E-5 )

-With h = 0.2  it yields   -14.34285066 and with h = 0.15  we get  -14.34265712

Notes:

-91 evaluations of the function are performed (!).

-Though the formula is of order 10 - due to roundoff-errors - 4th derivatives are difficult to get with a great precision.
-So, there are seldom more than 5 or 6 significant digits.

4°) Curl, Divergence, Gradient & Laplacian

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

-If the coordinates are rectangular, clear flags F01 & F02
-If the coordinates are cylindrical, set F01 and clear F02
-If the coordinates are spherical, set F02                                            ( if F02 & F01 are set, F01 is automatically cleared )

-In the last 2 cases, the HP41 sets the RAD mode automatically.

Formulae:

•  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 = h

R01 = x         R04-R05-R06 = Curl E       R08-R09-R10 = Grad (X)    R17 = Lap (X)
R02 = y         R07 = Div E                         R11-R12-R13 = Grad (Y)   R18 = Lap (Y)
R03 = z                                                     R14-R15-R16 = Grad (Z)    R19 = Lap (Z)

R23 to R38 contain the same results as R04 to R19

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 is in X-register, y is in Y-register and z is in Z-register upon entry.

XROM may be replaced by XEQ

 01 LBL "CDGL"  02 FC? 01  03 FS? 02  04 RAD  05 STO 13  06 RDN  07 STO 14  08 RDN  09 STO 15  10 X<>Y  11 STO 02   12 X<>Y  13 R^  14 R^  15 "X"  16 ASTO 00  17 XROM "dF3"  18 STO 26          19 STO 27         20 RDN  21 STO 28  22 CHS  23 STO 25  24 RDN  25 STO 29  26 STO 24  27 X<>Y  28 STO 36  29 FS? 02  30 GTO 02  31 FC? 01  32 GTO 00  33 XEQ 04  34 +  35 STO 36  36 RCL 12  37 RCL 13  38 ST/ 25 39 ST/ 28  40 /  41 ST+ 26  42 GTO 00  43 LBL 02  44 XEQ 03  45 STO 36  46 RCL 13  47 ST/ 25  48 ST/ 28  49 RCL 14   50 SIN  51 *  52 ST/ 24  53 ST/ 29  54 RCL 12  55 ST+ X  56 RCL 13         57 /  58 ST+ 26  59 LBL 00  60 "Y"  61 ASTO 00  62 RCL 02  63 RCL 15  64 RCL 14  65 RCL 13  66 XROM "dF3"  67 FC? 01  68 FS? 02  69 GTO 01  70 ST+ 25  71 STO 30  72 RDN  73 STO 31  74 ST+ 26  75 RDN  76 STO 32 77 CHS  78 STO 23  79 X<>Y  80 STO 37  81 GTO 00  82 LBL 01  83 FS? 02  84 GTO 02  85 XEQ 04  86 +  87 STO 37  88 RCL 03   89 STO 30  90 RCL 06  91 RCL 13  92 /  93 ST+ 26  94 STO 31         95 RCL 09  96 STO 32  97 CHS  98 STO 23  99 RCL 03 100 RCL 12 101 RCL 13 102 / 103 + 104 ST+ 25 105 GTO 00 106 LBL 02 107 XEQ 03 108 STO 37 109 RCL 03 110 STO 30 111 RCL 06 112 RCL 13 113 / 114 STO 31 115 RCL 09 116 LASTX 117 RCL 14 118 SIN 119 * 120 / 121 STO 32 122 RCL 12 123 RCL 14 124 TAN 125 / 126 RCL 06  127 + 128 RCL 13 129 / 130 ST+ 26 131 RCL 09        132 CHS 133 RCL 13 134 RCL 14 135 SIN 136 * 137 / 138 STO 23 139 RCL 12 140 RCL 13 141 / 142 RCL 03 143 + 144 ST+ 25 145 LBL 00 146 "Z" 147 ASTO 00 148 RCL 02 149 RCL 15 150 RCL 14 151 RCL 13 152 XROM "dF3" 153 FC? 01 154 FS? 02 155 GTO 01 156 STO 33 157 ST- 24 158 RDN 159 STO 34 160 ST+ 23 161 RDN 162 STO 35 163 ST+ 26 164 X<>Y 165 STO 38        166 GTO 00 167 LBL 01 168 FS? 02 169 GTO 02 170 XEQ 04 171 + 172 STO 38 173 RCL 03 174 STO 33 175 RCL 06 176 RCL 13 177 / 178 ST+ 23 179 STO 34 180 RCL 09 181 STO 35 182 ST+ 26 183 RCL 03 184 ST- 24 185 GTO 00 186 LBL 02 187 XEQ 03 188 STO 38 189 RCL 03 190 STO 33 191 RCL 06 192 RCL 13 193 / 194 STO 34 195 RCL 09 196 LASTX 197 RCL 14 198 SIN 199 * 200 / 201 STO 35 202 ST+ 26 203 RCL 06  204 RCL 12        205 RCL 14 206 TAN 207 / 208 + 209 RCL 13 210 / 211 ST+ 23 212 RCL 03 213 RCL 12 214 RCL 13 215 / 216 + 217 ST- 24 218 GTO 00 219 LBL 03 220 XEQ 04 221 RCL 14 222 SIN 223 X^2 224 / 225 RCL 06 226 RCL 14 227 TAN 228 / 229 + 230 RCL 13 231 X^2 232 / 233 + 234 RTN 235 LBL 04 236 RCL 07 237 RCL 13 238 / 239 RCL 03  240 FS? 02 241 ST+ X 242 + 243 RCL 13        244 / 245 RCL 04 246 + 247 RCL 10 248 RTN 249 LBL 00 250 RCL 02 251 STO 00 252 RCL 13 253 STO 01 254 RCL 14 255 STO 02 256 RCL 15 257 STO 03 258 23.004016 259 REGMOVE 260 4.019 261 SIGN 262 RCL 07 263 RCL 06 264 RCL 05 265 RCL 04 266 END

( 411 bytes / SIZE 039 )

 STACK INPUTS OUTPUTS T h Div E Z x3 Curl3 E Y x2 Curl2 E X x1 Curl1 E L / 4.019

4.019 = control number of all the results.

Example1 - Rectangular Coordinates:   CF 01  CF 02

E = ( f , g , h )  = [ exp(-x2) Ln(y2+z) , x2y2z2 , exp(x) y2z ]     With  x = 1 , y = 2 , z = 3

-Load the short routines:

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

CF 01   CF 02

-If you choose  h = 0.1

0.1  ENTER^
3    ENTER^
2    ENTER^
1    XEQ "CDGL"  >>>>     8.619381890  = R04                                                                    ---Execution time = 127s---
RDN    -32.56682777  = R05
RDN     71.78978318  = R06
RDN     45.44140669  = R07

-So,  Curl E = ( 8.619381890 ,  -32.56682777 ,  71.78978318 )

and    Div E = 45.44140669

-You also find in registers R08 to R19:

Grad f = ( -1.431720683 , 0.210216822 , 0.052554204 )
Grad g = ( 72 , 36 , 24 )
Grad h = ( 32.61938197 , 32.61938189 , 10.87312737 )

Lapl f = 1.409197689
Lapl g = 98
Lapl h = 48.92907013

Example2 - Cylindrical Coordinates:   SF 01  CF 02

E = ( f , g , h )  = (  r z2 sin2f , r2 z , r3 z cos f )                   r = 2 , f = PI/5 , z = 1

 01  LBL "X"  02  RDN  03  SIN  04  *  05  X^2  06  R^  07  *  08  RTN  09  LBL "Y"  10  X^2  11  RCL Z  12  *  13  RTN  14  LBL "Z"  15  X^2  16  LASTX  17  *  18  X<>Y  19  COS  20  *  21  *  22  RTN

SF 01   CF 02

-If you choose  h = 0.1

0.1   ENTER^
1
PI  5  /
2    R/S   >>>>    -6.351141008  = R04                                                                     ---Execution time = 135s---
RDN    -8.326237924  = R05
RDN     5.048943483  = R06
RDN     7.163118958  = R07

-So,  Curl E = ( -6.351141008 ,  -8.326237924 ,  5.048943483 )

and    Div E = 7.163118958

-You also find in registers R08 to R19:

Grad f = ( 0.345491503 , 0.951056518 , 1.381966009 )
Grad g = ( 4 , 0 , 4 )
Grad h = ( 9.708203933 , -2.351141008 , 6.472135952 )

Lapl f = 1.863728736
Lapl g = 4
Lapl h = 12.94427209

Example2 - Spherical Coordinates:    CF 01  SF 02

E = ( f , g , h )  = (  r sin2 q cos2 f , r2sin f  , r3 cos q cos2 f )                   r = 2 , q = PI/3 , f = PI/5

 01  LBL "X"  02  RDN  03  SIN  04  X<>Y  05  COS  06  *  07  X^2  08  R^  09  *  10  RTN  11  LBL "Y"  12  X^2  13  RCL Z  14  SIN  15  *  16  RTN  17  LBL "Z"  18  X^2  19  LASTX  20  *  21  X<>Y  22  COS  23  *  24  X<>Y  25  COS  26  X^2  27  *  28  RTN

CF 01    SF 02

-If you choose  h = 0.1

0.1
PI  5  /
PI  3  /
2    R/S   >>>>   -3.379867347  = R04                                                                     ---Execution time = 187s---
RDN    -6.059707086  = R05
RDN     2.959890531  = R06
RDN    -0.045010875  = R08

-So,  Curl E = ( -3.379867347 ,  -6.059707086 ,  2.959890531 )

and    Div E = -0.045010875

-You also find in registers R08 to R19:

Grad f = ( 0.490881375 , 0.566820985 , -0.823639103 )
Grad g = ( 2.351141010 , 0 , 1.868344717 )
Grad h = ( 3.927050988 , -2.267283945 , -2.196370944 )

Lapl f = 0.018237234
Lapl g = 2.742997604
Lapl h = 5.721039298

5°) Laplacian of a Function of N variables ( n < 10 )

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

"SD" is called as a subroutine( see below ), so the same formulas are used, of order 10 if h > 0 or of order 4 if h < 0

>>> 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 Rnn are to be initialized before executing "LAPL"  )

•  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn

R10 = h
R12-..........-R20: temp
Flag:  F10
Subroutines:   "SD" and a program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn

 01 LBL "FD"  02 LBL 01  03 STO 13  04 CLX  05 STO 14  06 RCL IND 13  07 STO 11  08 X<> Z  09 STO 10  10 FC? 03  11 XEQ IND 00  12 SIGN  13 X=0?  14 RTN  15 RCL 10  16 X<0?  17 GTO 00  18 XEQ 03  19 14  20 *  21 STO 12  22 XEQ 03  23 4  24 *  25 ST- 12  26 XEQ 03  27 ST+ 12  28 6  29 ST* 12  30 XEQ 03 31 ST- 12  32 25  33 ST* 12  34 XEQ 03  35 ST+ X  36 ST+ 12  37 2520  38 LBL 13  39 RCL 10  40 *  41 ST/ 12  42 RCL 11  43 STO IND 13   44 RCL 12  45 RTN  46 GTO 01  47 LBL 00  48 XEQ 03  49 8  50 *  51 STO 12  52 XEQ 03  53 ST- 12  54 12  55 GTO 13  56 LBL 03  57 RCL 10  58 ST+ 14  59 RCL 11  60 RCL 14 61 -  62 STO IND 13  63 XEQ IND 00  64 STO 15  65 RCL 11  66 RCL 14  67 +  68 STO IND 13  69 XEQ IND 00  70 RCL 15  71 -  72 RTN  73 LBL "SD"  74 LBL 02  75 CF 10  76 X=Y?  77 SF 10  78 STO 16  79 RDN  80 STO 17  81 CLX  82 STO 15  83 RCL IND T  84 STO 13  85 RCL IND 17  86 STO 14  87 R^  88 STO 10  89 XEQ IND 00  90 SIGN 91 X=0?  92 RTN  93 LASTX  94 STO 11  95 RCL 10  96 X<0?  97 GTO 00  98 XEQ 03  99 42 100 * 101 STO 12 102 XEQ 03 103 6 104 * 105 ST- 12 106 XEQ 03 107 ST+ 12 108 XEQ 03 109 8 110 / 111 ST- 12 112  E3 113 ST* 12 114 XEQ 03 115 8 116 * 117 ST+ 12 118 RCL 13 119 STO IND 16  120 25200 121 LBL 14 122 FC? 10 123 ST+ X 124 FC?C 10 125 CHS 126 RCL 10 127 X^2 128 * 129 ST/ 12 130 RCL 12 131 RTN 132 GTO 02 133 LBL 00 134 XEQ 03 135 16 136 * 137 STO 12 138 XEQ 03 139 ST- 12 140 RCL 13 141 STO IND 16  142 12 143 GTO 14 144 LBL "LAPL" 145 LBL 04 146 STO 19 147 X<>Y 148 STO 10 149 CLX 150 STO 20 151 LBL 05 152 RCL 10 153 RCL 19 154 ENTER 155 XROM "SD" 156 ST+ 20 157 DSE 19 158 GTO 05 159 RCL 20 160 RTN 161 GTO 04 162 LBL 03 163 RCL 10 164 ST+ 15 165 RCL 15 166 RCL 13 167 + 168 STO IND 16  169 XEQ IND 00 170 STO 18 171 RCL 13 172 RCL 15 173 - 174 STO IND 16 175 XEQ IND 00 176 RCL 11 177 ST+ X 178 - 179 ST+ 18 180 RCL 18 181 FS? 10 182 RTN 183 RCL 14 184 RCL 15 185 - 186 STO IND 17 187 XEQ IND 00 188 ST- 18 189 RCL 13 190 STO IND 16 191 XEQ IND 00 192 ST+ 18 193 RCL 14 194 RCL 15 195 + 196 STO IND 17  197 XEQ IND 00 198 ST+ 18 199 RCL 15 200 ST+ IND 16 201 XEQ IND 00 202 ST- 18 203 RCL 14 204 STO IND 17 205 RCL 18 206 RTN 207 END

( 345 bytes / SIZE 021 )

 STACK INPUTS OUTPUTS Y h / X n Laplacian(f)

where n is the number of variables

Example:     f(x) = exp(-x2.t).Ln(y2+z)    x = 1 , y = 1 , z = 1 , t = 1

LBL "T"       RCL 04      X^2           *                     "T"   ASTO 00
RCL 01       *                 RCL 03     RTN                 1    STO 01  STO 02  STO 03  STO 04
X^2             E^X            +
CHS            RCL 02      LN

-With h = 0.1

0.1   ENTER^
4     XEQ "LAPL"  >>>>    Laplacian(f) =  2f / x2 + 2f / y2 + 2f / z2 + 2f / t2 = 0.673013891                  ---Execution time = 54s---

( the last 3 decimals should be 932 )

-With  h < 0  we get more quickly:

-0.01  ENTER^
4        R/S        >>>>   Lapl(f) = 0.67300675                  ---Execution time = 25s---

First Derivatives of a Function of N variables

-If h > 0  "FD"  employs a formula of order 10
-If h < 0  "FD"  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 if CF 03

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

•  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn

R10 = h      R11: temp           R13 = i         R15 = xi
R12 =  f / xi       R14 = 0 , h , 2h , ....     R16: temp
Flag:  F03
Subroutine:   A program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn

 STACK INPUTS OUTPUTS Y h / X i f 'xi = ¶f / ¶xi

Example:     f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 1 , z = 1

LBL "T"        E^X               +                    "T"    ASTO 00
RCL 01        RCL 02          LN                  1     STO 01   STO 02   STO 03
X^2              X^2              *
CHS            RCL 03          RTN

-With  h = 0.1

0.1   ENTER^
1    XEQ "FD"   >>>>   f 'x1 = f / x1 = -0.509989198   ( The last decimal should be a 5 )               ---Execution time = 10s---

0.1   ENTER^
2       R/S     >>>>   f 'x2 = f / x2 = 0.367879440     ( The last decimal should be a 1 )                    ---Execution time = 11s---

0.1   ENTER^
3       R/S      >>>>   f 'x3 = f / x3 = 0.183939720     ( The last decimal should be a 1 )                    ---Execution time = 11s---

Notes:

-This program uses the same formula of order 10 as "dF"  if  h > 0
-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 Riemanian section.

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

-For instance, with h = -0.01

0.01  CHS  ENTER^
1           R/S          >>>>    f 'x1 = f / x1 = -0.509989196                     ---Execution time = 4.5s---

Second Derivatives of a Function of N variables

-If h > 0  "SD"  employs a formula of order 10
-If h < 0  "SD"  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 ) ]

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

>>> In both cases if the function places an alpha data in X-register, 0 is returned without calculating the whole formula

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

•  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn

R10 = h        R11: temp        R13 = xi           R15 = 0 , h , 2h , ....
R12 =  2f / xixj        R14 = xj            R16 =  i    R17 = j   R18: temp
Flag:  F10
Subroutine:   A program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn

 STACK INPUTS OUTPUTS Z h / Y j / X i f "ij = ¶2f / ¶xi¶xj

Example:       f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 1 , z = 1

LBL "T"        E^X               +                    "T"    ASTO 00
RCL 01        RCL 02          LN                  1     STO 01   STO 02   STO 03
X^2              X^2               *
CHS            RCL 03          RTN

-With  h = 0.1

0.1  ENTER^
1    ENTER^
XEQ "SD"  >>>>   f "xx = 2f / x2 =0.509989188   ( exact derivative = 0.509989195 )                   ---Execution time = 12s---

0.1   ENTER^
1    ENTER^
2       R/S     >>>>   f "xy = 2f / xy = -0.735758883  ( the last decimal should be a 2 )                         ---Execution time = 32s---

Notes:

-Like "FD" 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 = 0.509987083                   ---Execution time = 5.6s---

-0.01   ENTER^
1      ENTER^
2         R/S     >>>>   f "xy = 2f / xy = -0.735758125                   ---Execution time = 13.8s---

6°) Linear Systems

-We need a program to find the inverse of the matrix [ gij ]  i-e  [ gij ]
-But "LS" may also be used for another purpose

"LS" allows you to solve linear systems, including overdetermined and underdetermined systems.
You can also invert up to a 12x12 matrix.

The determinant of this matrix is also computed and stored in register R00.
( if there are more rows than columns, R00 is not always equal to the determinant of the upper left matrix because of rows exchanges )

Gaussian elimination with partial pivoting is used.

One advantage of this program is that you can choose the beginning data register - except R00 -
You store the first coefficient into Rbb , ... , up to the last one into Ree ( column by column )             ( bb > 00 )

Then you key in  bbb.eeerr  ENTER^
p          XQ "LS"    and the system will be solved.    ( bbb.eeerr ends up in L-register )

where r is the number of rows of the matrix
and  p  is a small non-negative number that you choose to determine tiny elements: if a pivot is smaller than p it will be considered to be zero.

Don't interrupt "LS" because registers P and Q are used ( there is no risk of crash, but their contents could be altered )

 01 LBL "LS"  02 STO M  03 X<>Y  04 ENTER  05 INT  06 LASTX  07 FRC  08 ISG X  09 INT  10 STO O         11 RCL Y   12 +  13 DSE X  14  E3  15 /  16 +  17 X<> O  18 .1  19 %  20 +  21 STO Q   22 SIGN  23 STO 00  24 X<>Y  25 STO P 26 LBL 01  27 STO N  28 INT  29 RCL O   30 FRC  31 +  32 ENTER  33 ENTER  34 CLX  35 LBL 02   36 RCL IND Z  37 ABS  38 X>Y?  39 STO Z   40 X>Y?  41 +  42 RDN  43 ISG Z  44 GTO 02  45 RCL N  46 ENTER  47 FRC  48 R^  49 INT  50 + 51 X=Y?  52 GTO 04  53 LBL 03  54 RCL IND X  55 X<> IND Z  56 STO IND Y  57 ISG Y  58 RDN  59 ISG Y  60 GTO 03  61 RCL 00  62 CHS  63 STO 00  64 LBL 04   65 RCL M   66 RCL IND N  67 ST* 00  68 ABS  69 X<=Y?  70 GTO 09  71 RCL N  72 LASTX  73 LBL 05  74 ST/ IND Y  75 ISG Y 76 GTO 05  77 LBL 06  78 RCL N  79 ENTER  80 FRC  81 RCL O   82 INT  83 +  84 X=Y?  85 GTO 08  86 RCL IND X  87 SIGN  88 RDN  89 LBL 07   90 RCL IND Y  91 LASTX  92 *  93 ST- IND Y  94 ISG Y  95 RDN  96 ISG Y  97 GTO 07  98 LBL 08  99 ISG O 100 GTO 06 101 RCL Q  102 INT 103 ST- O 104 LBL 09  105 RCL Q  106 ST+ O 107 X^2 108 RCL P        109 INT 110 ENTER 111 SIGN 112 ST+ N  113 - 114 + 115 RCL N  116 X>Y? 117 CLX 118 ISG X 119 GTO 01 120 X<> P 121 SIGN 122 RCL 00 123 CLA 124 END

( 185 bytes / SIZE var. )

 STACK INPUTS OUTPUTS Y bbb.eeerr 1 X p determinant L / bbb.eeerr

Example1:     Find the solution of the system:

2x + 3y + 5z + 4t = 39
-4x + 2y + z + 3t  = 15
3x  -  y + 2z + 3t = 19
5x + 7y - 3z + 2t = 18

If you choose to store the first element into R01, you have to store these 20 numbers:

2  3  5  4  39                  R01  R05  R09  R13  R17
-4  2  1  3  15        in        R02  R06  R10  R14  R18     respectively
3 -1  2  3  19                  R03  R07  R11  R15  R19
5  7 -3  2  18                  R04  R08  R12  R16  R20

The first register is R01, the last register is R20 and there are 4 rows, therefore the control number of the matrix is  1.02004
If you choose p = 10-7  key in:

1.02004 ENTER^
E-7    XEQ "LS"   and  you obtain  840 in X-register and in R00:  it is the determinant of the 4x4 matrix on the left.           ---Execution time = 31s---

Registers R01 thru R16 now contains the unit matrix and registers R17 thru R20 contains the solution  x = 1 ; y = 2 ; z = 3 ; t = 4
Thus, the array has been changed into:

1  0  0 0  1
0  1  0 0  2              the solution is the last column
0  0  1  0  3              of the new matrix.
0  0  0  1  4

Example2:    Solve the system:

5x + y + z  = 8
4x - y + 2z = 13
x + 2y - z  = -5
7x - 4y + 5z = 31

-Suppose we need to store the first element into  R11 , then we store these 16 numbers

5   1  1   8                       R11  R15  R19  R23
4  -1  2  13         in          R12  R16  R20  R24          respectively
1   2  -1  -5                     R13  R17  R21  R25
7  -4  5  31                      R14  R18  R22  R26

Here, the control number of this array is   11.02604   and if we take p = 10-7 again,

11.02604  ENTER^
E-7     XEQ "LS"   gives   -5.4  10-17  in X-register and in R00 = the determinant of the 4x4 matrix, which is of course 0   ---Execution time = 21s--

and registers R11 thru R27 are now:

1  0   0.333333333     2.333333333
0  1  -0.666666667   -3.666666669
0  0      5.10-10              -2.10-9
0  0          0                    4.10-9

Thus, the system is equivalent to:

x + z /3  =  7/3              and there is an infinite set of solutions:
y - 2z /3 = -11/3           z may be chosen at random and x and y are determined by
0  =  0                 x = -z /3 + 7/3
0  =  0                 y = 2z /3 - 11/3

-If the last number of the initial matrix is 32 instead of 31 the array is changed into:

1  0   0.333333333   0        i-e         x +  z /3 = 0
0  1  -0.666666667   0                     y - 2z /3 = 0
0  0        5.10-10        0                                 0 = 0
0  0           0              1                                 0 = 1        and there is no solution at all !

Example3:    To invert the matrix

2  3  4                                                                             R11  R14  R17
3  5  6      store these coefficients - for example - into       R12  R15  R18
4  6  7                                                                             R13  R16  R19

1  0  0            R20   R23  R26
Likewise, store the unit matrix   0  1  0   into    R21   R24  R27
0  0  1            R22   R25  R28

-The control number is  11.02803  and if you choose  p = 0

11.02803  ENTER^
0        XEQ "LS"   >>>>   det = -1 = R00                                           ---Execution time = 21s---

-The array is now

R11  R14  R17  R20  R23  R26         1    0    0    1   -3   2
R12  R15  R18  R21  R24  R27   =    0    1    0   -3   2    0
R13  R16  R19  R22  R25  R28         0    0    1    2    0   -1

-So the inverse matrix is

1    -3     2
-3     2     0
2     0    -1

Notes:

-Thus, this program can be used to invert an nxn matrix, but it requires 2n2 registers and so, it cannot invert a 13x13 matrix.

7°)   Riemannian Geometry

a) Initialization - Metric Tensor & Christoffel Symbols

-Always execute "INIGC" 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.

-However, if the metric is diagonal ( i-e orthogonal ), only create  LBL "G11"  LBL "G22" ......  provided you set flag  F00

-They must take  x1 , x2 , ........... , xn  in registers R01  R02 .............  Rnn   ( n < 7 ) and return  gij  in X-register

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

-The determinant g = det gij is stored in R40 and n(n+1)/2 is sored in R39

-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 "INIGC" , do not disturb register R09 , R39 & R41 to R40+n(n+1)(n+2)/2

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

•  R01 = x1                    R09 = n                              R41 to R40+n(n+1)/2                     =   gij   with  i <= j
•  R02 = x2                    R10 = h                               R41+n(n+1)/2 to R40+n(n+1)         =  gij    with i <= j
.............                     R39 =  n(n+1)/2                   R41+n(n+1) to R40+n(n+1)(n+2)/2 = Gkij with  i <= j
•  Rnn = xn                    R40 = g = det gij

Flag:   F04                       CF 04 = The HP41 calculates gij , gij  and Gkij
SF 04 =  Only gij and gij  are computed

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

 01 LBL "INIGC"  02 STO 09  03 STO 11  04 ENTER  05 X^2  06 +  07 2  08 /  09 STO 39   10 X<>Y  11 STO 10   12 CLX  13 40  14 +  15 STO 13         16 RCL 09  17 X^2  18 +  19 STO 26  20 RCL 09  21  E5  22 /  23 +  24 STO 14  25 CLX  26 STO 15  27 LBL 04  28 RCL 11  29 STO 12  30 LBL 05 31 RCL 11  32 RCL 12  33 Fµ  34 ASTO X  35 XEQ IND X  36 SIGN  37 X#0?  38 LASTX  39 STO IND 13  40 STO IND 14  41 STO IND 26  42 DSE 26   43 DSE 14  44 DSE 13  45 DSE 12  46 GTO 05  47 CLX  48 SIGN  49 ST+ 15  50 RCL 15  51 ST- 26  52 RCL 14   53 FRC  54 RCL 26  55 +  56 STO 14  57 DSE 11  58 GTO 04  59 RCL 09  60 X^2 61 ST+ X  62 RCL 39  63 +  64 40  65 +  66 .1  67 %  68 +  69 RCL 09  70 X^2  71 DSE X  72 -  73 CLRGX  74 RCL 09   75 1  76 +  77  E5  78 /  79 +  80 SIGN  81 LBL 06  82 STO IND L   83 ISG L  84 GTO 06  85 LASTX  86 FRC  87 41  88 +  89 RCL 39  90 + 91  E-5  92 -  93 0  94 XROM "LS"  95 STO 40  96 1  97 STO 12  98 RCL 09  99 STO 11  100 ST* X 101 LASTX 102 INT 103 ST+ Y 104 LBL 07 105 RCL IND Y 106 STO IND Y 107 CLX 108 SIGN 109 ST+ Z 110 + 111 DSE Z 112 GTO 07 113 ISG 12 114 CLX 115 X<>Y 116 RCL 11 117 + 118 DSE X 119 RCL 12 120 X<> Z 121 DSE 11 122 GTO 07 123 RCL 39 124 RCL 09 125 FS? 04 126 CLX 127 STO 23  128 2 129 + 130 * 131 RCL 13  132 + 133 STO 26        134 STO 29 135 FS? 04 136 GTO 14 137 LBL 10 138 RCL 09 139 STO 21 140 LBL 01 141 RCL 21 142 STO 22 143 LBL 02 144 CLX 145 STO 27 146 RCL 09 147 STO 24 148 LBL 03 149 RCL 23 150 RCL 24 151 ENTER 152 CLX 153 CIJK 154 X=0? 155 GTO 00 156 STO 28 157 CLX 158 STO 25  159 RCL 22 160 RCL 24  161 Fµ 162 ASTO 00 163 RCL 10 164 RCL 21 165 XROM "FD" 166 ST+ 25 167 RCL 21 168 RCL 24 169 Fµ 170 ASTO 00 171 RCL 10 172 RCL 22 173 XROM "FD" 174 ST+ 25 175 RCL 21 176 RCL 22 177 Fµ 178 ASTO 00 179 RCL 10 180 RCL 24 181 XROM "FD" 182 ST- 25 183 RCL 25 184 RCL 28 185 * 186 ST+ 27 187 LBL 00 188 DSE 24 189 GTO 03 190 RCL 27  191 2 192 / 193 STO IND 26 194 DSE 26 195 DSE 22 196 GTO 02 197 DSE 21 198 GTO 01 199 DSE 23 200 GTO 10 201 LBL 14 202 RCL 29 203  E3 204 / 205 41 206 + 207 END

( 322 bytes / SIZE 41+n(n+1)(n+2)/2 )

 STACK INPUTS OUTPUTS Y h / X n bbb.eee (g&G)

where n is the number of variables  n < 8

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

g11 = 1 + x2 y2z2
g22 = 1 + x2 + y2 + z2                       and      gij = 0   if  i # j
g33 = 1 + x2 y2 + x2 z2 + y2 z2

-Load for instance these 6 routines in main memory:

 01  LBL "G11"  02  CLX  03  SIGN  04  RCL 01  05  RCL 02  06  RCL 03  07  *  08  *  09  X^2  10  + 11  RTN  12  LBL "G22"  13  CLX   14  SIGN  15  RCL 01  16  X^2  17  +  18  RCL 02  19  X^2  20  + 21  RCL 03  22  X^2  23  +  24  RTN  25  LBL "G33"  26  CLX  27  SIGN  28  RCL 01      29  RCL 02  30  * 31  X^2  32  +  33  RCL 01   34  RCL 03      35  *  36  X^2  37  +  38  RCL 02  39  RCL 03  40  * 41  X^2  42  +  43  RTN  44  LBL "G12"  45  LBL "G13"  46  LBL "G23"  47  ASTO X  48  RTN  49  END

-At the point  x = y = z = 1

1  STO 01  STO 02  STO 03

-We have n = 3 and if you choose h = 0.1

SF 00

0.1   ENTER^
3     XEQ "INIGC"   >>>>   41.070                                               ---Execution time = 4m37s---

-And you find in registers R41 thru R52

R41 = g11 = 2     R42 = g12 = 0      R44 = g13 = 0              R47 = g11 = 1/2     R48 =  g12 = 0            R50 = g13 =  0
R43 = g22 = 4     R45 = g23 = 0                                             R49 = g22 = 1/4          R51 = g23 = 0
R46 = g33 = 4                                                                               R52 = g33 = 1/4

-And in registers R53 to R70

R53 = G111 = 1/2    R54 = G112 =  1/2     R56 = G113 =  1/2
R55 = G122 = -1/2     R57 = G123 =   0
R58 = G133 =  -1

R59 = G211 = -1/4     R60 = G212 = 1/4     R62 = G213 =  0
R61 = G222 = 1/4      R63 = G223 = 1/4
R64 = G233 = -1/2

R65 = G311 = -1/4      R66 = G312 =    0        R68 = G313 =  1/2
R67 = G322 =  -1/4      R69 = G323 =  1/2
R70 = G333 =  1/2

-We have also   R39 = 6  i-e  n(n+1)/2    and   R40 = 32 = g = det gij

Notes:

-Since all the gij functions are polynomials of degree < 5, we could get the same results with derivative-formulas of order 4
-So, if you choose h = -0.01, you get the same values in 2m57s

-Execute a SIZE at least

 n SIZE 2 053 3 071 4 101 5 146 6 209 7 293

-So, there are not enough data registers for n = 8 if CF 04.

b) Recalling Gij & Cijk

-CIJK is an M-Code routine that recalls  Gkij provided  R39 = n(n+1)/2  and  R09 = n
-These 2 registers have been initialized by "INIGC"

-CIJK checks that the data registers do exist and don't contain alpha data
-However, it does not check for alpha data in the stack.

-The formula to get the address of the register Rmm that contains Gkij  is  ( with  i < j )

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

08B  "K"
00A  "J"
009   "I"
003   "C"
378   C=c
03C  RCR 3
106  A=C S&X
130  LDI S&X
027  027h=39d
146  A=A+C S&X
130  LDI S&X
200  200h
306  ?A<C S&X
381  GOTO
00A  NEXIST
0A6  A<>C S&X
270  RAMSLCT
106  A=C S&X
038  READATA
0EE  B<>C ALL
130  LDI S&X
01E  1Eh=30d
246  C=A-C S&X
270  RAMSLCT
038  READATA
10E  A=C ALL
046  C=0 S&X
270  RAMSLCT
0EE B<>C ALL
355   chk A&C
050   and SETDEC
268  Q=C
0AE  A<>C ALL
128  L=C
10E  A=C ALL
0F8  C=X
135  C=
060  A*C
138  C=L
025  C=
060  AB*C
128  L=C
078  C=Z
10E  A=C ALL
0B8  C=Y
30E  ?A<C ALL
027  JC+04
068  Z=C
0AE  A<>C ALL
0A8  Y=C
10E   A=C ALL
135  C=
060  A*C
0B8  C=Y
2BE  C=-C
025  C=
061  AB+C
04E   C
35C   =
090   2
269   C=
060  AB/C
138  C=L
025  C=
060  AB+C
078  C=Z
025  C=
060  AB+C
028  T=C
260  SETHEX
38D  C
008  ->S&X
106  A=C S&X
130  LDI S&X
028  028h=040d
146  A=A+C S&X
378  C=c
03C  RCR 3
146  A=A+C S&X
130  LDI S&X
200  200h
306  ?A<C S&X
381  GOTO
00A  NEXIST
0A6  A<>C S&X
270  RAMSLCT
038  READATA
10E  A=C ALL
046  C=0 S&X
270  RAMSLCT
0AE  A<>C ALL
0E8  X=C
3E0  RTN

( 93 words )

 STACK INPUTS OUTPUTS T / address-40 Z i inf(i,j) Y j sup(i,j) X k Gkij

Example1:

-With the metric above and after executing "INIGC"

3  ENTER^
ENTER^
2  XEQ "CIJK"  >>>> G233 = -1/2

-Register T also contains 24 , so the data register which contains G233  is R64

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

-For instance,  3  ENTER^  ENTER^  0  gives  g33 = 1/4

and   2  ENTER^   ENTER^   1  CHS   yields   g22 = 4

Notes:

-These routines check that the data registers  R09 & R39 do exist and that they do contain numbers
-However, they do not check  X-Y-Z-registers  for alpha data.

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

c)  Curvature Tensor

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

Rijkl = (1/2) ( jk gil + il gjk - jl gik - kik gjl ) + gmp ( GmjkGpil - GmjlGpik )

-The following listing includes "RIJKL"  "RIJ"  and  "RIJK^L"

 01 LBL "RIJKL"  02 LBL 14  03 STO 24  04 RDN  05 STO 23  06 RDN  07 STO 22  08 X<>Y  09 STO 21  10 FS? 49  11 OFF  12 RCL 09  13 STO 25  14 CLX  15 STO 28             16 LBL 11  17 RCL 09   18 STO 26  19 LBL 13  20 RCL 25  21 RCL 26  22 ENTER  23 SIGN  24 CHS  25 CIJK  26 X=0?  27 GTO 00 28 STO 29  29 RCL 22  30 RCL 23  31 RCL 25  32 CIJK  33 STO 27  34 RCL 21  35 RCL 24  36 RCL 26  37 CIJK  38 ST* 27  39 RCL 22  40 RCL 24  41 RCL 25  42 CIJK  43 STO 19             44 RCL 21  45 RCL 23  46 RCL 26  47 CIJK  48 ST* 19  49 RCL 27  50 RCL 19  51 -  52 RCL 29  53 *  54 ST+ 28 55 LBL 00  56 DSE 26  57 GTO 13  58 DSE 25  59 GTO 11  60 RCL 21  61 RCL 24  62 Fµ  63 ASTO 00  64 RCL 10  65 RCL 22   66 RCL 23  67 XROM "SD"  68 ST+ 25  69 RCL 22   70 RCL 23             71 Fµ  72 ASTO 00  73 RCL 10  74 RCL 21  75 RCL 24  76 XROM "SD"  77 ST+ 25  78 RCL 21  79 RCL 23  80 Fµ  81 ASTO 00 82 RCL 10  83 RCL 22  84 RCL 24  85 XROM "SD"  86 ST- 25  87 RCL 22  88 RCL 24  89 Fµ  90 ASTO 00  91 RCL 10  92 RCL 21  93 RCL 23  94 XROM "SD"  95 ST- 25  96 RCL 25             97 2  98 /  99 RCL 28 100 + 101 SIGN 102 RCL 21 103 RCL 22 104 RCL 23 105 RCL 24 106 X<> L 107 RTN 108 GTO 14 109 LBL "RIJ" 110 LBL 10 111 STO 23 112 X<>Y 113 STO 21 114 RCL 09 115 STO 22 116 CLX 117 STO 20 118 LBL 01 119 RCL 09 120 STO 24 121 LBL 02 122 RCL 22 123 RCL 24           124 ENTER 125 CLX 126 CIJK 127 X=0? 128 GTO 00 129 STO 30 130 RCL 21 131 RCL 22 132 RCL 23 133 RCL 24 134XROM "RIJKL" 135 RCL 30 136 * 137 ST+ 20 138 LBL 00 139 DSE 24 140 GTO 02 141 DSE 22 142 GTO 01 143 RCL 23 144 SIGN 145 RCL 21 146 RCL 20 147 RTN 148 GTO 10 149 LBL "RIJK^L" 150 LBL 12  151 STO 30            152 RDN 153 STO 24 154 RDN 155 STO 23 156 X<>Y 157 STO 22 158 RCL 09 159 STO 21 160 CLX 161 STO 20 162 LBL 03 163 RCL 21 164 RCL 30 165 ENTER 166 CLX 167 CIJK 168 X=0? 169 GTO 00 170 STO 31 171 RCL 21 172 RCL 22 173 RCL 23 174 RCL 24 175 XROM "RIJKL" 176 RCL 31             177 * 178 ST+ 20 179 LBL 00 180 DSE 21 181 GTO 03 182 RCL 30 183 SIGN 184 RCL 22 185 RCL 23 186 RCL 24 187 RCL 20 188 RTN 189 GTO 12 190 END

( 348 bytes / SIZE 41+n(n+1)(n+2)/2 )

 STACK INPUTS OUTPUTS T i i Z j j Y k k X l Rijkl L / l

Examples:    The same diagonal metric and with x = y = z = 1

1  ENTER^
2  ENTER^
1  ENTER^
2  XEQ "RIJKL"  >>>>  R1212 = -3/4                                     ---Execution time = 28s---

-Likewise

2  ENTER^
3  ENTER^
1  ENTER^
3     R/S      >>>>  R2313 = 1/2                                      ---Execution time = 38s---   (  or  20s  if  h < 0 )

Notes:

-The curvature tensor has many components but it has 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

1-time contravariant 3-times covariant curvature tensor

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

Rijkl = gim Rjmkl

-"RIJK^L" calculates these components

 STACK INPUTS OUTPUTS T i i Z j j Y k k X l Rlijk L / l

Examples:    The same diagonal metric and with x = y = z = 1

2  ENTER^
1  ENTER^
2  ENTER^
1  XEQ "RIJK^L"  >>>>  R1212 = 3/8                                      ---Execution time = 30s---

-Likewise

3  ENTER^
1  ENTER^
3  ENTER^
2     R/S      >>>>  R2313 = -1/8                                      ---Execution time = 39s---   (  or  21s  if  h < 0 )

d)  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  ( or 0 )

 STACK INPUTS OUTPUTS Y i i X j Rij L / j

Examples:    The same diagonal metric and with x = y = z = 1

1  ENTER^
1  XEQ "RIJ"  >>>>  R11 = -5/16                                      ---Execution time = 1m49s---

-Likewise

2  ENTER^
3     R/S      >>>>  R23 = -3/8                                      ---Execution time = 2m34s---

-All the components - with i <= j - are

R11 = -5/16     R22 = -13/16      R33 = -11/16
R12 = +1/8      R23 =   -3/8
R13 = +5/16

e)  Scalar Riemannian Curvature

-Contracting the Ricci tensor, we get the scalar curvature    R = gij  Rij
-"RR" calculates and stores R in register R38

-The following listing contains also "EIJ" & "WIJKL" ( cf the 2 paragraphs below )

 01 LBL "RR"  02 CLX  03 STO 38  04 RCL 09  05 STO 21  06 LBL 01  07 RCL 21  08 STO 23  09 LBL 02  10 RCL 21  11 RCL 23               12 ENTER  13 CLX  14 CIJK  15 X=0?  16 GTO 00  17 STO 31  18 RCL 21  19 RCL 23  20 X=Y?  21 GTO 02  22 RCL 31   23 ST+ 31  24 RDN  25 LBL 02  26 XROM "RIJ"  27 RCL 31  28 * 29 ST+ 38  30 LBL 00  31 DSE 23  32 GTO 02  33 DSE 21  34 GTO 01  35 RCL 38  36 RTN  37 LBL "EIJ"  38 XEQ 12  39 LBL 10  40 STO 34   41 X<>Y  42 STO 33              43 X<>Y  44 1  45 CHS  46 CIJK  47 STO 31  48 RCL 33  49 RCL 34   50 XROM "RIJ"  51 RCL 38  52 RCL 31  53 *  54 2  55 /  56 - 57 SIGN  58 RCL 33  59 RCL 38  60 RCL 34  61 X<> L  62 CLD  63 RTN  64 GTO 10  65 LBL "WIJKL"  66 XEQ 12  67 LBL 13  68 STO 34   69 RDN  70 STO 33              71 RDN  72 STO 32  73 X<>Y  74 STO 31  75 CLX  76 STO 08  77 RCL 32  78 RCL 33   79 1  80 CHS  81 CIJK  82 STO 35  83 X=0?  84 GTO 00 85 RCL 31  86 RCL 34  87 XROM "RIJ"  88 RCL 35  89 *  90 ST+ 08  91 LBL 00  92 RCL 31  93 RCL 34              94 1  95 CHS  96 CIJK  97 ST* 35  98 X=0?  99 GTO 00 100 STO 36  101 RCL 32 102 RCL 33 103 XROM "RIJ" 104 RCL 36 105 * 106 ST+ 08 107 LBL 00  108 RCL 32 109 RCL 34 110 1 111 CHS 112 CIJK 113 STO 37 114 X=0? 115 GTO 00 116 RCL 31 117 RCL 33 118 XROM "RIJ" 119 RCL 37 120 * 121 ST- 08 122 LBL 00 123 RCL 31             124 RCL 33  125 1 126 CHS 127 CIJK 128 ST* 37 129 X=0? 130 GTO 00 131 STO 36 132 RCL 32  133 RCL 34 134 XROM "RIJ" 135 RCL 36 136 * 137 ST- 08 138 LBL 00 139 RCL 37 140 RCL 35 141 - 142 RCL 38 143 * 144 RCL 09 145 1 146 - 147 / 148 RCL 08 149 + 150 RCL 09 151 2 152 - 153 / 154 STO 08 155 RCL 31 156 RCL 32 157 RCL 33 158 RCL 34             159 XROM "RIJKL" 160 ST+ 08 161 X<> 08 162 CLD 163 RTN 164 GTO 13  165 LBL 12 166 "R38=R?" 167 AVIEW 168 END

( 293 bytes / SIZE 41+n(n+1)(n+2)/2 )

 STACK INPUT OUTPUT X / R

Examples:    The same diagonal metric and with x = y = z = 1

XEQ "RR"  >>>>   R = -17/32 = R38                              ---Execution time = 5m44s---

Notes:

-With a 2D-surface, the scalar curvature is twice the Gaussian curvature given by "KGS"

-"RR" uses 5 subroutine-levels, provided the "GIJ" do not call any subroutine.
-So, this program cannot be called as more than a first level subroutine.

-Always execute "RR" before computing Einstein tensor or Weyl tensor.
-Register R38 must contain R

f) Einstein Tensor

-Einstein tensor is a symmetric tensor defined as

Eij = Rij - (1/2) R gij

-Only after executing "RR" - which stores R into R38 - execute EIJ to get the components of Einstein tensor ( see listing in paragraph above )

 STACK INPUTS OUTPUTS Y i i X j Eij L / j

Examples:    The same diagonal metric and with x = y = z = 1

1  ENTER^
1  XEQ "EIJ"  >>>>   "R38=R?"  and finally     E11 = 7/32                                      ---Execution time = 1m49s---

-Likewise

2  ENTER^
3     R/S      >>>>  E23 = -3/8                                      ---Execution time = 2m35s---

-All the components - with i <= j - are

E11 = +7/32     E22 = +1/4      E33 = +3/8
E12 = +1/8       E23 = -3/8
E13 = +5/16

Notes:

-When you XEQ "EIJ" the HP41 displays "R38=R?" to remind that R38 must contain the scalar curvature R before executing "EIJ"
-This message is not displayed again if you press R/S to get the other components.

-"EIJ" uses 5 subroutine-levels, provided the "GIJ" do not call any subroutine.
-So, this program cannot be called as more than a first level subroutine.

-In our example, the metric is diagonal. So, if i#j , Eij = Rij

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

g) 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 "RR" first before executing "WIJKL" so that data register R38 = R    ( cf listing 2 paragraphs above )

 STACK INPUTS OUTPUTS T i i Z j j Y k k X l Wijkl L / l

Examples:    The same diagonal metric and with x = y = z = 1

1  ENTER^
2  ENTER^
1  ENTER^
2  XEQ "WIJKL"  >>>>  "R38=R?"  and finally    W1212 = 0                    ---Execution time = 4m11s---

-Likewise

2  ENTER^
3  ENTER^
1  ENTER^
3     R/S      >>>>  W2313 = 0                                      ---Execution time = 3m08s---

-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

Notes:

-When you XEQ "WIJKL" the HP41 displays "R38=R?" to remind that R38 must contain the scalar curvature R before executing "WIJKL"
-This message is not displayed again if you press R/S to get another component.

-When the program stops, register R08 contains  Rijkl

-"WIJKL" uses 5 subroutine-levels, provided the "GIJ" do not call any subroutine.
-So, this program cannot be called as more than a first level subroutine.

-With the same expression for gii  if  i = 3  and  g44 = 1 + x2 y2 z2 t2    at  x = y = z = t =1   ( n = 4 )
you will get for example - after executing "INIGC" -

R = -45/32  and   W1414 = 0    W2424 = +3/16    W3434 = -3/16

h) Covariant Derivative of a Tensor Field of order 2 or 1

-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 or the tensor is covariant or contravariant:

Covariant Vector - CF 01  CF 02

Dk Vi =    k Vi - Gmik Vm

Contravariant Vector - SF 01 CF 02

Dk Vi =    k Vi + Gikm Vm

Covariant Tensor of order 2 - CF 01  CF 02

Dk Tij =   k Tij - Gmik Tmj - Gmkj Tim

Contravariant Tensor of order 2 - SF 02

Dk Tij =    k Tij + Gikm Tmj + Gjkm Tim

Mixed Tensor of order 2 - SF 01  CF 02

Dk Tij =   k Tij + Gikm Tmj - Gmkj Tim

-"DKTIJ" will calculate the components of all these tensors, but not those of order 3 , 4 , ...

Data Registers:     R00: temp
R01 thru Rnn  coordinates of the point  ( n < 8 )

R09 = n       R11 to R18 are used by "FD"              R39 = n(n+1)/2   R41 thru R40+n(n+)(n+2)/2  =  gij ,  gijGkij
R10 = h       R20 to R25 are used by "DKTIJ"

Flags:  F10   and    F01 & F02    SF 01 if there is 1 upper index      ( contravariant )
SF 02 if there are 2 upper indices ( contravariant )

Subroutines:   CIJK & "FD"

-This listing also contains other programs below

 01 LBL "DKTIJ"  02 STO 23  03 RDN  04 X<>Y  05 FRC  06 CF 10  07 X#0?  08 SF 10  09 X<> L  10 X<>Y  11 FS? 10  12 1  13 STO 22             14 RDN  15 STO 21   16 X<>Y  17 STO 19  18 CLX  19 STO 24  20 RCL 09  21 STO 25  22 FS? 10  23 GTO 10  24 LBL 09  25 RCL 23  26 RCL 22  27 RCL 25  28 FS? 02  29 X<>Y  30 CIJK  31 X=0?  32 GTO 09  33 STO 20  34 RCL 21  35 RCL 25  36 RCL 09  37 ST* Y  38 -  39 +  40 ENTER  41 SIGN  42 -  43 RCL 19  44 +  45 RCL IND X  46 XEQ IND X  47 RCL 20  48 *  49 FC? 02  50 CHS 51 ST+ 24  52 LBL 09  53 DSE 25  54 GTO 09  55 RCL 09  56 STO 25  57 LBL 10  58 RCL 23  59 RCL 21  60 RCL 25              61 FC? 01  62 FS? 02  63 X<>Y  64 CIJK  65 X=0?  66 GTO 10  67 STO 20  68 RCL 25  69 RCL 22  70 RCL 09  71 ST* Y  72 -  73 +  74 RCL 19  75 +  76 ENTER  77 SIGN  78 -  79 RCL IND X  80 XEQ IND X  81 RCL 20  82 *  83 FC? 01  84 FS? 02  85 CHS  86 ST- 24  87 LBL 10  88 DSE 25  89 GTO 10  90 RCL 21  91 RCL 22  92 RCL 09  93 ST* Y  94 -  95 +  96 RCL 19  97 +  98 1  99 - 100 RCL IND X 101 STO 00 102 RCL 10 103 RCL 23 104 XROM "FD" 105 ST+ 24 106 RCL 23 107 SIGN 108 RCL 19 109 RCL 21 110 FC? 10 111 RCL 22           112 RCL 24  113 RTN 114 LBL "RCURL" 115 XEQ 14 116 ST- 26 117 LBL 06 118 RCL 21 119 STO 22 120 ISG 22 121 LBL 07 122 RCL 26 123 RCL 22 124 + 125 RCL IND X 126 STO 00 127 RCL 10 128 RCL 21 129 XROM "FD" 130 STO 25 131 RCL 26 132 RCL 21 133 + 134 RCL IND X 135 STO 00 136 RCL 10 137 RCL 22 138 XROM "FD" 139 RCL 25 140 - 141 STO IND 23 142 CLX 143 SIGN 144 ST+ 22 145 ST+ 23 146 RCL 09 147 RCL 22 148 X<=Y? 149 GTO 07 150 SIGN 151 ST+ 21 152 CLX 153 RCL 21 154 XV^" 288 LBL 16 289 XEQ 14 290 LBL 11 291 CLX 292 STO 14 293 SIGN 294 STO 12 295 LBL 12 296 RCL 21 297 RCL 12 298 1 299 CHS 300 FC? 01 301 CLX 302 CIJK 303 RCL IND 26 304 * 305 ST+ 14 306 CLX 307 SIGN 308 ST+ 26 309 ST+ 12 310 RCL 09            311 RCL 12  312 X<=Y? 313 GTO 12 314 RCL 14 315 STO IND 23 316 RCL 09 317 ST- 26 318 SIGN 319 ST+ 21 320 ST+ 23 321 LASTX 322 RCL 21 323 X<=Y? 324 GTO 11 325 XEQ 13 326 RTN 327 GTO 16 328 LBL 14 329 STO 19 330 INT 331 STO 26 332 RCL 39 333 RCL 09 334 ST+ Z 335 2 336 + 337 * 338 41 339 + 340 XY 342 STO 23 343 STO 24 344 SIGN 345 STO 21 346 END

( 589 bytes / SIZE 41+n(n+1)(n+2)/2 )

 STACK T  INPUTS OUTPUTS V  INPUTS OUTPUTS T bbb.eee bbb.eee / / Z i i bbb.eee bbb.eee Y j j i i X k Dk Tij or Dk Tij or Dk Tij k Dk Vi   or   Dk Vi L / k / k

Examples:   With the same metric as above and still at the same point  ( x , y , z ) = ( 1 , 1 , 1 )

Covariant Tensor Field of order 2 - CF 01  CF 02

T11 = x2 + y z3      T12 = x2 y z3     T13 = x2 + y2 + z2
T21 = x2 y z           T22 = x y z        T23 = x2 + y z
T31 = x + y z         T32 = x2 y2 z2    T33 = x + y + z

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

 01  LBL "T11"  02  RCL 01   03  X^2  04  RCL 03  05  ENTER^  06  X^2  07  *  08  RCL 02  09  *  10  +  11  RTN 12  LBL "T21"  13  RCL 01   14  X^2  15  RCL 02  16  *  17  RCL 03  18  *  19  RTN  20  LBL "T31"  21  RCL 01  22  RCL 02 23  RCL 03   24  *  25  +  26  RTN  27  LBL "T12"  28  RCL 01  29  X^2  30  RCL 02  31  *  32  RCL 03  33  ST* Y 34  X^2  35  *  36  RTN  37  LBL "T22"  38  RCL 01   39  RCL 02  40  RCL 03  41  *  42  *  43  RTN  44  LBL "T32" 45  RCL 01  46  RCL 02  47  RCL 03   48  *  49  *  50  X^2  51  RTN  52  LBL "T13"  53  RCL 01  54  X^2  55  RCL 02 56  X^2  57  RCL 03   58  X^2  59  +  60  +  61  RTN  62  LBL "T23"  63  RCL 01  64  X^2  65  RCL 02  66  RCL 03 67  *  68  +  69  RTN  70  LBL "T33"  71  RCL 01  72  RCL 02   73  RCL 03  74  +  75  +  76  RTN  77  END

-Store the names of theses subroutines in 9 consecutive registers, without disturbing R41 to R70 , or R39 , R09 , R10 , R21 to R25 which are used by "DKTIJ"

-For instance, in R27 thru R35 ->  control number = 27.035

T11      T12     T13                      R27     R30     R33
T21      T22     T23    stored in     R28     R31     R34
T31      T32     T33                      R29     R32     R35

-If you want to get D3 T12

Covariant Tensor Field of order 2 - CF 01  CF 02

CF 01   CF 02       ( and still  SF 00  since the metric is orthogonal )

27.035    ENTER^
1        ENTER^
2        ENTER^
3        XEQ "DKTIJ"    >>>>    D3 T12  =  1/4                                             ---Execution time = 16s---

Contravariant Tensor Field of order 2 - CF 01  SF 02

-Assuming we keep the same expressions but the tensor is Tij instead of Tij

CF 01    SF 02

27.035    ENTER^
1        ENTER^
2        ENTER^
3        XEQ "DKTIJ"    >>>>    D3 T12  =  5/4                                             ---Execution time = 16s---

Mixed Tensor Field of order 2 - SF 01  CF 02

-Assuming again the same expressions with the tensor Tij

SF 01    CF 02

27.035    ENTER^
1        ENTER^
2        ENTER^
3        XEQ "DKTIJ"    >>>>    D3 T12  =  3/4                                             ---Execution time = 16s---

Covariant Vector Field - CF 01  CF 02

-Let's keep the first 3 functions as the components of the covariant vector field

CF 01    CF 02

27.029    ENTER^
2        ENTER^
3        XEQ "DKTIJ"    >>>>    D3 V2  =  -1/4                                             ---Execution time = 12s---

Contravariant Vector Field - SF 01  CF 02

-Let's take again the first 3 functions as the components of the contravariant vector field

SF 01    CF 02

27.029    ENTER^
2        ENTER^
3        XEQ "DKTIJ"    >>>>    D3 V2  =  +1/4                                             ---Execution time = 12s---

i)  Curl of a Covariant Vector Field, Divergence of a Contravariant
Vector Field, Laplacian & Gradient of a Scalar 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

-So, to save space, "RCURL" only calculates and stores in consecutive registers  Curlij V with i < j  &  i = 1 , 2 , ... , n-1   i-e  n(n-1)/2 components

Data Registers:          R00: temp
•  R01 thru Rnn  coordinates of the point  ( n < 8 )

•  R09 = n      R20 to R26 are used by "RCURL"       •  R39 = n(n+1)/2   R41 thru R40+n(n+)(n+2)/2  =  gij ,  gijGkij
•  R10 = h

Flags: /
Subroutine:         "FD"

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

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

X = x2 + y z3
Y = x2 y z2
Z = x + y z

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

-If you store the function names in R31 R32 R33 -> control number = 31.033

31.033   XEQ "RCURL"  >>>>   71.073                                  ---Execution time = 53s---

-And you get

R71 = -1                                                      0    +1   -2
R72 = 2       so, the complete result is          -1     0    -1
R73 = 1                                                        2     1      0

Notes

-The HP41 stores the results in a block of consecutive registers starting at RBB
where  BB = Max ( bb+n , 41+n(n+1)(n+2)/2 )
-So, the data are not disturb.

-This definition will not give the same results as "CDGL" 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 Rnn  coordinates of the point  ( n < 8 )

•   R09 = n      R19 to R25 are used by "RDIV"       •  R39 = n(n+1)/2   R41 thru R40+n(n+)(n+2)/2  =  gij ,  gijGkij
•  R10 = h

Flags: /
Subroutines:    CIJK     "FD"

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

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

X = x2 + y z3
Y = x2 y z2
Z = x + y z

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

-If you store the function names in R31 R32 R33 -> control number = 31.033

31.033   XEQ "RDIV"  >>>>    41/4                                 ---Execution time = 33s---

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

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

•  R09 = n      R20 to R26 are used by "RLAPGR"       •  R39 = n(n+1)/2   R41 thru R40+n(n+)(n+2)/2  =  gij ,  gijGkij
•  R10 = h      R11 thru R19 are used by "SD"

Flags: /
Subroutines:    CIJK     "FD"  "SD"  and a function that computes f assuming x1 is in R01 ..................  xn is in Rnn

 STACK INPUTS OUTPUTS Y / bbb.eeeGRAD X / Lapl (f)

With bbb = 41+n(n+1)(n+2)/2

Example:   Again the same metric at the same point ( 1 , 1 , 1 ) and the following scalar field:

f(x,y,z) = x2 + y2 z3

 01  LBL "FF"  02  RCL 01  03  X^2  04  RCL 02  05  RCL 03  06  *  07  X^2  08  RCL 03  09  *  10  +  11  RTN

FF  ASTO 00
XEQ "RLAPGR"  >>>>     61/16                                                                      ---Execution time = 56s---
X<>Y    71.073

And you get the covariant components of the gradient in R71-R72-R73 i-e  ( 2 , 2 , 3 )

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

j) Covariant Vector <> Contravariant Vector

-If flag F01 is clear  "V<>V^"  transforms the covariant components of a vector into its contravariant components
-If flag F01 is set, the reverse operation is performed.

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^"       •  R39 = n(n+1)/2   R41 thru R40+n(n+)(n+2)/2  =  gij ,  gijGkij
•  R10 = h

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

 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 , 3 , 4 )
-Let's calculate the contravariant components.

-If you store 2  3  4  into  R31  R32  R33  respectively

CF 01
31.033   XEQ "V<>V^"  >>>>  71.073                                                                       ---Execution time = 8s---

-And you get the contravariant components  ( 1 , 3/4 , 1 )  in R71-R72-R73

-To recalculate the covariant components

SF 01   71.073   XEQ "V<>V^"   >>>>   74.076    and the covariant components ( 2 , 3 , 4 )  are in R74-R75-R76

Note:

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

8°)  An Example of a non-diagonal metric Tensor ( 3D-manifolfd )

g11 = 1 + x2 y2 z2                            g12 = x y
g22 = 1 + x2 + y2 + z2                     g13 = x z
g33 = 1 + x2 y2 + x2 z2 + y2 z2         g23 = y z

-Load for instance the 6 routines in main memory:

 01  LBL "G11"  02  CLX  03  SIGN  04  RCL 01   05  RCL 02  06  RCL 03   07  *  08  *  09  X^2  10  +  11  RTN  12  LBL "G22" 13  CLX   14  SIGN  15  RCL 01   16  X^2  17  +  18  RCL 02       19  X^2  20  +  21  RCL 03  22  X^2  23  +  24  RTN 25  LBL "G33"  26  CLX  27  SIGN  28  RCL 01      29  RCL 02   30  *  31  X^2  32  +  33  RCL 01  34  RCL 03  35  *  36  X^2 37  +  38  RCL 02  39  RCL 03   40  *  41  X^2  42  +  43  RTN  44  LBL"G12"  45  RCL 01  46  RCL 02  47  *  48  RTN 49  LNL "G13"  50  RCL 01  51  RCL 03   52  *  53  RTN  54  LBL "G23"  55  RCL 02  56  RCL 03  57  *  58  RTN  59  END

-At the point  x = y = z = 1

Christoffel Symbols:

1  STO 01  STO 02  STO 03    ( SIZE 071 at least )   and if you choose  h = 0.1

•  SF 03   to avoid the test with alpha data - unuseful here since all the gijare different from 0
•  CF 04   to get the Gkij

0.1  ENTER^
3    XEQ "INIGC"   >>>>   41.070                                                               ---Execution time = 21m24s---

-And you find in registers R41 thru R52

R41 = g11 = 2     R42 = g12 = 1      R44 = g13 = 1              R47 = g11 = +5/8      R48 =  g12 = -1/8           R50 = g13 = -1/8
R43 = g22 = 4     R45 = g23 = 1                                                R49 = g22 = +7/24         R51 =  g23 = -1/24
R46 = g33 = 4                                                                                      R52 = g33 = +7/24

>>>    R40 = g = det gij = 24

-And in registers R53 to R70

R53 = G111 = 5/8    R54 = G112 =  1/2     R56 = G113 =  3/8
R55 = G122 = -1/8     R57 = G123 = -3/8
R58 = G133 = -3/4

R59 = G211 = -1/8     R60 = G212 = 1/6     R62 = G213 = -5/24
R61 = G222 = 7/24    R63 = G223 = 5/24
R64 = G233 = -1/4

R65 = G311 = -1/8          R66 = G312 = -1/6           R68 = G313 =  11/24
R67 = G322 = -1/24          R69 = G323 =  13/24
R70 = G333 =   3/4

-We have also   R39 = 6  i-e  n(n+1)/2

Curvature Tensor:

2  ENTER^
3  ENTER^
1  ENTER^
3  XEQ "RIJKL"  >>>>  R2313 = -7/24                                       ---Execution time = 86s---

2  ENTER^
1  ENTER^
2  ENTER^
1  XEQ "RIJK^L"  >>>>  R1212 = -13/96                                     ---Execution time = 232s---

Ricci Tensor:

1  ENTER^
3  XEQ "RIJ"      >>>>  R13 = -17/96                                         ---Execution time = 12m26s---

-All the components - with i <= j - are

R11 = +7/96       R22 = -79/288      R33 = -5/144
R12 = -28/96      R23 = -277/288
R13 = -17/96

Scalar Curvature:

XEQ "RR"  >>>>   R = +11/72 = R38                                                   ---Execution time = 59mn36s---        ( 40m55s  with  h < 0 )

>>> You can change the h-value in R10 to get faster results, for example,  h = - 0.01  STO 10 )

Einstein Tensor: ( after executing "RR" )

2  ENTER^
3  XEQ "EIJ"      >>>>  "R38=R?"  and   E23 = -299/288                               ---Execution time = 6m59s---   ( with  h < 0 )

-All the components - with i <= j - are

E11 = -23/288         E22 = -167/288      E33 = -49/144
E12 = -53/144         E23 = -299/288
E13 = -73/288

Weyl Tensor:   ( after executing "RR" )

1  ENTER^
2  ENTER^
1  ENTER^
2  XEQ "WIJKL"  >>>>  "R38=R?"  and finally    W1212 = 0

-Not a surprise since W = 0 if n < 4.

Differential Operators:

-With the same fields as above ( FF  X  Y  Z  T11 ..... T33 )  you will get:

Lap(f)   = 317/96                                      ---Execution time = 140s---

Curlij & Gradk are unchanged because they don't depend on the metric.

Div(X Y Z ) = 10.5                                 ---Execution time = 33s---

D3T12 = 31/24                                          ---Execution time = 18s---

-Likewise,      D3T33 = 1         D1T23 = -5/12          D2T22 = 19/24     ... and so on ...

Notes:

-As you can see, though dim = 3 , several routines are very slow, especially "RR" to get the scalar curvature.
-That's why a good emulator is recommended.

-Practically, there are often several gij = 0, even if the metric is not diagonal and execution time could be much smaller.
-Remember also that the programs run faster if you choose h < 0.

9°) Schwarzchild Metric ( 4D-Spacetime )

-Expressed in spherical coordinates ( r , u , v , ct ), we have:

ds2 = ( 1 - a / r ) -1 dr2 + r2 ( du2 + Sin2 u dv2 ) - ( 1 - a / r )  c2 dt2        where    a = 2GM/c2

-Let's store a in register R07 and load the following routines ( we must add the 6 functions LBL "G12" ........... LBL "G34" )

 01  LBL "G11"  02  CLX  03  SIGN  04  RCL 07  05  RCL 01  06  /  07  -  08  1/X 09  RTN  10  LBL "G22"  11  RCL 01  12  X^2  13  RTN  14  LBL "G33"  15  RCL 01  16  RCL 02 17  SIN  18  *  19  X^2  20  RTN  21  LBL "G44"  22  RCL 07  23  RCL 01  24  / 25  ENTER^  26  CLX  27  SIGN  28  -  29  RTN  30  LBL "G12"  31  LBL "G13"  32  LBL "G14" 33  LBL "G23"  34  LBL "G24"  35  LBL "G34"  36  ASTO X  37  RTN  38  END

-Let's take, for example   r = 2    u = PI/6   v = 0   ct = 0    and   a = 1

2  STO 01    PI  6  /  STO 02   CLX   STO 03   STO 04    1  STO 07

Initialization:

XEQ "RAD"     CF 04   CF 03

-Dimension = n = 4  and if you choose  h = 0.1

0.1   ENTER^
4    XEQ INIGC    >>> >   41.100                                                        ---Execution time = 8m18s---

And the Christoffel symbols are  ( after rounding the results given by the HP41 )

R61 = G111 = -1/4    R62 = G112 =  0        R64 = G113 =    0      R67 = G144 =   0
R63 = G122 = -1       R65 = G123 =   0      R68 = G134 =   0
R66 = G133 = -1/4     R69 = G124 =   0
R70 = G114 = 1/16

R71 = G211 =   0     R72 = G212 = 1/2     R74 = G213 =   0              R77 = G244 =  0
R73 = G222 =   0      R75 = G223 =   0              R78 = G234 =   0
R76 = G233 = -sqrt(3)/4    R79 = G224 =    0
R80 = G214 =   0

R81 = G311 =   0     R82 = G312 =  0        R84 = G313 =  1/2      R87 = G314 =  0
R83 = G322 =  0        R85 = G323 = sqrt(3)  R88 = G324 =  0
R86 = G333 =   0        R89 = G334 =   0
R90 = G344 =   0

R91 = G411 =   0    R92 = G412 =   0    R94 = G413 =  0    R97 = G414 = 1/4
R93 = G422 =   0    R95 = G423 =   0    R98 = G424 =  0
R96 = G433 =   0    R99 = G434 =   0
R100 = G444 = 0

Curvature Tensor:

1  ENTER^
2  ENTER^
1  ENTER^
2  XEQ "RIJKL"  >>>>  R1212 = -1/2       ( HP41 gives  -0.500001321 )                        ---Execution time = 30s---

1  ENTER^
4  ENTER^
1  ENTER^
4     R/S     >>>>  R1414 = -0.124999921      ( exact result = -1/8 )

3  ENTER^
4  ENTER^
3  ENTER^
4     R/S     >>>>  R3434 =  1/32                    ... and so on ...

Ricci Tensor:

2  ENTER^
3  XEQ "RIJ"      >>>>  R23 = 0                                                                ---Execution time = 3m14s---

-All the components are actually zero.
-However, roundoff-errors are unavoidable and the HP41 sometime returns small numbers instead of 0, for example:

R11 = -0.000000820

Scalar Curvature:

XEQ "RR"  >>>>   R =  -0.000000846    ( exact result = 0 )

Einstein Tensor:

2  ENTER^
3  XEQ "EIJ"      >>>>  E23 = 0

-Here again, all the components are 0

Weyl Tensor:

1  ENTER^
2  ENTER^
1  ENTER^
2  XEQ "WIJKL"  >>>>  "R38=R?"  and finally    W1212 = -0.500000095   ( exact result = -1/2 )

2  ENTER^
3  ENTER^
2  ENTER^
3     R/S     >>>>      W2323 = 0.500000097   ( exact result = +1/2 )

1  ENTER^
4  ENTER^
1  ENTER^
4     R/S     >>>>      W1414 = -0.125000024   ( exact result = -1/8 )

3  ENTER^
4  ENTER^
3  ENTER^
4     R/S     >>>>      W3434 = 0.031250006   ( exact result = +1/32 )     ....  and so on ...

10°) Non-Riemannian Geometry

-In Riemannian manifolds, the metric tensor determines the connexion by the Christoffel symbols.
-In more general differential manifolds, the linear connexions Gkij are independant functions and may be non-symmetric.

>>>  So, you have to load in main memory  n3  functions called   LBL "C111"  LBL "C112"  LBL "C121" and so on...

-Registers R09 & R10 are also to be initialized before executing the programs below:

•  R09 = n < 9
•  R10 = h

-And of course,

•  R01 = x1  .................  •  Rnn = xn    ( n < 9 )

>>>  If the tensor depends on the metric tensor, execute "INIGC" first - with SF 04 to avoid calculating the Christoffel symbols.

a) Curvature Tensor (1-3)

-If  CF 02,  "2RIJK^L" calculates the curvature tensor defined as

Rlijk = jGlik - kGlij + GlmjGmik - GlmkGmij

-If  SF 02   this routine calculates the curvature tensor defined by the transposed connexion  Gkji

Rlijk = jGlki - kGlji + GljmGmki - GlkmGmji

-This listing also contains other programs below.

 01 LBL "2RIJK^L"  02 LBL 10  03 STO 14  04 RDN  05 STO 13  06 RDN  07 STO 12  08 X<>Y  09 STO 11  10 CLX  11 STO 20  12 RCL 09  13 STO 15                    14 LBL 01  15 RCL 15  16 RCL 13   17 FS? 02  18 X<>Y  19 RCL 14   20 -DIFF MANIF   21 ASTO X  22 XEQ IND X  23 SIGN  24 X#0?  25 LASTX  26 STO 26  27 RCL 11  28 RCL 12  29 FS? 02  30 X<>Y  31 RCL 15  32 -DIFF MANIF  33 ASTO X  34 XEQ IND X  35 SIGN  36 X#0? 37 LASTX  38 ST* 26  39 RCL 15  40 RCL 12  41 FS? 02  42 X<>Y  43 RCL 14  44 -DIFF MANIF  45 ASTO X  46 XEQ IND X  47 SIGN  48 X#0?  49 LASTX  50 STO 19   51 RCL 11                   52 RCL 13  53 FS? 02  54 X<>Y  55 RCL 15  56 -DIFF MANIF   57 ASTO X  58 XEQ IND X  59 SIGN  60 X#0?  61 LASTX  62 RCL 19  63 *  64 RCL 26  65 -  66 ST+ 20  67 DSE 15  68 GTO 01  69 RCL 11  70 STO 21  71 RCL 12  72 STO 22 73 FS? 02  74 X<>Y  75 RCL 14  76 STO 24  77 -DIFF MANIF  78 ASTO 00  79 RCL 10  80 RCL 13  81 STO 23  82 XROM "FD"  83 ST- 20  84 RCL 21  85 RCL 23                   86 FS? 02  87 X<>Y  88 RCL 24  89 -DIFF MANIF   90 ASTO 00  91 RCL 10  92 RCL 22  93 XROM "FD"  94 ST+ 20  95 RCL 24  96 SIGN  97 RCL 21  98 RCL 22  99 RCL 23 100 RCL 20 101 RTN 102 GTO 10 103 LBL "RQIJ" 104 LBL 12 105 STO 23 106 X<>Y 107 STO 28 108 RCL 09 109 STO 24 110 CLX 111 STO 27 112 LBL 02 113 RCL 28 114 RCL 24 115 FS? 00 116 X<>Y 117 RCL 23 118 RCL 24 119 XROM "2RIJK^L" 120 ST+ 27 121 DSE 24 122 GTO 02 123 RCL 23                 124 SIGN 125 RCL 28 126 RCL 27 127 RTN 128 GTO 12 129 LBL "RQ" 130 LBL 09 131 RCL 09 132 STO 28 133 CLX 134 STO 38 135 LBL 03 136 RCL 09 137 STO 29 138 LBL 04 139 RCL 28 140 RCL 29 141 ENTER 142 CLX 143 CIJK 144 X=0? 145 GTO 00 146 STO 30 147 RCL 28 148 RCL 29 149 XROM "RQIJ" 150 RCL 30 151 * 152 ST+ 38 153 LBL 00 154 DSE 29 155 GTO 04 156 DSE 28 157 GTO 03 158 RCL 38                 159 RTN 160 GTO 09 161 LBL "SIJK" 162 LBL 05 163 STO 13 164 RDN 165 STO 12 166 X<>Y 167 STO 11 168 R^ 169 -DIFF MANIF 170 ASTO X 171 XEQ IND X 172 SIGN 173 X#0? 174 LASTX 175 STO 14 176 RCL 11 177 RCL 12 178 RCL 13 179 -DIFF MANIF 180 ASTO X 181 XEQ IND X 182 SIGN 183 X#0? 184 LASTX 185 RCL 14 186 - 187 STO 14 188 RCL 13 189 SIGN 190 RCL 11 191 RCL 12                  192 RCL 14  193 RTN 194 GTO 05 195 LBL "SJ" 196 LBL 06 197 STO 11 198 RCL 09  199 STO 12 200 CLX 201 STO 15 202 LBL 07 203 RCL 11 204 RCL 12 205 ENTER 206 XROM "SIJ" 207 ST+ 15 208 DSE 12 209 GTO 07 210 RCL 11 211 SIGN 212 RCL 15 213 RTN 214 GTO 06 215 END

( 349 bytes /  SIZE var. )

 STACK INPUTS OUTPUTS T i i Z j j Y k k X l Rlijk L / l

Example:   Let's take the 3D-manifold with the same metric as above at the point  x = y = z = 1

g11 = 1 + x2 y2 z2                            g12 = 0
g22 = 1 + x2 + y2 + z2                     g13 = 0
g33 = 1 + x2 y2 + x2 z2 + y2 z2         g23 = 0

-Load the 6 routines in main memory:

 01  LBL "G11"  02  CLX  03  SIGN  04  RCL 01  05  RCL 02  06  RCL 03  07  *  08  *  09  X^2  10  + 11  RTN  12  LBL "G22"  13  CLX   14  SIGN  15  RCL 01  16  X^2  17  +  18  RCL 02  19  X^2  20  + 21  RCL 03  22  X^2  23  +  24  RTN  25  LBL "G33"  26  CLX  27  SIGN  28  RCL 01  29  RCL 02  30  * 31  X^2  32  +  33  RCL 01   34  RCL 03       35  *  36  X^2  37  +  38  RCL 02  39  RCL 03  40  * 41  X^2  42  +  43  RTN  44  LBL "G12"  45  LBL "G13"  46  LBL "G23"  47  ASTO X  48  RTN  49  END

-But let's change the connexion as follows:

G111 = x y2 z2  /  g11           G122 = - x / g11
G112 = x2 y z2  /  g11  = G121G123 = G132 = 0
G113 = x2 y2 z  /  g11  = G131G133 = - x ( y2 + z2 ) / g11

G211 = - x2 y z2   /  g22       G222 =  y / g22
G212 = x  /  g22  =  G221      G223 = G232 = z / g22
G213 = 0  =  G231               G233 = - y ( x2 + z2 ) / g22

G311 = - x2 y2 z  /  g33                G322 =  - z / g33                                      so, it is the same as the Riemannian connexion
G312 =  0  =  G321                      G323 = - G332 = y ( x2 + z2 ) / g33            except that  G323 = - G332
G313 = x ( y2 + z2 ) / g33  =  G331G333 = z ( x2 + y2 ) / g33

-We have to load these 27 functions in main memory:

 01 LBL "C111"  02 XEQ "G11"  03 RCL 02  04 RCL 03  05 *  06 X^2  07 RCL 01  08 *  09 X<>Y  10 /  11 RTN  12 LBL "C121"  13 LBL "C211"  14 XEQ "G11"  15 RCL 01  16 RCL 03  17 *  18 X^2  19 RCL 02  20 *  21 X<>Y  22 /  23 RTN  24 LBL "C131"  25 LBL "C311"  26 XEQ "G11"  27 RCL 01  28 RCL 02  29 *  30 X^2 31 RCL 03  32 *  33 X<>Y  34 /  35 RTN  36 LBL "C221"  37 XEQ "G11"  38 RCL 01  39 CHS  40 X<>Y  41 /  42 RTN  43 LBL "C231"  44 LBL "C321"  45 CLX  46 RTN  47 LBL "C331"  48 XEQ "G11"  49 RCL 02  50 X^2  51 RCL 03  52 X^2  53 +  54 RCL 01  55 *  56 CHS  57 X<>Y  58 /  59 RTN  60 LBL "C112" 61 XEQ "G22"  62 RCL 01  63 RCL 03  64 *  65 X^2  66 RCL 02  67 *  68 CHS  69 X<>Y  70 /  71 RTN  72 LBL "C122"  73 LBL "C212"  74 XEQ "G22"  75 RCL 01  76 X<>Y  77 /  78 RTN  79 LBL "C132"  80 LBL "C312"  81 CLX  82 RTN  83 LBL "C222"  84 XEQ "G22"  85 RCL 02  86 X<>Y  87 /  88 RTN  89 LBL "C232"  90 LBL "C322" 91 XEQ "G22"  92 RCL 03  93 X<>Y  94 /  95 RTN  96 LBL "C332"  97 XEQ "G22"  98 RCL 01  99 X^2 100 RCL 03 101 X^2 102 + 103 RCL 02 104 * 105 CHS 106 X<>Y 107 / 108 RTN 109 LBL "C113" 110 XEQ "G33" 111 RCL 01 112 RCL 02 113 * 114 X^2 115 RCL 03 116 * 117 CHS 118 X<>Y 119 / 120 RTN 121 LBL "C123" 122 LBL "C213" 123 CLX 124 RTN 125 LBL "C133" 126 LBL "C313" 127 XEQ "G33" 128 RCL 02 129 X^2 130 RCL 03 131 X^2 132 + 133 RCL 01 134 * 135 X<>Y 136 / 137 RTN 138 LBL "C223" 139 XEQ "G33" 140 RCL 03 141 CHS 142 X<>Y 143 / 144 RTN 145 LBL "C323" 146 SF 09 147 GTO 00 148 LBL "C233" 149 CF 09 150 LBL 00 151 XEQ "G33" 152 RCL 01 153 X^2 154 RCL 03 155 X^2 156 + 157 RCL 02 158 * 159 X<>Y 160 / 161 FC? 09 162 RTN 163 CHS 164 RTN 165 LBL "C333" 166 XEQ "G33" 167 RCL 01 168 X^2 169 RCL 02 170 X^2 171 + 172 RCL 03 173 * 174 X<>Y 175 / 176 RTN 177 END

-Though this is not necessary to calculate Rlijk  let's execute "INIGC"  with SF04  to calculate the metric tensor ( covariant & contravariant ) only

1  STO 01  STO 02  STO 03   CF 03  SF 04

0.1  ENTER^
3   XEQ "INIGC"   >>>>   41.052                                                   ---Execution time = 33s---

-We get the same values as before:

R41 = g11 = 2     R42 = g12 = 0      R44 = g13 = 0              R47 = g11 = 1/2     R48 =  g12 = 0            R50 = g13 =  0
R43 = g22 = 4     R45 = g23 = 0                                             R49 = g22 = 1/4          R51 =  g23 = 0
R46 = g33 = 4                                                                               R52 = g33 = 1/4

-Then,  with  CF 02

3  ENTER^
2  ENTER^
3  ENTER^
XEQ "2RIJK^L"  >>>>  R3323 = 0                                      ---Execution time = 62s---

-Likewise with SF 02  ( transposed connexion )

3  ENTER^
2  ENTER^
3  ENTER^
R/S      >>>>  R3323 = +1/4                                      ---Execution time = 59s---

b) Ricci Tensors & Segmental Curvatures

-Ricci tensor is  Rij = Rmimj

-There is also a segmental curvature tensor  Qij = Rmmij   which equals zero in Riemannian manifolds.

Qij = iGmmj - jGmmi  is actually a curl.

If  CF 00  "RQIJ"  returns the components of  Rij
-If  CF 02 with the connexion Gkij
-If  SF 02 with the connexion Gkji

If  SF 00  "RQIJ" returns the components of Qij
-If  CF 02 with the connexion Gkij
-If  SF 02 with the connexion Gkji

 STACK INPUTS OUTPUTS Y i i X j Rij or Qij L / j

Rij  if CF 00   ( connexion Gkij if CF 02 ,  transpose connexion Gkji  if SF 02 )
Qij  if SF 00   ( connexion Gkij if CF 02 ,  transpose connexion Gkji  if SF 02 )

Examples:

CF 00   CF 02

2  ENTER^
3  XEQ "RQIJ"   >>>>   R23 = - 3/8                                      ---Execution time = 2m44s---

-Remark that  R32 = - 7/8   Ricci tensor is no more symmetric

CF 00   SF 02

2  ENTER^
3     R/S       >>>>   R23 = - 9/8

SF 00   CF 02

2  ENTER^
3     R/S       >>>>   Q23 =  0

SF 00   SF 02

2  ENTER^
3      R/S      >>>>   Q23 =  0

c)  Scalar Curvatures

-The contracted tensors  R = gij Rij  and  Q = gij Qij  gives 2 scalars.

-To obtain correct results, "INIGC" must be executed before "RQ" to store  gij  in the proper registers.

 STACK INPUT OUTPUT X / R or Q

R if CF 00   ( connexion Gkij if CF 02 ,  transpose connexion Gkji  if SF 02 )
Q if SF 00   ( connexion Gkij if CF 02 ,  transpose connexion Gkji  if SF 02 )

Examples:

CF 00   CF 02

2  ENTER^
3  XEQ "RQ"   >>>>   R = - 17/32                               ---Execution time = 9m40s---

CF 00   SF 02

2  ENTER^
3     R/S     >>>>   R = - 9/16

SF 00   CF 02

2  ENTER^
3     R/S     >>>>   Q =  0

SF 00   SF 02

2  ENTER^
3     R/S     >>>>   Q =  0

Notes:

-A more important modification of the connection would be necessary to get Q # 0
-"RQ" uses 5 subroutine levels, so the functions "CIJK" must use less than 2 subroutine levels...

d) Torsion Tensor

-Since the Christoffel symbols are symmetric in Riemannian manifolds, GkijGkji = 0
-In the general case, we get a tensor of order 3 that is called the torsion tensor

Skij = GkijGkji                ( some authors define  Skij = (1/2) ( GkijGkji )    others by   Skij = ( GkjiGkij )  )

-Though Gkij is not a tensor, Skij is a tensor !

 STACK INPUTS OUTPUTS Z i i Y j j X k Skij L / k

Examples:

2  ENTER^
3  ENTER^
XEQ "SIJK"  >>>>   S323 = +1

3  ENTER^
2  ENTER^
3     R/S      >>>>  S332 = -1

e) Torsion Vector

-Contracting  Skij  we get a covariant vector  Sj = Skjk

 STACK INPUT OUTPUT X j Sj L / j

Examples:

1  XEQ "SJ"  >>>>  S1 =  0
2      R/S       >>>>  S2 = +1
3      R/S       >>>>  S3 =  0

Note:

-I called this program "SJ" instead of "SI" to avoid a conflict with the Sine Integral...

11°)  Differential Geometry QRG

 Functions DESCRIPTION REMARKS KHT Curvature & Torsion 3D-parametric Curve            X(t) Y(t) Z(t)               h ­ t   ->   Khi(t)  X<>Y  Tau(t) / KGS Gaussian Curvature of a surface z = f(x,y)           •  R00 = f name        h ­ y ­ x   ->  K(Gauss) / dF 1st- 2nd- 3rd-derivative of a function of 1 var.    •  R00 = f name               h ­ x   ->  f'(x)  RDN  f"(x)  RDN  f'''(x) / dF2 Derivatives of a function of 2 variables    •  R00 = f name     h ­ y ­ x   ->   f 'x  RDN  f 'y  RDN  f ''xx  RDN  f ''yy / dF3 Derivatives of a function of 3 variables   •  R00 = f name   h ­ z ­ y ­ x  ->  f 'x   RDN  f 'y    RDN   f 'z   RDN  Df / MDV Mixed Derivative of a function of 2 var.                •  R00 = f name                 h ­ y ­ x    ->  f ''xy / BHRM Biharmonic Operator of a function of 3 var.          •  R00 = f name           h ­ z ­ y ­ x    ->  D2f / CDGL Curl-Diverg-Grad-Lapl                           h ­ z ­ y ­ x  ->  Curl1  RDN Curl2 RDN Curl3 RDN  Div E     E = [ X(x,y,z) Y(x,y,z) Z(x,y,z) ]                                           L = 4.019 ( contr. numb. all results ) CF01 CF 02 = Rect Coord  SF 01 CF 02 = Cyl Coord  CF 01 SF 02 = Sph. Coord. FD First Deriv. funct n var. ( n < 10 )   •  R00 = f name  •  R01 = x1 , ..... , •  Rnn = xn        h ­ i  ->   ¶f / ¶xi SF03 -> X'=0 if X = alpha             h > 0 = order10            h < 0 = order 4 SD 2nd Deriv.  funct n var. ( n < 10 )   •  R00 = f name  •  R01 = x1 , ..... , •  Rnn = xn   h ­ j ­ i  ->  ¶2f / ¶xi¶xj X'=0 if X = alpha             h > 0 = order10            h < 0 = order 4 LAPL Laplacian   funct n var. ( n < 10 )   •  R00 = f name  •  R01 = x1 , ..... , •  Rnn = xn         h ­ n  ->   D f X'=0 if X = alpha            h > 0 = order10            h < 0 = order 4 LS Solves linear system,  p = small number for tiny elements, rr = nb rows               bbb.eee,rr ­ p  ->  det M / / N-DIMENSIONAL RIEMANNIAN  MANIFOLDS / INIGC Initialization:  •  R01 = x1 , ..... , •  Rnn = xn        h ­ n  ->  41.eee = control number of gij , gij  and Gkij SF 04 -> gij , gij  only CIJK Recalls  Gkij  assuming INIGC has been executed      i ­ j ­ k  -> Gkij / RIJKL Covariant Curvature Tensor      i ­ j ­ k ­ l  ->  Rijkl = (1/2) ( ¶jk gil + ¶il gjk - ¶jl gik - ¶ik gjl )                                                                                           + gmp ( GmjkGpil - GmjlGpik ) / RIJK^L Curvature Tensor ( 1-3 )           i ­ j ­ k ­ l  ->  Rlijk = glm Rmijk / RIJ Ricci Tensor                                         i ­ j   ->  Rij  = gkm  Rikjm / RR Scalar Riemannian Curvature            no input  ->  R =  gij  Rij / EIJ Einstein Tensor                                     i ­ j   ->  Eij = Rij - (1/2) R gij XEQ "RR" first to get R38 = R WIJKL Weyl Tensor                            i ­ j ­ k ­ l  ->  Wijkl = Rijkl + [ 1/(n-2) ] ( Ril gjk - Rik gjl + Rjk gil - Rjl gik )                                                                                               + [ R/(n-1)/(n-2) ] ( gik gjl - gil gjk ) XEQ "RR" first to get R38 = R DKTIJ Covariant Derivative of a Tensor    bbb.eee ­  i ­ j ­ k ->  Dk Tij or Dk Tij or Dk Tij              or a Vector                            bbb.eee ­  i ­ k  ->       Dk Vi   or   Dk Vi CF 01  CF 02  covariant          SF 01      1 contrav          SF 02      2 contrav RCURL Curl of a Covariant Vector Field                bbb.eee  ->  Curlij V = ¶i Vj - ¶j Vi / RDIV Divergence of a Contrav. Vector Fied        bbb.eee  ->  Div V =  Di Vi =     ¶i Vi + Giik Vk / RLAPGR Laplacian & Gradient of a Scalar Field   •  R00 = f name   no input ->  Lapl(f) = gjk ( ¶jk f  -  Gijk ¶i f  )                                                                                                    X<>Y  bbb.eee(Gradient) / V<>V^ Exchanging Covariant & Contravariant Coord. of a Vector             bbb.eee(Vi) ->  bbb.eee(Vi) SF 01 = inverse operation / GENERAL DIFFERENTIAL MANIFOLDS including NON-RIEMANNIAN MANIFOLDS / 2RIJK^L 2 Curvature Tensors: with the connexion  Gkij  & with the transpose connexion  Gkji            i ­ j ­ k ­ l  ->  Rlijk   Rlijk = ¶jGlik - ¶kGlij + GlmjGmik - GlmkGmij   OR   Rlijk = ¶jGlki - ¶kGlji + GljmGmki - GlkmGmji CF 02 = first Tensor   SF 02 = second Tensor RQIJ Ricci Tensor(s)                           i ­ j   ->  Rij = Rmimj   & Segmental Curvatue Tensor(s)     i ­ j   ->  Qij = Rmmij   = ¶iGmmj - ¶jGmmi CF 00 = Ricci Tensor   SF 00 = Segmental Tensor   CF 02 & SF 02 as above RQ Scalar Curvatures      no input ->  R = gij Rij                                   no input -> Q = gij Qij CF 00 = R              SF 00 = Q   CF 02 & SF 02 as above SIJK Torsion Tensor          i ­ j ­ k ->   Skij = Gkij -  Gkji / SJ Torsion Vector                     j  ->   Sj = Skjk /

References:

[1]  Elie Cartan - "Leçons sur la géométrie des espaces de Riemann" - Gauthier-Villars, Paris   ( in French )
[2]  Denis-Papin & Kaufmann - "Cours de calcul Tensoriel Appliqué" - Albin Michel    ( in French )
[3]  List of Formulas in Riemannian Geometry
[4]  Nikodem J. Poplawski - "Spacetime and fields"- Department of Physics, Indiana University, Bloomington, IN 47405, USA
[5]  Marie-Antoinette Tonnelat - Theorie unitaire affine du champ physique. J. Phys. Radium, 1951, 12 (2),
pp.81-88. <10.1051/jphysrad:0195100120208100>. <jpa-00234360>