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 / ¶xk¶xl
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
/ ¶xk¶xm
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
)
t = 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 / ¶x¶y
)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
/ ¶x¶y
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
/ ¶x¶y
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
/ ¶x¶y
= -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 "zz = ¶2f
/ ¶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
/ ¶x¶y
with the formula:
f "xy = ¶2f
/ ¶x¶y
= (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 / ¶x¶y
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
/ ¶x¶y
= -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 / ¶x2¶y2 + ¶4f / ¶x2¶z2 + ¶4f / ¶y2¶z2 )
-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 / ¶x2¶y2
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
-Load the short routines:
| 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 / ¶x¶y = (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 / ¶xi¶xj
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
/ ¶x¶y
= -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
/ ¶x¶y
= -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) gkm ( ¶i gmj + ¶j gim - ¶m 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
/ ¶xk¶xl
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
038 READATA
0EE B<>C ALL
130 LDI S&X
01E 1Eh=30d
246 C=A-C S&X
270 RAMSLCT
038 READATA
10E A=C ALL
046 C=0 S&X
270 RAMSLCT
0EE B<>C ALL
355 chk A&C
050 and SETDEC
268 Q=C
0AE A<>C ALL
128 L=C
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
038 READATA
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
, gij , Gkij
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 , gij , Gkij
• 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 X<Y? 42 GTO 06 43 LBL 13 44 RCL 23 45 DSE X 46 E3 47 / 48 RCL 24 49 + 50 RCL 19 51 SIGN 52 RDN 53 RTN 54 LBL "RDIV" 55 STO 19 56 INT 57 RCL 09 58 STO 21 59 + 60 STO 25 61 CLX 62 STO 24 63 DSE 25 64 LBL 04 65 RCL 09 66 STO 22 67 CLX 68 STO 23 69 LBL 05 70 RCL 21 71 RCL 22 72 ENTER |
73 CIJK
74 ST+ 23 75 DSE 22 76 GTO 05 77 RCL IND 25 78 STO 00 79 XEQ IND 00 80 ST* 23 81 RCL 10 82 RCL 21 83 XROM "FD" 84 RCL 23 85 + 86 ST+ 24 87 DSE 25 88 DSE 21 89 GTO 04 90 RCL 19 91 SIGN 92 RCL 24 93 RTN 94 LBL "RLAPGR" 95 CLX 96 STO 24 97 RCL 39 98 RCL 09 99 STO 22 100 STO 23 101 2 102 + 103 * 104 41 105 + 106 STO 26 107 ST+ 23 108 DSE 23 |
109 LBL 08
110 RCL 10 111 RCL 22 112 XROM "FD" 113 STO IND 23 114 DSE 23 115 DSE 22 116 GTO 08 117 RCL 09 118 STO 21 119 ST+ 23 120 LBL 01 121 RCL 21 122 STO 22 123 LBL 02 124 CLX 125 STO 20 126 RCL 09 127 STO 25 128 LBL 03 129 RCL 21 130 RCL 22 131 RCL 25 132 CIJK 133 X=0? 134 GTO 00 135 RCL IND 23 136 * 137 ST+ 20 138 LBL 00 139 DSE 23 140 DSE 25 141 GTO 03 142 RCL 09 143 ST+ 23 144 RCL 10 |
145 RCL 21
146 RCL 22 147 XROM "SD" 148 ST- 20 149 RCL 21 150 RCL 22 151 X=Y? 152 GTO 00 153 RCL 20 154 ST+ 20 155 RDN 156 LBL 00 157 ENTER 158 CLX 159 CIJK 160 RCL 20 161 * 162 ST- 24 163 FS? 00 164 GTO 00 165 DSE 22 166 GTO 02 167 LBL 00 168 DSE 21 169 GTO 01 170 RCL 23 171 E3 172 / 173 RCL 26 174 + 175 RCL 24 176 RTN 177 LBL "V<>V^" 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 X<Y? 247 X<>Y 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 , gij
, Gkij
• 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 , gij , Gkij
• 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 , gij , Gkij
• 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:
XEQ "RAD"
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