DIFF GEOM

# Differential Geometry for the HP41

Overview

0°)  An M-Code Routine
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^ + Another M-Code Routine
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
k) Dot Product

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

-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, especially if the metric is no diagonal.

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

>>> Last update:

-I've slightly improved the 1st version of several programs.

-The 2nd version of some programs calculate and store or recall the partial derivatives gij / xkxl  in a DATA file ( in X-Memory )
-Thus, the initialization is slower, but many routines - like RIJKL  RIJ  RR - become much faster.

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°) A Test involving F00 + GIJ name in Alpha register

-The M-code routine Fµ compares the contents of X- and Y-registers i & j
-If i # j and if flag F00 is set ( diagonal metric tensor ) the next line in program memory is executed
-Otherwise, it is skipped.

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

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

 FS? 00   X=Y?   X=0?   GTO 00   X>Y?   X<>Y   "G"   ARCL X   ARCL Y

are replaced by

 Fµ   GTO 00

0CC  "µ"
006   "F"
04E  C=0
168  M=C
1A8  N=C
1E8  O=C
228  P=C
0B8  C=Y
10E  A=C ALL
0F8  C=X
30E  ?A<C ALL
017  JC+02
0AE  A<>C
06E  A<>B ALL
0AE  A<>C ALL
04E  C=0
09C  PT=5
110  LD@PT- 4
1D0  LD@PT- 7
0D0  LD@PT- 3
010  LD@PT- 0
0D0  LD@PT- 3
0AE  A<>C ALL
37C  RCR 12
102  A=C @PT
21C  PT=2
0EE  B<>C ALL
0FC  RCR 10
102  A=C @PT
0AE  A<>C ALL
168  M=C
3B8  C=d
10E  A=C ALL
2DC  PT=13
210   LD@PT- 8
3B0  C=AandC
2FE  ?C#0 MS
033  JNC+06
0F8  C=X
10E  A=C ALL
0B8  C=Y
36E  ?A#C ALL
360  ?C RTN
141   PC
0A4
3E5   =
0A8
0BD
08E   PC+1 line

( 49 words )

-This function is only used in the "Riemannian programs"

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' |                                ( Registers R01 thru R03 are to be initialized before executing "KHT" )

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

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

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

"X"  ASTO 01
"Y"  ASTO 02
"Z"  ASTO 03

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 02  68 LBL 01  69 RCL 04  70 ST- 06  71 RCL 03  72 RCL 06  73 +  74 RCL 02  75 RCL 01  76 XEQ IND 00  77 STO 07  78 RCL 03  79 RCL 06  80 -  81 RCL 02  82 RCL 01  83 XEQ IND 00  84 ST+ 07  85 RCL 03  86 RCL 02  87 RCL 06  88 +  89 RCL 01  90 XEQ IND 00  91 ST+ 07  92 RCL 03  93 RCL 02  94 RCL 06  95 -  96 RCL 01 97 XEQ IND 00  98 ST+ 07  99 RCL 03 100 RCL 02 101 RCL 01 102 RCL 06 103 + 104 XEQ IND 00 105 ST+ 07 106 RCL 03 107 RCL 02 108 RCL 01 109 RCL 06 110 - 111 XEQ IND 00 112 ST+ 07 113 RCL 03 114 RCL 02 115 RCL 06 116 ST+ Z 117 + 118 RCL 01 119 XEQ IND 00 120 STO 08 121 RCL 03 122 RCL 02 123 RCL 06 124 ST+ Z 125 - 126 RCL 01 127 XEQ IND 00 128 ST+ 08 129 RCL 03 130 RCL 02 131 RCL 06 132 ST- Z 133 + 134 RCL 01 135 XEQ IND 00 136 ST+ 08 137 RCL 03 138 RCL 02 139 RCL 06 140 ST- Z 141 - 142 RCL 01 143 XEQ IND 00 144 ST+ 08 145 RCL 03 146 RCL 02 147 RCL 01 148 RCL 06 149 ST+ T 150 + 151 XEQ IND 00 152 ST+ 08 153 RCL 03 154 RCL 02 155 RCL 01 156 RCL 06 157 ST+ T 158 - 159 XEQ IND 00 160 ST+ 08 161 RCL 03 162 RCL 02 163 RCL 01 164 RCL 06 165 ST- T 166 + 167 XEQ IND 00 168 ST+ 08 169 RCL 03 170 RCL 02 171 RCL 01 172 RCL 06 173 ST- T 174 - 175 XEQ IND 00 176 ST+ 08 177 RCL 03 178 RCL 02 179 RCL 01 180 RCL 06 181 ST+ Z 182 + 183 XEQ IND 00 184 ST+ 08 185 RCL 03 186 RCL 02 187 RCL 01 188 RCL 06 189 ST+ Z 190 - 191 XEQ IND 00 192 ST+ 08 193 RCL 03 194 RCL 02 195 RCL 01 196 RCL 06 197 ST- Z 198 + 199 XEQ IND 00 200 ST+ 08 201 RCL 03 202 RCL 02 203 RCL 01 204 RCL 06 205 ST- Z 206 - 207 XEQ IND 00 208 RCL 08 209 + 210 RCL 09 211 ST- 07 212 ST+ X 213 - 214 RTN 215 LBL 02 216 RCL 05 217 + 218 180 219 RCL 04 220 X^2 221 X^2 222 * 223 / 224 STO 05 225 END

( 331 bytes / SIZE 009 )

 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.34278111 = R05             ---Execution time = 1m39s---

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

-With h = 0.2  it yields   -14.34286490

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 )

-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                                ( Registers R01-R02-R03 are to be initialized before executing "CDGL" )

•  R01 = f name         R04 = x    R07-R08-R09 = Curl E    R11-R12-R13 = Grad (f)    R20 = Lap (f)
•  R02 = g name        R05 = y    R10 = Div E                      R14-R15-R16 = Grad (g)   R21 = Lap (g)
•  R03 = h name        R06 = z                                             R17-R18-R19 = Grad (h)    R22 = Lap (h)

R23 to R38 contain the same results as R07 to R22

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

( 442 bytes / SIZE 039 )

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

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

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

U  ASTO 01
V  ASTO 02
W  ASTO 03     CF 01   CF 02

-If you choose  h = 0.1

0.1  ENTER^
3    ENTER^
2    ENTER^
1    XEQ "CDGL"  >>>>     8.619381890  = R07                                                                     ---Execution time = 127s---
RDN    -32.56682777  = R08
RDN     71.78978318  = R09
RDN     45.44140669  = R10

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

and    Div E = 45.44140669

-You also find in registers R11 to R22:

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 "U"  02  RDN  03  SIN  04  *  05  X^2  06  R^  07  *  08  RTN  09  LBL "V"  10  X^2  11  RCL Z  12  *  13  RTN  14  LBL "W"  15  X^2  16  LASTX  17  *  18  X<>Y  19  COS  20  *  21  *  22  RTN

U  ASTO 01
V  ASTO 02
W  ASTO 03              SF 01   CF 02

-If you choose  h = 0.1

0.1   ENTER^
1
PI  5  /
2    R/S   >>>>    -6.351141008  = R07                                                                     ---Execution time = 135s---
RDN    -8.326237924  = R08
RDN     5.048943483  = R09
RDN     7.163118958  = R10

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

and    Div E = 7.163118958

-You also find in registers R11 to R22:

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

U  ASTO 01
V  ASTO 02
W  ASTO 03                 SF 02

-If you choose  h = 0.1

0.1
PI  5  /
PI  3  /
2    R/S   >>>>    -3.379867347  = R07                                                                     ---Execution time = 187s---
RDN    -6.059707086  = R08
RDN     2.959890531  = R09
RDN    -0.045010875  = R10

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

and    Div E = -0.045010875

-You also find in registers R11 to R22:

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

Notes:

-There are a few comments in the listing itself
-So you can take a look to get some info if need be:

 01  LBL "CDGL"  02  GTO 01  03  "SF1=CYL"  04  "SF2=SPH"  05  "CDGL=7-10-11-20  06  LBL 01  07  .....

Curl begins at R07
Div E =  R10
Grad  are in R11 to R19
Laplacians begin at R20

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

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

( 333 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) ]

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

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:   F00                       CF 00 = general case
SF 00 = diagonal metric

Subroutines:   The n(n+1)/2 functions gij if CF 00    ( i <= j )          and     "LS"  CIJK  "FD"   Fµ
The n functions gii            if  SF 00

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

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

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

g11 = 1 + x2 y2z2
g22 = 1 + x2 + y2 + z2                            SET FLAG F00
g33 = 1 + x2 y2 + x2 z2 + y2 z2

-Load for instance the 3 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  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 = 3m52s---   ( 9mn03s with the 2nd version below )

-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 2m15s instead of  3m55s

-There are a few comments in the listing itself

 01  LBL "INIGC"  02  "H^N"  03  SF0=DIAG  04  ...........

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

Second Version

-This version also stores the partial derivatives gij / xkxl  in a DATA file ( in X-Memory ) named "D2G"

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

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

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

-Same instructions, but execution time = 9m03s

b) Recalling Gij & Cijk + Another M-Code Routine

-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

-To recall the partial derivatives stored in "D2G" datafile, another routine may also be used.
-In the 2nd module, it's "hidden" in the header -RIEMANNIAN

-RIEMANNIAN  takes  i , j , k , l  in the stack, calculates the position N of  kl gij in the file and executes SEEKPT

-The formula is  N = ( n - i ) + A + B [ ( n2 - n - l2 + l ) / 2 + ( n - k ) }

with   A = 0   if  F00  is set ( diagonal metric )  or  A = ( n2 - n - j2 + j ) / 2    if  F00 is clear ( general case )
and    B = n   if  F00  is set  ------------------  or  B = n ( n + 1 ) / 2             if  F00 is clear  ---------------

08E  "N"
001  "A"
009   "I"
00E  "N"
00E  "N"
001  "A"                 This routine is not needed if you use the first version of "INIGC"
00D  "M"
005  "E"
009  "I"
012  "R"
02D  "-"
248   SETF 9
033   JNC+06
08B  "K"
00A  "J"
009   "I"
003   "C"
244   CLRF 9
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
0EE  B<>C ALL
130  LDI S&X
01E  1Eh=30d
246  C=A-C S&X
270  RAMSLCT
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
24C ?FSET 9
1E7  JC+3Ch
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
10E  A=C ALL
046  C=0 S&X
270  RAMSLCT
0AE  A<>C ALL
0E8  X=C
3E0  RTN
0B8  C=Y
10E  A=C ALL
0F8  C=X
30E  ?A<C ALL
027  JC+04
0A8  Y=C
0AE  A<>C ALL
0E8  X=C
3B8  C=d
10E  A=C ALL
2DC  PT=13
210  LD@PT- 8
3B0  C=AandC
2FE  ?C#0 MS
017  JC+02
244  CLRF 9
278  C=Q
2BE  C=-C
10E  A=C ALL
138  C=L
01D  C=
060  A+C
168  M=C
0F8  C=X
10E  A=C ALL
2BE  C=-C
135  C=
061  A*C
0F8  C=X
025  C=
060  AB+C
04E  C
35C   =
090   2
1A8  N=C
269   C=
060  AB/C
178  C=M
025  C=
060  AB+C
1E8  O=C
278  C=Q
10E  A=C ALL
0B8  C=Y
1CE  A=A-C ALL
1F8  C=O
01D  C=
060  A+C
138  C=L
24C  ?FSET 9
013  JNC+02
278  C=Q
13D  C=
060  AB*C
1E8  O=C
046  C
270  =
038  T
10E  A=C ALL
078  C=Z
30E  ?A<C ALL
01F  JC+03
0AE  A<>C ALL
068  Z=C
278  C=Q
0AE  A<>C ALL
1CE  A=A-C ALL
1F8  C=O
01D C=
060  A+C
0E8  X=C
24C  ?FSET 9
09F  JC+13h=19d
078  C=Z
10E  A=C ALL
2BE  C=-C
135   C=
061  A*C
078  C=Z
025  C=
060  AB+C
1B8  C=N
269  C=
060  AB/C
178  C=M
025  C=
060  AB+C
0F8  C=X
025  C=
060  AB+C
0E8  X=C
345  ?NCXQ
040  CLA
260  SETHEX
0B1  ?NCGO
0FE  SEEKPT

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

 01 LBL "RIJKL"  02 STO 24  03 RDN  04 STO 23   05 RDN  06 STO 22   07 X<>Y  08 STO 21         09 RCL 09  10 STO 25  11 CLX  12 STO 28   13 LBL 01  14 RCL 09  15 FS? 00  16 RCL 25  17 STO 26  18 LBL 02  19 RCL 25  20 RCL 26  21 ENTER  22 SIGN  23 CHS  24 CIJK 25 X=0?  26 GTO 00  27 STO 29  28 RCL 22  29 RCL 23   30 RCL 25  31 CIJK  32 STO 27         33 RCL 21  34 RCL 24  35 RCL 26   36 CIJK  37 ST* 27  38 RCL 22  39 RCL 24  40 RCL 25  41 CIJK  42 STO 19  43 RCL 21  44 RCL 23  45 RCL 26  46 CIJK  47 ST* 19  48 RCL 27 49 RCL 19  50 -  51 RCL 29   52 *  53 ST+ 28  54 LBL 00  55 FS? 00  56 GTO 00  57 DSE 26  58 GTO 02   59 LBL 00   60 DSE 25  61 GTO 01  62 RCL 21  63 RCL 24   64 Fµ  65 GTO 00  66 ASTO 00  67 RCL 10  68 RCL 22  69 RCL 23  70 XROM "SD"  71 ST+ 25  72 LBL 00 73 RCL 22  74 RCL 23   75 Fµ  76 GTO 00  77 ASTO 00  78 RCL 10  79 RCL 21  80 RCL 24   81 XROM "SD"  82 ST+ 25  83 LBL 00  84 RCL 21   85 RCL 23         86 Fµ  87 GTO 00  88 ASTO 00  89 RCL 10  90 RCL 22  91 RCL 24  92 XROM "SD"  93 ST- 25  94 LBL 00  95 RCL 22  96 RCL 24 97 Fµ  98 GTO 00  99 ASTO 00 100 RCL 10 101 RCL 21  102 RCL 23  103 XROM "SD" 104 ST- 25 105 LBL 00 106 RCL 25  107 2 108 / 109 RCL 28 110 + 111 SIGN 112 RCL 21 113 RCL 22 114 RCL 23 115 RCL 24 116 X<> L 117 END

( 214 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 = 25s---   (  or  16s  if  h < 0 )

-Likewise

2  ENTER^
3  ENTER^
1  ENTER^
3     R/S      >>>>  R2313 = 1/2                                      ---Execution time = 34s---   (  or  18s  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

Second Version

-The advantage of the 2nd version of "INIGC" is that we just have to recall the partial derivatives in the "D2G" data file.

 01 LBL "RIJKL"  02 STO 14  03 RDN  04 STO 13   05 RDN  06 STO 12   07 X<>Y  08 STO 11  09 RCL 09  10 STO 15                11 CLX  12 STO 18  13 LBL 01   14 RCL 09  15 FS? 00  16 RCL 15  17 STO 16  18 LBL 02  19 RCL 15  20 RCL 16  21 ENTER  22 SIGN  23 CHS 24 CIJK  25 X=0?  26 GTO 00  27 STO 19  28 RCL 12   29 RCL 13   30 RCL 15  31 CIJK  32 STO 17  33 RCL 11  34 RCL 14  35 RCL 16                36 CIJK  37 ST* 17  38 RCL 12            39 RCL 14  40 RCL 15  41 CIJK  42 STO 20  43 RCL 11  44 RCL 13  45 RCL 16  46 CIJK 47 ST* 20  48 RCL 17  49 RCL 20  50 -  51 RCL 19   52 *  53 ST+ 18  54 LBL 00  55 FS? 00  56 GTO 00  57 DSE 16  58 GTO 02  59 LBL 00  60 DSE 15  61 GTO 01  62 RCL 11  63 RCL 14  64 Fµ  65 GTO 00  66 RCL 12  67 RCL 13  68 -RIEMANNIAN  69 GETX 70 ST+ 15  71 LBL 00  72 RCL 12  73 RCL 13  74 Fµ  75 GTO 00  76 RCL 11  77 RCL 14  78 -RIEMANNIAN  79 GETX  80 ST+ 15  81 LBL 00  82 RCL 11   83 RCL 13               84 Fµ  85 GTO 00  86 RCL 12  87 RCL 14  88 -RIEMANNIAN  89 GETX  90 ST- 15  91 LBL 00  92 RCL 12 93 RCL 14  94 Fµ  95 GTO 00  96 RCL 11  97 RCL 13   98 -RIEMANNIAN  99 GETX 100 ST- 15 101 LBL 00 102 RCL 15 103 2 104 / 105 RCL 18 106 + 107 SIGN 108 RCL 11 109 RCL 12 110 RCL 13 111 RCL 14 112 X<> L 113 END

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

-Same instructions, but execution time = 7 seconds instead of 26s

1-time contravariant 3-times covariant curvature tensor

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

Rijkl = gim Rjmkl

-"RI^JKL" calculates these components

 01 LBL "RI^JKL"  02 STO 24  03 RDN  04 STO 23               05 RDN  06 STO 21   07 X<>Y  08 STO 30  09 RCL 09 10 FS? 00  11 X<>Y  12 STO 22               13 CLX  14 STO 20   15 LBL 03   16 RCL 22  17 RCL 30  18 ENTER 19 CLX  20 CIJK  21 X=0?  22 GTO 00  23 STO 31   24 RCL 21               25 RCL 22   26 RCL 23  27 RCL 24 28 XROM "RIJKL"  29 RCL 31               30 *  31 ST+ 20  32 LBL 00  33 FS? 00  34 GTO 00  35 DSE 22  36 GTO 03 37 LBL 00  38 RCL 24               39 SIGN  40 RCL 30   41 RCL 21  42 RCL 23  43 RCL 20  44 END

( 83 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 "RI^JKL"  >>>>  R1212 = 3/8                                      ---Execution time = 26s---   (  or  17s  if  h < 0 )

-Likewise

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

Second Version

-Here, we can replace R21 by R14 and so on , thus saving a few bytes.

 01 LBL "RI^JKL"  02 STO 14  03 RDN  04 STO 13               05 RDN  06 STO 11   07 X<>Y  08 STO 21  09 RCL 09 10 FS? 00  11 X<>Y  12 STO 12   13 CLX  14 STO 22               15 LBL 03  16 RCL 12  17 RCL 21  18 ENTER 19 CLX  20 CIJK  21 X=0?  22 GTO 00  23 STO 23               24 RCL 11   25 RCL 12  26 RCL 13  27 RCL 14 28 XROM "RIJKL"  29 RCL 23               30 *  31 ST+ 22  32 LBL 00  33 FS? 00  34 GTO 00  35 DSE 12  36 GTO 03 37 LBL 00  38 RCL 14               39 SIGN  40 RCL 21  41 RCL 11   42 RCL 13  43 RCL 22  44 END

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

-Same instructions but for example, R2313 is calculated in 8 seconds instead of 35.

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 )

 01 LBL "RIJ"  02 STO 23  03 X<>Y  04 STO 21              05 RCL 09   06 STO 22   07 CLX  08 STO 20  09 LBL 01 10 RCL 09  11 FS? 00  12 RCL 22  13 STO 24              14 LBL 02   15 RCL 22  16 RCL 24  17 ENTER  18 CLX 19 CIJK  20 X=0?  21 GTO 00  22 STO 30  23 RCL 21  24 RCL 22  25 RCL 23   26 RCL 24              27 XROM "RIJKL" 28 RCL 30   29 *  30 ST+ 20  31 LBL 00  32 FS? 00  33 GTO 00   34 DSE 24   35 GTO 02              36 LBL 00 37 DSE 22  38 GTO 01  39 RCL 23              40 SIGN  41 RCL 21   42 RCL 20  43 END

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

 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 = 1m37s---   (  or  59s  if  h < 0 )

-Likewise

2  ENTER^
3     R/S      >>>>  R23 = -3/8                                      ---Execution time = 2m20s---   (  or  75s  if  h < 0 )

-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

Second Version

 01 LBL "RIJ"  02 STO 13               03 X<>Y  04 STO 11   05 RCL 09   06 STO 12  07 CLX  08 STO 21  09 LBL 01 10 RCL 09  11 FS? 00  12 RCL 12               13 STO 14   14 LBL 02  15 RCL 12  16 RCL 14  17 ENTER  18 CLX 19 CIJK  20 X=0?  21 GTO 00  22 STO 22     23 RCL 11   24 RCL 12  25 RCL 13  26 RCL 14  27 XROM "RIJKL" 28 RCL 22               29 *  30 ST+ 21  31 LBL 00  32 FS? 00  33 GTO 00  34 DSE 14  35 GTO 02   36 LBL 00 37 DSE 12  38 GTO 01  39 RCL 13               40 SIGN   41 RCL 11   42 RCL 21  43 END

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

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

-Same instructions, but  R23 is calculated in 23 seconds instead of 2m20s.

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

( 304 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 = 5m06s---   (  or  3m06s  if  h < 0 )

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

Second Version

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

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

 STACK INPUT OUTPUT X / R

-With the same example, R is computed in 75 seconds instead of 5m06s

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"  >>>>   "XEQ RR FIRST"  and finally     E11 = 7/32                                      ---Execution time = 1m37s---   (  or  59s  if  h < 0 )

-Likewise

2  ENTER^
3     R/S      >>>>  E23 = -3/8                                      ---Execution time = 2m20s---  or  24 seconds with the 2nd version !

-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 "XEQ RR FIRST" 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"  >>>>  "XEQ RR FIRST"  and finally    W1212 = 0                    ---Execution time = 3m42s---   (  or  2m17s  if  h < 0 )

-Likewise

2  ENTER^
3  ENTER^
1  ENTER^
3     R/S      >>>>  W2313 = 0                                      ---Execution time = 2m50s---    or  45s with the 2nd version !

-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 "XEQ RR FIRST" 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"

 01 LBL "DKTIJ"  02 FS? 30  03 "SF1=^ SF2=^^"  04 STO 23               05 RDN  06 X<>Y  07 FRC  08 CF 10  09 X#0?  10 SF 10  11 X<> L  12 X<>Y  13 FS? 10  14 1  15 STO 22  16 RDN  17 STO 21  18 X<>Y  19 STO 19  20 CLX  21 STO 24  22 RCL 09  23 STO 25 24 FS? 10  25 GTO 02  26 LBL 01  27 RCL 23   28 RCL 22               29 RCL 25   30 FS? 02  31 X<>Y  32 CIJK  33 X=0?  34 GTO 01  35 STO 20  36 RCL 21  37 RCL 25  38 RCL 09  39 ST* Y  40 -  41 +  42 ENTER  43 SIGN  44 -  45 RCL 19  46 + 47 RCL IND X  48 XEQ IND X  49 RCL 20   50 *  51 FC? 02  52 CHS  53 ST+ 24  54 LBL 01  55 DSE 25   56 GTO 01  57 RCL 09               58 STO 25   59 LBL 02  60 RCL 23  61 RCL 21  62 RCL 25  63 FC? 01  64 FS? 02  65 X<>Y  66 CIJK  67 X=0?  68 GTO 02  69 STO 20 70 RCL 25  71 RCL 22  72 RCL 09   73 ST* Y  74 -  75 +  76 RCL 19               77 +  78 ENTER  79 SIGN  80 -  81 RCL IND X  82 XEQ IND X  83 RCL 20   84 *  85 FC? 01  86 FS? 02  87 CHS  88 ST- 24  89 LBL 02  90 DSE 25  91 GTO 02  92 RCL 21 93 RCL 22  94 RCL 09   95 ST* Y  96 -  97 +  98 RCL 19               99 + 100 1 101 - 102 RCL IND X 103 STO 00 104 RCL 10 105 RCL 23  106 XROM "FD" 107 ST+ 24 108 RCL 23  109 SIGN 110 RCL 19 111 RCL 21 112 FC? 10 113 RCL 22 114 RCL 24 115 END

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

Note:

-There is some info in the listing itself:

 01  LBL "DKTIJ"  02  FS? 30  03  SF1=^ SF2=^^  04  ..................

-Thus set flag F01 if there is only 1 upstairs index and set flag 2 if there are 2 upstairs index
-Clear these flags otherwise

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"

-This listing also contains the programs in the next paragraphs

 01 LBL "RCURL"  02 XEQ 14  03 ST- 26  04 LBL 06  05 RCL 21   06 STO 22  07 ISG 22  08 LBL 07  09 RCL 26              10 RCL 22  11 +  12 RCL IND X  13 STO 00  14 RCL 10  15 RCL 21  16 XROM "FD"  17 STO 25  18 RCL 26   19 RCL 21   20 +  21 RCL IND X  22 STO 00  23 RCL 10  24 RCL 22  25 XROM "FD"  26 RCL 25  27 -  28 STO IND 23  29 CLX  30 SIGN  31 ST+ 22  32 ST+ 23  33 RCL 09  34 RCL 22  35 X<=Y?  36 GTO 07 37 SIGN  38 ST+ 21  39 CLX  40 RCL 21   41 XV^" 178 FS? 30 179 "CF1=>V^" 180 XEQ 14 181 LBL 11 182 CLX 183 STO 14 184 SIGN 185 STO 12 186 LBL 12 187 RCL 21  188 RCL 12            189 1 190 CHS 191 FC? 01 192 CLX 193 CIJK 194 RCL IND 26 195 * 196 ST+ 14 197 CLX 198 SIGN 199 ST+ 26 200 ST+ 12 201 RCL 09  202 RCL 12  203 X<=Y? 204 GTO 12 205 RCL 14  206 STO IND 23 207 RCL 09 208 ST- 26 209 SIGN 210 ST+ 21 211 ST+ 23 212 LASTX 213 RCL 21 214 X<=Y? 215 GTO 11 216 XEQ 13 217 RTN 218 LBL "V*W^" 219 STO M 220 CLX 221 LBL 10 222 RCL IND Y 223 RCL IND M 224 * 225 + 226 ISG M 227 TEXT0 228 ISG Y 229 GTO 10            230 0 231 X<> M 232 X<>Y 233 RTN 234 LBL 14 235 STO 19  236 INT 237 STO 26  238 RCL 39 239 RCL 09 240 ST+ Z 241 2 242 + 243 * 244 41 245 + 246 XY 248 STO 23 249 STO 24 250 SIGN 251 STO 21 252 END

( 443 bytes / SIZE var. )

 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

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 )

-There is some info in the listing about flag F01:

 189  LBL "V<>V^"  190  FS? 30  191  "CF1=>V^"

-Flag F30 is used to skip the next line when the program is executed, so the alpha "register" is unchanged.

k) Dot Product of a Covariant & a Contravariant Vector

-"V*W^"  calculates  vk wk = v1 w1 + v2 w2 + ............... + vn wn

-Store the n components of each vector in contiguous registers

 STACK INPUTS OUTPUTS Y bbb.eee(V) / X bbb.eee(W^) vk wk

Example:   V = ( 2 , 3 , 4 )  &  W^ = ( 7 , 2 , 8 )

-If you store V into R31-R32-R33 and W^ into R71-R72-R73

31.033  ENTER^
71.073  XEQ "V*W^"  >>>>    52                                              --Execution time = 1s---

Notes:

-If V & W are covariant, execute "V<>V^" - CF 01 - with V or W ( but not both ) before executing "V*W^"
-If V & W are contravariant, execute "V<>V^" - SF 01 -with V or W ( but not both ) before executing "V*W^"

-We have   V.W = vk wk = gij vi wj = gij vi wj

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

CF 00

0.1  ENTER^
3    XEQ "INIGC"   >>>>   41.070                                                               ---Execution time = 20m59s---     ( 30m14s with the 2nd module )

-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 = 98s---   ( 18s with the 2nd version )

1  ENTER^
2  ENTER^
1  ENTER^
2  XEQ "RI^JKL"  >>>>  R1212 = 13/96                                     ---Execution time = 269s---   ( 58s with the 2nd version )

Ricci Tensor:

1  ENTER^
3  XEQ "RIJ"      >>>>  R13 = -17/96                                         ---Execution time = 14m12s---   ( 2m54s with the 2nd version )

-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             ( probably the slowest program )      ---Execution time = 68mn---   ( 17m30s with the 2nd version )

Einstein Tensor:   ( after executing "RR" )

2  ENTER^
3  XEQ "EIJ"      >>>>  "XEQ RR FIRST"  and   E23 = -299/288

-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"  >>>>  "XEQ RR FIRST"  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

 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

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

SF 00       ( since this metric is orthogonal ).

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

0.1   ENTER^
4    XEQ INIGC    >>> >   41.100                                            ---Execution time = 6m25s---   ( 17m24s with the 2nd version )

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 )

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

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

Ricci Tensor:

2  ENTER^
3  XEQ "RIJ"      >>>>  R23 = 0

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

Scalar Curvature:

XEQ "RR"  >>>>   R =  -0.000000846    ( exact result = 0 )                       ---Execution time = 8m32s---   ( 2m40s with the 2nd version )

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"  >>>>  "XEQ RR FIRST"  and finally    W1212 = -0.500000095   ( exact result = -1/2 )              ---Execution time = 5m10s---
but  90s with the 2nd version

2  ENTER^
3  ENTER^
2  ENTER^
3     R/S     >>>>      W2323 = 0.500000096   ( 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 ...

Notes:

-Though it's not very fast, the execution times are more "acceptable" when the metric is diagonal.
-And the second version of "INIGC" makes several routines about 3 or 4 times as fast.

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