Differential Geometry(2) for the HP41
Overview
0°) 2 M-Code Routines
1°) Curvature & Torsion of a Parametric
Curve
2°) Gaussian Curvature of a Surface
3°) Biharmonic Operator
4°) Curl, Divergence, Gradient &
Laplacian
5°) Laplacian of a function of n variables
( n < 10 )
6°) Linear Systems
7°) Riemannian Geometry
a) Initialization
b) Recalling Gij & Cij^k
c) Curvature Tensor
d) Ricci Tensor
e) Scalar Curvature
f) Einstein Tensor
g) Weyl Tensor
h) Covariant Derivative of
a Vector or a Tensor of order 2
i) Curl, Divergence, Laplacian
& Gradient
j) Covariant Vector <>
Contravariant Vector
8°) An Example of a Non-Diagonal Metric
9°) Schwarzchild Metric
10°) Non-Riemannian Geometry
a) Curvature Tensor(s)
b) Ricci Tensor & Segmental
Curvature
c) Scalar Curvature(s)
d) Torsion Tensor
e) Torsion Vector
11°) Differential Geometry QRG
-Almost all these programs require the evaluation of derivatives
-Formulas of order 10 are usually used, but the execution time may
be prohibitive to compute all tensors..
-So, you can also choose 4th-order methods, to get faster - but less
accurate - results.
-However, a good emulator like V41 or a 41CL is recommended.
-Several programs come from "Numerical Differentiation for the HP41",
sometime with minor modifications...
-All are included in a module published here.
>>> Last update:
-I've slightly modified several programs.
-Instead of using F00 for diagonal metrics, I've used a trick to get
the derivative faster when the function = 0
-Simply store an alpha string in X-register and "FD" ( if CF 03 ) or
"SD" will return 0 without evaluating the function 10 times.
-Several routines have been added to calculate a few tensors in non-Riemannian
manifolds.
Notations:
-In the comments below, I've used Einstein summation convention:
-When an index appears as an upper index & a lower index in a monomial
and is not otherwise defined,
it means summation of that term over all the values of the index.
-For example, in a 3 dimensional space,
ak bk = a1 b1 + a2 b2 + a3 b3
gij ai
aj = g11 a1 a1 + g12
a1 a2 + g13 a1 a3
+ g21 a2 a1 + g22 a2
a2 + g23 a2 a3
+ g31 a3 a1 + g32 a3
a2 + g33 a3 a3
-If the tensor gij is symmetric ( gij = gji ) the last expression may be simplified and becomes:
gij ai
aj = g11 a1 a1 + g22
a2 a2 + g33 a3 a3
+ 2 ( g12 a1 a2 + g13 a1
a3 + g23 a2 a3 )
-The partial derivative ¶m
gij means ¶gij
/ ¶xm
-Likewise, ¶km
gij = ¶gij
/ ¶xk¶xm
0°) GIJ name & CIJK name in Alpha register
-The M-code routine Fµ places "GIJ" in the alpha register
where I = Inf ( i , j ) & J = Sup ( i , j )
if X-register = i and Y-register = j ,
-So, the program lines ( with FIX 0 CF 29 ):
| X>Y?
X<>Y "G" ARCL X ARCL Y |
are replaced by
| Fµ |
-A similar function ( called "DIFF MANIF" - a section header in the
module )
places CIJK in alpha provided Z = i Y = j X = k
086 "F"
009 "I"
00E "N"
001 "A"
00D "M"
020 " "
006 "F"
006 "F"
009 "I"
004 "D"
02D "-"
288 SETF 7
023 JNC +04
0CC "µ"
006 "F"
284 CLRF 7
04E C=0
168 M=C
1A8 N=C
1E8 O=C
228 P=C
28C ?FSET 7
033 JNC +06
29C PT=7
110 LD@PT- 4
0D0 LD@PT- 3
0D0 LD@PT- 3
023 JNC+04
09C PT=5
110 LD@PT- 4
1D0 LD@PT- 7
01C PT=3
0D0 LD@PT- 3
010 LD@PT- 0
0D0 LD@PT- 3
28C ?FSET 7
07B JNC+15d
10E A=C ALL
078 C=Z
13C RCR 8
05C PT=4
102 A=C @PT
0B8 C=Y
0FC RCR 10
21C PT=2
102 A=C @PT
0F8 C=X
37C RCR 12
39C PT=0
102 A=C @PT
083 JNC+16d
0EE B<>C ALL
0B8 C=Y
10E A=C ALL
0F8 C=X
30E ?A<C ALL
017 JC+02
0AE A<>C ALL
06E A<>B ALL
37C RCR 12
39C PT=0
102 A=C @PT
0EE B<>C ALL
21C PT=2
0FC RCR 10
102 A=C @PT
0AE A<>C ALL
168 M=C
3E0 RTN
( 69 words )
1°) Curvature & Torsion of a Parametric Curve
-A curve may be defined by 3 functions X(t) Y(t) Z(t) in
a 3-dimensional Euclidean space.
-"KHT" returns the curvature c ( khi ) and
the torsion t ( tau ) of the curve at a point
defined by t
-The curvature is calculated by c
= | r' x r'' | / | r' |3
where the vector r = [ X(t) , Y(t) , Z(t) ]
and the torsion by t = ( r'
x r'' ) .r''' / | r' x r'' |2
x = cross-product and . = dot-product
Data Registers: R00 = | r' |
R04 = z' R07 = y'
R10 = | r' x
r'' |2
R13 = x' R16 = x
R05 = z'' R08 = y''
R11 = t
R14 = x''
R06 = z''' R09 = y'''
R12 = h
R15 = x'''
Flags: /
Subroutines: "dF" and
3 functions named "X" "Y" "Z" that compute
x(t) , y(t) , z(t) assuming t is in X-register upon entry
| 01 LBL "KHT"
02 LBL 13 03 "Z" 04 ASTO 00 05 XROM "dF" 06 STO 12 07 RDN 08 STO 13 09 X<>Y 10 STO 14 11 "Y" 12 ASTO 00 13 RCL 02 14 RCL 01 15 XROM "dF" 16 STO 09 17 RDN 18 STO 10 19 X<>Y 20 STO 11 21 "X" 22 ASTO 00 23 RCL 02 24 RCL 01 25 XROM "dF" 26 RCL 12 27 RCL 09 28 R-P 29 RCL 03 30 R-P 31 STO 00 32 RCL 13 33 RCL 09 34 * 35 RCL 12 36 RCL 10 37 * 38 - 39 X^2 40 STO 06 41 LASTX 42 RCL 05 43 STO 08 44 * 45 RCL 12 46 RCL 04 47 STO 07 48 * 49 RCL 13 50 RCL 03 51 * 52 - 53 ENTER 54 X^2 55 ST+ 06 56 CLX 57 RCL 11 58 * |
59 +
60 RCL 10 61 RCL 03 62 * 63 RCL 09 64 RCL 07 65 * 66 - 67 ENTER 68 X^2 69 ST+ 06 70 CLX 71 RCL 14 72 * 73 + 74 RCL 06 75 ST/ Y 76 SQRT 77 RCL 00 78 3 79 Y^X 80 / 81 X<>Y 82 STO 04 83 X<>Y 84 STO 03 85 RTN 86 GTO 13 87 LBL "dF" 88 LBL 10 89 STO 01 90 X<>Y 91 STO 02 92 X<>Y 93 XEQ IND 00 94 STO 06 95 CLX 96 STO 08 97 XEQ 01 98 STO 04 99 RCL 07 100 STO 03 101 STO 05 102 42 103 ST* 04 104 ST+ X 105 ST* 03 106 70098 107 CHS 108 ST* 05 109 XEQ 01 110 6 111 * 112 ST- 04 113 RCL 07 114 24 115 * 116 ST- 03 |
117 2184.5
118 * 119 ST+ 05 120 XEQ 01 121 ST+ 04 122 RCL 07 123 6 124 * 125 ST+ 03 126 2434.5 127 * 128 ST- 05 129 XEQ 01 130 8 131 / 132 ST- 04 133 RCL 07 134 ST- 03 135 2522 136 * 137 ST+ 05 138 25 139 ST* 03 140 E3 141 ST* 04 142 XEQ 01 143 8 144 * 145 ST+ 04 146 RCL 07 147 ST+ X 148 ST+ 03 149 RCL 07 150 205 151 * 152 ST- 05 153 RCL 02 154 ST/ 04 155 2520 156 * 157 ST/ 03 158 10 159 * 160 ST/ 04 161 1.2 162 * 163 RCL 02 164 X^2 165 * 166 ST/ 05 167 RCL 05 168 RCL 04 169 RCL 03 170 RTN 171 GTO 10 172 LBL 01 173 RCL 02 174 ST+ 08 |
175 RCL 01
176 RCL 08 177 + 178 XEQ IND 00 179 STO 07 180 RCL 01 181 RCL 08 182 - 183 XEQ IND 00 184 RCL 07 185 X<>Y 186 ST- 07 187 + 188 RCL 06 189 ST+ X 190 - 191 RTN 192 LBL 02 193 RCL 14 194 RCL 15 195 X<> Z 196 GTO IND 16 197 LBL 03 198 RCL 15 199 X<>Y 200 RCL 13 201 GTO IND 16 202 LBL 04 203 RCL 14 204 RCL 13 205 GTO IND 16 206 LBL 07 207 RCL 11 208 X<>Y 209 GTO IND 12 210 LBL 08 211 RCL 10 212 GTO IND 12 213 LBL "dF2" 214 LBL 11 215 XROM "MDV" 216 STO 09 217 RCL 01 218 STO 10 219 RCL 02 220 STO 11 221 8 222 X<> 00 223 STO 12 224 RCL 03 225 RCL 11 226 XROM "dF" 227 STO 13 228 RDN 229 STO 14 230 X<>Y 231 STO 15 232 DSE 00 |
233 RCL 02
234 RCL 10 235 XROM "dF" 236 RCL 12 237 STO 00 238 RCL 15 239 STO 08 240 RCL 09 241 SIGN 242 RCL 14 243 STO 07 244 RCL 04 245 RCL 13 246 STO 06 247 RCL 03 248 RTN 249 GTO 11 250 LBL "dF3" 251 LBL 05 252 STO 13 253 RDN 254 STO 14 255 RDN 256 STO 15 257 X<>Y 258 STO 02 259 4 260 X<> 00 261 STO 16 262 RCL 02 263 RCL 15 264 XROM "dF" 265 STO 09 266 RDN 267 STO 10 268 X<>Y 269 STO 11 270 DSE 00 271 RCL 02 272 RCL 14 273 XROM "dF" 274 STO 17 275 RDN 276 STO 18 277 X<>Y 278 STO 19 279 DSE 00 280 RCL 02 281 RCL 13 282 XROM "dF" 283 RCL 16 284 STO 00 285 RCL 19 286 STO 08 287 RCL 06 288 STO 12 289 RCL 04 290 RCL 18 |
291 STO 07
292 + 293 RCL 10 294 + 295 RCL 09 296 RCL 17 297 STO 06 298 RCL 03 299 RTN 300 GTO 05 301 LBL "KGS" 302 LBL 06 303 XROM "dF2" 304 X^2 305 X<>Y 306 X^2 307 + 308 1 309 + 310 X^2 311 X<> Z 312 * 313 RCL 09 314 X^2 315 - 316 X<>Y 317 / 318 RTN 319 GTO 06 320 LBL "MDV" 321 LBL 12 322 STO 01 323 RDN 324 STO 02 325 X<>Y 326 STO 03 327 X<> T 328 XEQ IND 00 329 STO 04 330 CLX 331 STO 06 332 XEQ 01 333 42 334 * 335 STO 05 336 XEQ 01 337 6 338 * 339 ST- 05 340 XEQ 01 341 ST+ 05 342 XEQ 01 343 8 344 / 345 ST- 05 346 E3 347 ST* 05 348 XEQ 01 349 8 |
350 *
351 RCL 05 352 + 353 50400 354 RCL 03 355 X^2 356 * 357 / 358 STO 05 359 RTN 360 GTO 12 361 LBL 01 362 RCL 03 363 ST+ 06 364 RCL 02 365 RCL 01 366 RCL 06 367 + 368 XEQ IND 00 369 STO 07 370 RCL 02 371 RCL 01 372 RCL 06 373 - 374 XEQ IND 00 375 ST+ 07 376 RCL 02 377 RCL 06 378 + 379 RCL 01 380 XEQ IND 00 381 ST+ 07 382 RCL 02 383 RCL 06 384 - 385 RCL 01 386 XEQ IND 00 387 ST+ 07 388 RCL 02 389 RCL 01 390 RCL 06 391 ST+ Z 392 + 393 XEQ IND 00 394 ST- 07 395 RCL 02 396 RCL 01 397 RCL 06 398 ST- Z 399 - 400 XEQ IND 00 401 RCL 07 402 - 403 RCL 04 404 ST+ X 405 + 406 RTN 407 GTO 12 408 END |
( 613 bytes
/ SIZE var. )
| STACK | INPUTS | OUTPUTS |
| Y | h | t (t) |
| X | t | c (t) |
Example: The curve defined by
x(t) = et cos t , y(t) = et sin t
, z(t) = Ln(t) ( t expressed
in radians )
-Evaluate the curvature c (t) and the torsion
t
(t) for t = 1.3
-Load the following routines in main memory
| 01 LBL "X"
02 E^X 03 LASTX 04 COS 05 * 06 RTN 07 LBL "Y" 08 E^X 09 LASTX 10 SIN 11 * 12 RTN 13 LBL "Z" 14 LN 15 RTN |
XEQ "RAD" and if you choose h = 0.1
0.1 ENTER^
1.3 XEQ "KHT" >>>>
c
(1.3)
= 0.194807823
---Execution time = 45s---
X<>Y t (1.3)
= 0.123665258
-The exact values are c (1.3) = 0.194807823 & t (1.3) = 0.123665529
Notes:
-If you want to estimate the curvature and torsion with other parameters, simply press h ENTER^ t R/S
-With an helicoidal curve: x(t) = r Cos t , y(t) =
r Sin t , z(t) = a t
where r & a are positive constants,
c = r / ( a2 + r2
)
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 04 68 LBL 01 69 RCL 04 70 ST- 06 71 CLX 72 STO 07 73 STO 08 74 LBL 02 75 RCL 03 76 RCL 06 77 + 78 RCL 02 |
79 RCL 01
80 XEQ IND 00 81 ST+ 07 82 RCL 03 83 RCL 02 84 RCL 06 85 + 86 RCL 01 87 XEQ IND 00 88 ST+ 07 89 RCL 03 90 RCL 02 91 RCL 01 92 RCL 06 93 + 94 XEQ IND 00 95 ST+ 07 96 RCL 06 97 CHS 98 STO 06 99 X<0? 100 GTO 02 101 STO 10 102 LBL 03 103 RCL 03 104 RCL 06 |
105 +
106 RCL 02 107 RCL 10 108 + 109 RCL 01 110 XEQ IND 00 111 ST+ 08 112 RCL 03 113 RCL 06 114 + 115 RCL 02 116 RCL 01 117 RCL 10 118 + 119 XEQ IND 00 120 ST+ 08 121 RCL 03 122 RCL 02 123 RCL 06 124 + 125 RCL 01 126 RCL 10 127 + 128 XEQ IND 00 129 ST+ 08 130 RCL 10 |
131 CHS
132 STO 10 133 X<0? 134 GTO 03 135 RCL 06 136 CHS 137 STO 06 138 X<0? 139 GTO 03 140 RCL 08 141 RCL 09 142 ST- 07 143 ST+ X 144 - 145 RTN 146 LBL 04 147 RCL 05 148 + 149 180 150 RCL 04 151 X^2 152 X^2 153 * 154 / 155 STO 05 156 END |
( 231 bytes / SIZE 010 )
| STACK | INPUTS | OUTPUTS |
| T | h | / |
| Z | z | / |
| Y | y | / |
| X | x | D2 f |
Example: f(x) = exp(-x2).Ln(y2+z)
x = 1 , y = 2 , z = 3
| 01 LBL "T"
02 X^2 03 CHS 04 E^X 05 RDN 06 X^2 07 + 08 LN 09 R^ 10 * 11 RTN |
T ASTO 00
-If you choose h = 0.1
0.1 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "BHRM" >>>>
D2
f = -14.34282056 = R05
---Execution time = 1m41s---
-The exact result is -14.34264116 ( relative error about 2 E-5 )
-With h = 0.2 it yields -14.34285066 and with h = 0.15 we get -14.34265712
Notes:
-91 evaluations of the function are performed (!).
-Though the formula is of order 10 - due to roundoff-errors - 4th derivatives
are difficult to get with a great precision.
-So, there are seldom more than 5 or 6 significant digits.
4°) Curl, Divergence, Gradient & Laplacian
"CDGL" computes the Curl , divergence, gradient and Laplacian of a 3D-Vector Field E= ( f , g , h ) = ( X , Y , Z )
-If the coordinates are rectangular, clear flags F01 & F02
-If the coordinates are cylindrical, set F01 and clear F02
-If the coordinates are spherical, set F02
( if F02 & F01 are set, F01 is automatically cleared )
-In the last 2 cases, the HP41 sets the RAD mode automatically.
Formulae:
• Rectangular Coordinates x , y , z
Curl E = ( ¶h/¶y - ¶g/¶z , ¶f/¶z - ¶h/¶x , ¶g/¶x - ¶f/¶y )
Div E = ¶f/¶x + ¶g/¶y + ¶h/¶z
Grad f = ( ¶f/¶x , ¶f/¶y , ¶f/¶z ) and similar formulae for g & h
Lapl f = ¶2f
/ ¶x2 + ¶2f
/ ¶y2 + ¶2f
/ ¶z2
and similar formulae for g & h
• Cylindrical Coordinates r , f , z
Curl E = [ (1/r) ¶h/¶f - ¶g/¶z , ¶f/¶z - ¶h/¶r , (1/r) ( g + r ¶g/¶r - ¶f/¶f ) ]
Div E = ¶f/¶r + (1/r) f + (1/r) ¶g/¶f + ¶h/¶z
Grad f = ( ¶f/¶r , (1/r) ¶f/¶f , ¶f/¶z ) and similar formulae for g & h
Lapl f = ¶2f
/ ¶r2 + (1/r) ¶f/¶r
+ (1/r2) ¶2f / ¶f2
+ ¶2f / ¶z2
and similar formulae for g & h
• Spherical Coordinates r , q , f
Curl E = [ (1/(r.sin q) ) ( h cos q + sin q ¶h/¶q - ¶g/¶f ) , (1/(r.sin q) ) ¶f/¶f - ¶h/¶r - h / r , (1/r) ( g + r ¶g/¶r - ¶f/¶q ) ]
Div E = ¶f/¶r + (2/r) f + g (cos q)/(r sin q) + (1/r) ¶g/¶q + (1/(r sin q)) ¶h/¶f
Grad f = ( ¶f/¶r , (1/r) ¶f/¶q , (1/(r sin q)) ¶f/¶f ) and similar formulae for g & h
Lapl f = ¶2f
/ ¶r2 + (2/r) ¶f/¶r
+ (1/r2) ¶2f / ¶q2
+ (1/(r2tan q)) ¶f
/ ¶q + (1/(r2sin2q))
¶2f
/ ¶f2
and similar formulae for g & h
Data Registers: R00 = h
R01 = x R04-R05-R06 = Curl
E R08-R09-R10 = Grad (X)
R17 = Lap (X)
R02 = y R07 = Div E
R11-R12-R13 = Grad (Y) R18 = Lap (Y)
R03 = z
R14-R15-R16 = Grad (Z) R19 = Lap (Z)
R23 to R38 contain the same results as R04 to R19
Flags: /
Subroutines: 3 programs
named "X" "Y" "Z" which compute X(x,y,z) , Y(x,y,z) and
Z(x,y,z)
assuming x is in X-register, y is in Y-register and z is in Z-register
upon entry.
XROM may be replaced by XEQ
| 01 LBL "CDGL"
02 FC? 01 03 FS? 02 04 RAD 05 STO 13 06 RDN 07 STO 14 08 RDN 09 STO 15 10 X<>Y 11 STO 02 12 X<>Y 13 R^ 14 R^ 15 "X" 16 ASTO 00 17 XROM "dF3" 18 STO 26 19 STO 27 20 RDN 21 STO 28 22 CHS 23 STO 25 24 RDN 25 STO 29 26 STO 24 27 X<>Y 28 STO 36 29 FS? 02 30 GTO 02 31 FC? 01 32 GTO 00 33 XEQ 04 34 + 35 STO 36 36 RCL 12 37 RCL 13 38 ST/ 25 |
39 ST/ 28
40 / 41 ST+ 26 42 GTO 00 43 LBL 02 44 XEQ 03 45 STO 36 46 RCL 13 47 ST/ 25 48 ST/ 28 49 RCL 14 50 SIN 51 * 52 ST/ 24 53 ST/ 29 54 RCL 12 55 ST+ X 56 RCL 13 57 / 58 ST+ 26 59 LBL 00 60 "Y" 61 ASTO 00 62 RCL 02 63 RCL 15 64 RCL 14 65 RCL 13 66 XROM "dF3" 67 FC? 01 68 FS? 02 69 GTO 01 70 ST+ 25 71 STO 30 72 RDN 73 STO 31 74 ST+ 26 75 RDN 76 STO 32 |
77 CHS
78 STO 23 79 X<>Y 80 STO 37 81 GTO 00 82 LBL 01 83 FS? 02 84 GTO 02 85 XEQ 04 86 + 87 STO 37 88 RCL 03 89 STO 30 90 RCL 06 91 RCL 13 92 / 93 ST+ 26 94 STO 31 95 RCL 09 96 STO 32 97 CHS 98 STO 23 99 RCL 03 100 RCL 12 101 RCL 13 102 / 103 + 104 ST+ 25 105 GTO 00 106 LBL 02 107 XEQ 03 108 STO 37 109 RCL 03 110 STO 30 111 RCL 06 112 RCL 13 113 / 114 STO 31 |
115 RCL 09
116 LASTX 117 RCL 14 118 SIN 119 * 120 / 121 STO 32 122 RCL 12 123 RCL 14 124 TAN 125 / 126 RCL 06 127 + 128 RCL 13 129 / 130 ST+ 26 131 RCL 09 132 CHS 133 RCL 13 134 RCL 14 135 SIN 136 * 137 / 138 STO 23 139 RCL 12 140 RCL 13 141 / 142 RCL 03 143 + 144 ST+ 25 145 LBL 00 146 "Z" 147 ASTO 00 148 RCL 02 149 RCL 15 150 RCL 14 151 RCL 13 152 XROM "dF3" |
153 FC? 01
154 FS? 02 155 GTO 01 156 STO 33 157 ST- 24 158 RDN 159 STO 34 160 ST+ 23 161 RDN 162 STO 35 163 ST+ 26 164 X<>Y 165 STO 38 166 GTO 00 167 LBL 01 168 FS? 02 169 GTO 02 170 XEQ 04 171 + 172 STO 38 173 RCL 03 174 STO 33 175 RCL 06 176 RCL 13 177 / 178 ST+ 23 179 STO 34 180 RCL 09 181 STO 35 182 ST+ 26 183 RCL 03 184 ST- 24 185 GTO 00 186 LBL 02 187 XEQ 03 188 STO 38 189 RCL 03 190 STO 33 |
191 RCL 06
192 RCL 13 193 / 194 STO 34 195 RCL 09 196 LASTX 197 RCL 14 198 SIN 199 * 200 / 201 STO 35 202 ST+ 26 203 RCL 06 204 RCL 12 205 RCL 14 206 TAN 207 / 208 + 209 RCL 13 210 / 211 ST+ 23 212 RCL 03 213 RCL 12 214 RCL 13 215 / 216 + 217 ST- 24 218 GTO 00 219 LBL 03 220 XEQ 04 221 RCL 14 222 SIN 223 X^2 224 / 225 RCL 06 226 RCL 14 227 TAN 228 / |
229 +
230 RCL 13 231 X^2 232 / 233 + 234 RTN 235 LBL 04 236 RCL 07 237 RCL 13 238 / 239 RCL 03 240 FS? 02 241 ST+ X 242 + 243 RCL 13 244 / 245 RCL 04 246 + 247 RCL 10 248 RTN 249 LBL 00 250 RCL 02 251 STO 00 252 RCL 13 253 STO 01 254 RCL 14 255 STO 02 256 RCL 15 257 STO 03 258 23.004016 259 REGMOVE 260 4.019 261 SIGN 262 RCL 07 263 RCL 06 264 RCL 05 265 RCL 04 266 END |
( 411 bytes / SIZE 039 )
| STACK | INPUTS | OUTPUTS |
| T | h | Div E |
| Z | x3 | Curl3 E |
| Y | x2 | Curl2 E |
| X | x1 | Curl1 E |
| L | / | 4.019 |
4.019 = control number of all the results.
Example1 - Rectangular Coordinates: CF 01 CF 02
E = ( f , g , h ) = [ exp(-x2) Ln(y2+z) , x2y2z2 , exp(x) y2z ] With x = 1 , y = 2 , z = 3
-Load the short routines:
| 01 LBL "X"
02 X^2 03 CHS 04 E^X 05 RDN 06 X^2 07 + 08 LN 09 R^ 10 * 11 RTN 12 LBL "Y" 13 * 14 * 15 X^2 16 RTN 17 LBL "Z" 18 E^X 19 RDN 20 X^2 21 * 22 R^ 23 * 24 RTN |
CF 01 CF 02
-If you choose h = 0.1
0.1 ENTER^
3 ENTER^
2 ENTER^
1 XEQ "CDGL" >>>>
8.619381890 = R04
---Execution time = 127s---
RDN -32.56682777 = R05
RDN 71.78978318 = R06
RDN 45.44140669 = R07
-So, Curl E = ( 8.619381890 , -32.56682777 , 71.78978318 )
and Div E = 45.44140669
-You also find in registers R08 to R19:
Grad f = ( -1.431720683 , 0.210216822 , 0.052554204
)
Grad g = ( 72 , 36 , 24 )
Grad h = ( 32.61938197 , 32.61938189 , 10.87312737
)
Lapl f = 1.409197689
Lapl g = 98
Lapl h = 48.92907013
Example2 - Cylindrical Coordinates: SF 01 CF 02
E = ( f , g , h ) = ( r z2 sin2f
, r2 z , r3 z cos f
)
r = 2 , f = PI/5 , z = 1
| 01 LBL "X"
02 RDN 03 SIN 04 * 05 X^2 06 R^ 07 * 08 RTN 09 LBL "Y" 10 X^2 11 RCL Z 12 * 13 RTN 14 LBL "Z" 15 X^2 16 LASTX 17 * 18 X<>Y 19 COS 20 * 21 * 22 RTN |
SF 01 CF 02
-If you choose h = 0.1
0.1 ENTER^
1
PI 5 /
2 R/S
>>>> -6.351141008 = R04
---Execution time = 135s---
RDN -8.326237924 = R05
RDN 5.048943483 = R06
RDN 7.163118958 = R07
-So, Curl E = ( -6.351141008 , -8.326237924 , 5.048943483 )
and Div E = 7.163118958
-You also find in registers R08 to R19:
Grad f = ( 0.345491503 , 0.951056518 , 1.381966009
)
Grad g = ( 4 , 0 , 4 )
Grad h = ( 9.708203933 , -2.351141008 , 6.472135952
)
Lapl f = 1.863728736
Lapl g = 4
Lapl h = 12.94427209
Example2 - Spherical Coordinates: CF
01 SF 02
E = ( f , g , h ) = ( r sin2 q
cos2 f , r2sin f
, r3 cos q cos2 f
)
r = 2 , q = PI/3 , f
= PI/5
| 01 LBL "X"
02 RDN 03 SIN 04 X<>Y 05 COS 06 * 07 X^2 08 R^ 09 * 10 RTN 11 LBL "Y" 12 X^2 13 RCL Z 14 SIN 15 * 16 RTN 17 LBL "Z" 18 X^2 19 LASTX 20 * 21 X<>Y 22 COS 23 * 24 X<>Y 25 COS 26 X^2 27 * 28 RTN |
CF 01 SF 02
-If you choose h = 0.1
0.1
PI 5 /
PI 3 /
2 R/S
>>>> -3.379867347 = R04
---Execution time = 187s---
RDN -6.059707086 = R05
RDN 2.959890531 = R06
RDN -0.045010875 = R08
-So, Curl E = ( -3.379867347 , -6.059707086 , 2.959890531 )
and Div E = -0.045010875
-You also find in registers R08 to R19:
Grad f = ( 0.490881375 , 0.566820985 , -0.823639103
)
Grad g = ( 2.351141010 , 0 , 1.868344717 )
Grad h = ( 3.927050988 , -2.267283945 , -2.196370944
)
Lapl f = 0.018237234
Lapl g = 2.742997604
Lapl h = 5.721039298
5°) Laplacian of a Function of N variables ( n <
10 )
"LAPL" may be used to evaluate the Laplacian of a function of n variables, provided n < 10
"SD" is called as a subroutine( see below ), so the same formulas
are used, of order 10 if h > 0 or of order 4 if h < 0
>>> In both cases if the function returns an alpha data, 0 is returned
without calculating the whole formula
Data Registers: • R00 = Function name ( Registers R00 thru Rnn are to be initialized before executing "LAPL" )
• R01 = x1 , • R02 = x2 , .......... , • Rnn = xn
R10 = h
R12-..........-R20: temp
Flag: F10
Subroutines: "SD" and a program which
computes f(x1, ..... ,xn) assuming
x1 is in R01 , ............ , xn is in Rnn
| 01 LBL "FD"
02 LBL 01 03 STO 13 04 CLX 05 STO 14 06 RCL IND 13 07 STO 11 08 X<> Z 09 STO 10 10 FC? 03 11 XEQ IND 00 12 SIGN 13 X=0? 14 RTN 15 RCL 10 16 X<0? 17 GTO 00 18 XEQ 03 19 14 20 * 21 STO 12 22 XEQ 03 23 4 24 * 25 ST- 12 26 XEQ 03 27 ST+ 12 28 6 29 ST* 12 30 XEQ 03 |
31 ST- 12
32 25 33 ST* 12 34 XEQ 03 35 ST+ X 36 ST+ 12 37 2520 38 LBL 13 39 RCL 10 40 * 41 ST/ 12 42 RCL 11 43 STO IND 13 44 RCL 12 45 RTN 46 GTO 01 47 LBL 00 48 XEQ 03 49 8 50 * 51 STO 12 52 XEQ 03 53 ST- 12 54 12 55 GTO 13 56 LBL 03 57 RCL 10 58 ST+ 14 59 RCL 11 60 RCL 14 |
61 -
62 STO IND 13 63 XEQ IND 00 64 STO 15 65 RCL 11 66 RCL 14 67 + 68 STO IND 13 69 XEQ IND 00 70 RCL 15 71 - 72 RTN 73 LBL "SD" 74 LBL 02 75 CF 10 76 X=Y? 77 SF 10 78 STO 16 79 RDN 80 STO 17 81 CLX 82 STO 15 83 RCL IND T 84 STO 13 85 RCL IND 17 86 STO 14 87 R^ 88 STO 10 89 XEQ IND 00 90 SIGN |
91 X=0?
92 RTN 93 LASTX 94 STO 11 95 RCL 10 96 X<0? 97 GTO 00 98 XEQ 03 99 42 100 * 101 STO 12 102 XEQ 03 103 6 104 * 105 ST- 12 106 XEQ 03 107 ST+ 12 108 XEQ 03 109 8 110 / 111 ST- 12 112 E3 113 ST* 12 114 XEQ 03 115 8 116 * 117 ST+ 12 118 RCL 13 119 STO IND 16 120 25200 |
121 LBL 14
122 FC? 10 123 ST+ X 124 FC?C 10 125 CHS 126 RCL 10 127 X^2 128 * 129 ST/ 12 130 RCL 12 131 RTN 132 GTO 02 133 LBL 00 134 XEQ 03 135 16 136 * 137 STO 12 138 XEQ 03 139 ST- 12 140 RCL 13 141 STO IND 16 142 12 143 GTO 14 144 LBL "LAPL" 145 LBL 04 146 STO 19 147 X<>Y 148 STO 10 149 CLX 150 STO 20 |
151 LBL 05
152 RCL 10 153 RCL 19 154 ENTER 155 XROM "SD" 156 ST+ 20 157 DSE 19 158 GTO 05 159 RCL 20 160 RTN 161 GTO 04 162 LBL 03 163 RCL 10 164 ST+ 15 165 RCL 15 166 RCL 13 167 + 168 STO IND 16 169 XEQ IND 00 170 STO 18 171 RCL 13 172 RCL 15 173 - 174 STO IND 16 175 XEQ IND 00 176 RCL 11 177 ST+ X 178 - 179 ST+ 18 180 RCL 18 |
181 FS? 10
182 RTN 183 RCL 14 184 RCL 15 185 - 186 STO IND 17 187 XEQ IND 00 188 ST- 18 189 RCL 13 190 STO IND 16 191 XEQ IND 00 192 ST+ 18 193 RCL 14 194 RCL 15 195 + 196 STO IND 17 197 XEQ IND 00 198 ST+ 18 199 RCL 15 200 ST+ IND 16 201 XEQ IND 00 202 ST- 18 203 RCL 14 204 STO IND 17 205 RCL 18 206 RTN 207 END |
( 345 bytes / SIZE 021 )
| STACK | INPUTS | OUTPUTS |
| Y | h | / |
| X | n | Laplacian(f) |
where n is the number of variables
Example: f(x) = exp(-x2.t).Ln(y2+z) x = 1 , y = 1 , z = 1 , t = 1
LBL "T" RCL 04
X^2 *
"T" ASTO 00
RCL 01 *
RCL 03 RTN
1 STO 01 STO 02 STO 03 STO 04
X^2
E^X +
CHS
RCL 02 LN
-With h = 0.1
0.1 ENTER^
4 XEQ "LAPL" >>>>
Laplacian(f) = ¶2f / ¶x2
+ ¶2f / ¶y2
+ ¶2f / ¶z2
+ ¶2f / ¶t2
= 0.673013891
---Execution time = 54s---
( the last 3 decimals should be 932 )
-With h < 0 we get more quickly:
-0.01 ENTER^
4
R/S >>>> Lapl(f)
= 0.67300675
---Execution time = 25s---
First Derivatives of a Function of N variables
-If h > 0 "FD" employs a formula of order 10
-If h < 0 "FD" employs a formula of order 4 (
faster ) i-e
df/dx = (1/12h).[ f(x-2h) - 8.f(x-h)
+ 8.f(x+h) - f(x+2h) ]
>>> In both cases if the function returns an alpha data, 0 is returned
without calculating the whole formula if CF 03
Data Registers: • R00 = Function name ( Registers R00 thru Rnn are to be initialized before executing "FD" )
• R01 = x1 , • R02 = x2 , .......... , • Rnn = xn
R10 = h R11: temp
R13 = i R15 = xi
R12 = ¶f / ¶xi
R14 = 0 , h , 2h , .... R16: temp
Flag: F03
Subroutine: A program which computes
f(x1, ..... ,xn) assuming x1
is in R01 , ............ , xn is in Rnn
| STACK | INPUTS | OUTPUTS |
| Y | h | / |
| X | i | f 'xi = ¶f / ¶xi |
Example: f(x) = exp(-x2).Ln(y2+z)
x = 1 , y = 1 , z = 1
LBL "T" E^X
+
"T" ASTO 00
RCL 01 RCL 02
LN
1 STO 01 STO 02 STO 03
X^2
X^2
*
CHS
RCL 03 RTN
-With h = 0.1
0.1 ENTER^
1 XEQ "FD" >>>>
f 'x1 = ¶f / ¶x1
= -0.509989198 ( The last decimal should be a 5 )
---Execution time = 10s---
0.1 ENTER^
2 R/S
>>>> f 'x2 = ¶f
/ ¶x2 = 0.367879440
( The last decimal should be a 1 )
---Execution time = 11s---
0.1 ENTER^
3 R/S
>>>> f 'x3 = ¶f
/ ¶x3 = 0.183939720
( The last decimal should be a 1 )
---Execution time = 11s---
Notes:
-This program uses the same formula of order 10 as "dF" if
h > 0
-Though it's relatively fast, it may take a long time when this routine
is called many times as it is the case in the Riemanian section.
-So, if you choose h < 0 , "FD" will use a formula of order 4 , often less accurate but faster
-For instance, with h = -0.01
0.01 CHS ENTER^
1
R/S >>>>
f 'x1 = ¶f / ¶x1
= -0.509989196
---Execution time = 4.5s---
Second Derivatives of a Function of N variables
-If h > 0 "SD" employs a formula of order 10
-If h < 0 "SD" employs a formula of order 4 (
faster ) i-e
d2f/dx2 = (1/12h2).[ - f(x-2h) + 16.f(x-h) - 30.f(x) +16.f(x+h) - f(x+2h) ]
and, for the crossed derivative
f "xy = ¶2f / ¶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 )
>>> In both cases if the function places an alpha data in X-register,
0 is returned without calculating the whole formula
Data Registers: • R00 = Function name ( Registers R00 thru Rnn are to be initialized before executing "SD" )
• R01 = x1 , • R02 = x2 , .......... , • Rnn = xn
R10 = h R11: temp
R13 = xi
R15 = 0 , h , 2h , ....
R12 = ¶2f / ¶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: F04
CF 04 = The HP41 calculates gij , gij and Gkij
SF 04 = Only gij and gij are computed
Subroutines: The n(n+1)/2
functions gij with
i <= j
| 01 LBL "INIGC"
02 STO 09 03 STO 11 04 ENTER 05 X^2 06 + 07 2 08 / 09 STO 39 10 X<>Y 11 STO 10 12 CLX 13 40 14 + 15 STO 13 16 RCL 09 17 X^2 18 + 19 STO 26 20 RCL 09 21 E5 22 / 23 + 24 STO 14 25 CLX 26 STO 15 27 LBL 04 28 RCL 11 29 STO 12 30 LBL 05 |
31 RCL 11
32 RCL 12 33 Fµ 34 ASTO X 35 XEQ IND X 36 SIGN 37 X#0? 38 LASTX 39 STO IND 13 40 STO IND 14 41 STO IND 26 42 DSE 26 43 DSE 14 44 DSE 13 45 DSE 12 46 GTO 05 47 CLX 48 SIGN 49 ST+ 15 50 RCL 15 51 ST- 26 52 RCL 14 53 FRC 54 RCL 26 55 + 56 STO 14 57 DSE 11 58 GTO 04 59 RCL 09 60 X^2 |
61 ST+ X
62 RCL 39 63 + 64 40 65 + 66 .1 67 % 68 + 69 RCL 09 70 X^2 71 DSE X 72 - 73 CLRGX 74 RCL 09 75 1 76 + 77 E5 78 / 79 + 80 SIGN 81 LBL 06 82 STO IND L 83 ISG L 84 GTO 06 85 LASTX 86 FRC 87 41 88 + 89 RCL 39 90 + |
91 E-5
92 - 93 0 94 XROM "LS" 95 STO 40 96 1 97 STO 12 98 RCL 09 99 STO 11 100 ST* X 101 LASTX 102 INT 103 ST+ Y 104 LBL 07 105 RCL IND Y 106 STO IND Y 107 CLX 108 SIGN 109 ST+ Z 110 + 111 DSE Z 112 GTO 07 113 ISG 12 114 CLX 115 X<>Y 116 RCL 11 117 + 118 DSE X 119 RCL 12 120 X<> Z |
121 DSE 11
122 GTO 07 123 RCL 39 124 RCL 09 125 FS? 04 126 CLX 127 STO 23 128 2 129 + 130 * 131 RCL 13 132 + 133 STO 26 134 STO 29 135 FS? 04 136 GTO 14 137 LBL 10 138 RCL 09 139 STO 21 140 LBL 01 141 RCL 21 142 STO 22 143 LBL 02 144 CLX 145 STO 27 146 RCL 09 147 STO 24 148 LBL 03 149 RCL 23 150 RCL 24 |
151 ENTER
152 CLX 153 CIJK 154 X=0? 155 GTO 00 156 STO 28 157 CLX 158 STO 25 159 RCL 22 160 RCL 24 161 Fµ 162 ASTO 00 163 RCL 10 164 RCL 21 165 XROM "FD" 166 ST+ 25 167 RCL 21 168 RCL 24 169 Fµ 170 ASTO 00 171 RCL 10 172 RCL 22 173 XROM "FD" 174 ST+ 25 175 RCL 21 176 RCL 22 177 Fµ 178 ASTO 00 179 RCL 10 180 RCL 24 |
181 XROM "FD"
182 ST- 25 183 RCL 25 184 RCL 28 185 * 186 ST+ 27 187 LBL 00 188 DSE 24 189 GTO 03 190 RCL 27 191 2 192 / 193 STO IND 26 194 DSE 26 195 DSE 22 196 GTO 02 197 DSE 21 198 GTO 01 199 DSE 23 200 GTO 10 201 LBL 14 202 RCL 29 203 E3 204 / 205 41 206 + 207 END |
( 322 bytes / SIZE 41+n(n+1)(n+2)/2
)
| STACK | INPUTS | OUTPUTS |
| Y | h | / |
| X | n | bbb.eee (g&G) |
where n is the number of variables n < 8
Example: A 3D-Riemannian manifold, with a diagonal metric defined by
g11 = 1 + x2 y2z2
g22 = 1 + x2 + y2
+ z2
and gij = 0 if
i # j
g33 = 1 + x2 y2
+ x2 z2 + y2 z2
-Load for instance these 6 routines in main memory:
| 01 LBL "G11"
02 CLX 03 SIGN 04 RCL 01 05 RCL 02 06 RCL 03 07 * 08 * 09 X^2 10 + |
11 RTN
12 LBL "G22" 13 CLX 14 SIGN 15 RCL 01 16 X^2 17 + 18 RCL 02 19 X^2 20 + |
21 RCL 03
22 X^2 23 + 24 RTN 25 LBL "G33" 26 CLX 27 SIGN 28 RCL 01 29 RCL 02 30 * |
31 X^2
32 + 33 RCL 01 34 RCL 03 35 * 36 X^2 37 + 38 RCL 02 39 RCL 03 40 * |
41 X^2
42 + 43 RTN 44 LBL "G12" 45 LBL "G13" 46 LBL "G23" 47 ASTO X 48 RTN 49 END |
-At the point x = y = z = 1
1 STO 01 STO 02 STO 03
-We have n = 3 and if you choose h = 0.1
SF 00
0.1 ENTER^
3 XEQ "INIGC"
>>>> 41.070
---Execution time = 4m37s---
-And you find in registers R41 thru R52
R41 = g11 = 2
R42 = g12 = 0 R44 = g13
= 0
R47 = g11 = 1/2 R48 = g12
= 0 R50
= g13 = 0
R43 = g22 = 4 R45 = g23 =
0
R49 = g22 = 1/4
R51 = g23 = 0
R46 = g33 = 4
R52 = g33 = 1/4
-And in registers R53 to R70
R53 = G111
= 1/2 R54 = G112
= 1/2 R56 = G113
= 1/2
R55 = G122 = -1/2
R57 = G123 =
0
R58 = G133 = -1
R59 = G211
= -1/4 R60 = G212
= 1/4 R62 = G213
= 0
R61 = G222 = 1/4
R63 = G223 = 1/4
R64 = G233 = -1/2
R65 = G311
= -1/4 R66 = G312
= 0 R68 = G313
= 1/2
R67 = G322 = -1/4
R69 = G323 = 1/2
R70 = G333 = 1/2
-We have also R39 = 6 i-e n(n+1)/2 and R40 = 32 = g = det gij
Notes:
-Since all the gij functions are polynomials of degree <
5, we could get the same results with derivative-formulas of order 4
-So, if you choose h = -0.01, you get the same values in 2m57s
-Execute a SIZE at least
| n | SIZE |
| 2 | 053 |
| 3 | 071 |
| 4 | 101 |
| 5 | 146 |
| 6 | 209 |
| 7 | 293 |
-So, there are not enough data registers for n = 8 if CF 04.
b) Recalling Gij & Cijk
-CIJK is an M-Code routine that recalls Gkij
provided R39 = n(n+1)/2 and R09 = n
-These 2 registers have been initialized by "INIGC"
-CIJK checks that the data registers do exist and don't contain alpha
data
-However, it does not check for alpha data in the stack.
-The formula to get the address of the register Rmm that contains Gkij is ( with i < j )
m = 40 + i + ( j2 - j ) / 2 + ( k + 1
) n ( n + 1 ) / 2
08B "K"
00A "J"
009 "I"
003 "C"
378 C=c
03C RCR 3
106 A=C S&X
130 LDI S&X
027 027h=39d
146 A=A+C S&X
130 LDI S&X
200 200h
306 ?A<C S&X
381 GOTO
00A NEXIST
0A6 A<>C S&X
270 RAMSLCT
106 A=C S&X
038 READATA
0EE B<>C ALL
130 LDI S&X
01E 1Eh=30d
246 C=A-C S&X
270 RAMSLCT
038 READATA
10E A=C ALL
046 C=0 S&X
270 RAMSLCT
0EE B<>C ALL
355 chk A&C
050 and SETDEC
268 Q=C
0AE A<>C ALL
128 L=C
10E A=C ALL
0F8 C=X
135 C=
060 A*C
138 C=L
025 C=
060 AB*C
128 L=C
078 C=Z
10E A=C ALL
0B8 C=Y
30E ?A<C ALL
027 JC+04
068 Z=C
0AE A<>C ALL
0A8 Y=C
10E A=C ALL
135 C=
060 A*C
0B8 C=Y
2BE C=-C
025 C=
061 AB+C
04E C
35C =
090 2
269 C=
060 AB/C
138 C=L
025 C=
060 AB+C
078 C=Z
025 C=
060 AB+C
028 T=C
260 SETHEX
38D C
008 ->S&X
106 A=C S&X
130 LDI S&X
028 028h=040d
146 A=A+C S&X
378 C=c
03C RCR 3
146 A=A+C S&X
130 LDI S&X
200 200h
306 ?A<C S&X
381 GOTO
00A NEXIST
0A6 A<>C S&X
270 RAMSLCT
038 READATA
10E A=C ALL
046 C=0 S&X
270 RAMSLCT
0AE A<>C ALL
0E8 X=C
3E0 RTN
( 93 words )
| STACK | INPUTS | OUTPUTS |
| T | / | address-40 |
| Z | i | inf(i,j) |
| Y | j | sup(i,j) |
| X | k | Gkij |
Example1:
-With the metric above and after executing "INIGC"
3 ENTER^
ENTER^
2 XEQ "CIJK" >>>> G233
= -1/2
-Register T also contains 24 , so the data register which contains G233 is R64
Example2: CIJK may also be used to recall gij if k = -1 and gij if k = 0
-For instance, 3 ENTER^ ENTER^ 0 gives g33 = 1/4
and 2 ENTER^ ENTER^ 1 CHS yields g22 = 4
Notes:
-These routines check that the data registers R09 & R39 do
exist and that they do contain numbers
-However, they do not check X-Y-Z-registers for alpha data.
-Since n is a positive integer smaller than 10, CIJK only take
into account the first digit of the mantissas
c) Curvature Tensor
-The curvature tensor Rijkl ( 4 times covariant ) may be computed by the formula:
Rijkl = (1/2) ( ¶jk gil + ¶il gjk - ¶jl gik - ¶kik gjl ) + gmp ( GmjkGpil - GmjlGpik )
-The following listing includes "RIJKL" "RIJ" and
"RIJK^L"
| 01 LBL "RIJKL"
02 LBL 14 03 STO 24 04 RDN 05 STO 23 06 RDN 07 STO 22 08 X<>Y 09 STO 21 10 FS? 49 11 OFF 12 RCL 09 13 STO 25 14 CLX 15 STO 28 16 LBL 11 17 RCL 09 18 STO 26 19 LBL 13 20 RCL 25 21 RCL 26 22 ENTER 23 SIGN 24 CHS 25 CIJK 26 X=0? 27 GTO 00 |
28 STO 29
29 RCL 22 30 RCL 23 31 RCL 25 32 CIJK 33 STO 27 34 RCL 21 35 RCL 24 36 RCL 26 37 CIJK 38 ST* 27 39 RCL 22 40 RCL 24 41 RCL 25 42 CIJK 43 STO 19 44 RCL 21 45 RCL 23 46 RCL 26 47 CIJK 48 ST* 19 49 RCL 27 50 RCL 19 51 - 52 RCL 29 53 * 54 ST+ 28 |
55 LBL 00
56 DSE 26 57 GTO 13 58 DSE 25 59 GTO 11 60 RCL 21 61 RCL 24 62 Fµ 63 ASTO 00 64 RCL 10 65 RCL 22 66 RCL 23 67 XROM "SD" 68 ST+ 25 69 RCL 22 70 RCL 23 71 Fµ 72 ASTO 00 73 RCL 10 74 RCL 21 75 RCL 24 76 XROM "SD" 77 ST+ 25 78 RCL 21 79 RCL 23 80 Fµ 81 ASTO 00 |
82 RCL 10
83 RCL 22 84 RCL 24 85 XROM "SD" 86 ST- 25 87 RCL 22 88 RCL 24 89 Fµ 90 ASTO 00 91 RCL 10 92 RCL 21 93 RCL 23 94 XROM "SD" 95 ST- 25 96 RCL 25 97 2 98 / 99 RCL 28 100 + 101 SIGN 102 RCL 21 103 RCL 22 104 RCL 23 105 RCL 24 106 X<> L 107 RTN 108 GTO 14 |
109 LBL "RIJ"
110 LBL 10 111 STO 23 112 X<>Y 113 STO 21 114 RCL 09 115 STO 22 116 CLX 117 STO 20 118 LBL 01 119 RCL 09 120 STO 24 121 LBL 02 122 RCL 22 123 RCL 24 124 ENTER 125 CLX 126 CIJK 127 X=0? 128 GTO 00 129 STO 30 130 RCL 21 131 RCL 22 132 RCL 23 133 RCL 24 134XROM "RIJKL" 135 RCL 30 |
136 *
137 ST+ 20 138 LBL 00 139 DSE 24 140 GTO 02 141 DSE 22 142 GTO 01 143 RCL 23 144 SIGN 145 RCL 21 146 RCL 20 147 RTN 148 GTO 10 149 LBL "RIJK^L" 150 LBL 12 151 STO 30 152 RDN 153 STO 24 154 RDN 155 STO 23 156 X<>Y 157 STO 22 158 RCL 09 159 STO 21 160 CLX 161 STO 20 162 LBL 03 |
163 RCL 21
164 RCL 30 165 ENTER 166 CLX 167 CIJK 168 X=0? 169 GTO 00 170 STO 31 171 RCL 21 172 RCL 22 173 RCL 23 174 RCL 24 175 XROM "RIJKL" 176 RCL 31 177 * 178 ST+ 20 179 LBL 00 180 DSE 21 181 GTO 03 182 RCL 30 183 SIGN 184 RCL 22 185 RCL 23 186 RCL 24 187 RCL 20 188 RTN 189 GTO 12 190 END |
( 348 bytes / SIZE 41+n(n+1)(n+2)/2
)
| STACK | INPUTS | OUTPUTS |
| T | i | i |
| Z | j | j |
| Y | k | k |
| X | l | Rijkl |
| L | / | l |
Examples: The same diagonal metric
and with x = y = z = 1
1 ENTER^
2 ENTER^
1 ENTER^
2 XEQ "RIJKL" >>>> R1212 = -3/4
---Execution time = 28s---
-Likewise
2 ENTER^
3 ENTER^
1 ENTER^
3 R/S
>>>> R2313 = 1/2
---Execution time = 38s--- ( or 20s if
h < 0 )
Notes:
-The curvature tensor has many components but it has many symmetries
-So there are in fact N = n2 ( n2 - 1 )
/ 12 independant components i-e
N = 1 if n = 2
N = 6 if n = 3
N =20 if n = 4
1-time contravariant 3-times covariant
curvature tensor
-The 1-time contravariant 3-times covariant curvature tensor is given by:
Rijkl = gim Rjmkl
-"RIJK^L" calculates these components
| STACK | INPUTS | OUTPUTS |
| T | i | i |
| Z | j | j |
| Y | k | k |
| X | l | Rlijk |
| L | / | l |
Examples: The same diagonal metric
and with x = y = z = 1
2 ENTER^
1 ENTER^
2 ENTER^
1 XEQ "RIJK^L" >>>> R1212
= 3/8
---Execution time = 30s---
-Likewise
3 ENTER^
1 ENTER^
3 ENTER^
2 R/S
>>>> R2313 = -1/8
---Execution time = 39s--- ( or 21s if
h < 0 )
d) Ricci Tensor
-If we contract 2 indices of the curvature tensor ( the 2nd and the 4th ), we get the Ricci tensor which is symmetric
Rij = gkm Rikjm
-We can also contract other indices, but it would yield -Rij
( or 0 )
| STACK | INPUTS | OUTPUTS |
| Y | i | i |
| X | j | Rij |
| L | / | j |
Examples: The same diagonal metric
and with x = y = z = 1
1 ENTER^
1 XEQ "RIJ" >>>> R11 = -5/16
---Execution time = 1m49s---
-Likewise
2 ENTER^
3 R/S
>>>> R23 = -3/8
---Execution time = 2m34s---
-All the components - with i <= j - are
R11 = -5/16 R22
= -13/16 R33 = -11/16
R12 = +1/8 R23
= -3/8
R13 = +5/16
e) Scalar Riemannian Curvature
-Contracting the Ricci tensor, we get the scalar curvature
R = gij Rij
-"RR" calculates and stores R in register R38
-The following listing contains also "EIJ" & "WIJKL" ( cf the 2
paragraphs below )
| 01 LBL "RR"
02 CLX 03 STO 38 04 RCL 09 05 STO 21 06 LBL 01 07 RCL 21 08 STO 23 09 LBL 02 10 RCL 21 11 RCL 23 12 ENTER 13 CLX 14 CIJK 15 X=0? 16 GTO 00 17 STO 31 18 RCL 21 19 RCL 23 20 X=Y? 21 GTO 02 22 RCL 31 23 ST+ 31 24 RDN 25 LBL 02 26 XROM "RIJ" 27 RCL 31 28 * |
29 ST+ 38
30 LBL 00 31 DSE 23 32 GTO 02 33 DSE 21 34 GTO 01 35 RCL 38 36 RTN 37 LBL "EIJ" 38 XEQ 12 39 LBL 10 40 STO 34 41 X<>Y 42 STO 33 43 X<>Y 44 1 45 CHS 46 CIJK 47 STO 31 48 RCL 33 49 RCL 34 50 XROM "RIJ" 51 RCL 38 52 RCL 31 53 * 54 2 55 / 56 - |
57 SIGN
58 RCL 33 59 RCL 38 60 RCL 34 61 X<> L 62 CLD 63 RTN 64 GTO 10 65 LBL "WIJKL" 66 XEQ 12 67 LBL 13 68 STO 34 69 RDN 70 STO 33 71 RDN 72 STO 32 73 X<>Y 74 STO 31 75 CLX 76 STO 08 77 RCL 32 78 RCL 33 79 1 80 CHS 81 CIJK 82 STO 35 83 X=0? 84 GTO 00 |
85 RCL 31
86 RCL 34 87 XROM "RIJ" 88 RCL 35 89 * 90 ST+ 08 91 LBL 00 92 RCL 31 93 RCL 34 94 1 95 CHS 96 CIJK 97 ST* 35 98 X=0? 99 GTO 00 100 STO 36 101 RCL 32 102 RCL 33 103 XROM "RIJ" 104 RCL 36 105 * 106 ST+ 08 107 LBL 00 108 RCL 32 109 RCL 34 110 1 111 CHS 112 CIJK |
113 STO 37
114 X=0? 115 GTO 00 116 RCL 31 117 RCL 33 118 XROM "RIJ" 119 RCL 37 120 * 121 ST- 08 122 LBL 00 123 RCL 31 124 RCL 33 125 1 126 CHS 127 CIJK 128 ST* 37 129 X=0? 130 GTO 00 131 STO 36 132 RCL 32 133 RCL 34 134 XROM "RIJ" 135 RCL 36 136 * 137 ST- 08 138 LBL 00 139 RCL 37 140 RCL 35 |
141 -
142 RCL 38 143 * 144 RCL 09 145 1 146 - 147 / 148 RCL 08 149 + 150 RCL 09 151 2 152 - 153 / 154 STO 08 155 RCL 31 156 RCL 32 157 RCL 33 158 RCL 34 159 XROM "RIJKL" 160 ST+ 08 161 X<> 08 162 CLD 163 RTN 164 GTO 13 165 LBL 12 166 "R38=R?" 167 AVIEW 168 END |
(
293 bytes / SIZE 41+n(n+1)(n+2)/2 )
| STACK | INPUT | OUTPUT |
| X | / | R |
Examples: The same diagonal metric
and with x = y = z = 1
XEQ "RR" >>>> R = -17/32 = R38
---Execution time = 5m44s---
Notes:
-With a 2D-surface, the scalar curvature is twice the Gaussian curvature given by "KGS"
-"RR" uses 5 subroutine-levels, provided the "GIJ" do not call any subroutine.
-So, this program cannot be called as more than a first level subroutine.
-Always execute "RR" before computing Einstein tensor or Weyl tensor.
-Register R38 must contain R
f) Einstein Tensor
-Einstein tensor is a symmetric tensor defined as
Eij = Rij - (1/2) R gij
-Only after executing "RR" - which stores R into R38 -
execute EIJ to get the components of Einstein tensor ( see listing in paragraph
above )
| STACK | INPUTS | OUTPUTS |
| Y | i | i |
| X | j | Eij |
| L | / | j |
Examples: The same diagonal metric
and with x = y = z = 1
1 ENTER^
1 XEQ "EIJ" >>>> "R38=R?" and
finally E11 = 7/32
---Execution time = 1m49s---
-Likewise
2 ENTER^
3 R/S
>>>> E23 = -3/8
---Execution time = 2m35s---
-All the components - with i <= j - are
E11 = +7/32 E22
= +1/4 E33 = +3/8
E12 = +1/8
E23 = -3/8
E13 = +5/16
Notes:
-When you XEQ "EIJ" the HP41 displays "R38=R?" to remind that R38 must
contain the scalar curvature R before executing "EIJ"
-This message is not displayed again if you press R/S to get the other
components.
-"EIJ" uses 5 subroutine-levels, provided the "GIJ" do not call any
subroutine.
-So, this program cannot be called as more than a first level subroutine.
-In our example, the metric is diagonal. So, if i#j , Eij = Rij
-The cosmolical constant L is omitted here: the complete tensor is Eij = Rij - (1/2) R gij + L gij
-Elie Cartan proved that this tensor is the unique tensor T of order
2 involving derivatives of order < 3
and linear with respect to the 2nd derivatives that satisfies
the conservative equations
Di Tij = 0 for all j where Di is the covariant differentiation.
-They are the equations of General Relativity !
g) Weyl Tensor
-Weyl tensor is the covariant tensor of order 4 defined as
Wijkl = Rijkl + [ 1/(n-2) ] ( Ril gjk - Rik gjl + Rjk gil - Rjl gik ) + [ R/(n-1)/(n-2) ] ( gik gjl - gil gjk )
-Like Einstein tensor, execute "RR" first before executing "WIJKL"
so that data register R38 = R ( cf listing 2 paragraphs
above )
| STACK | INPUTS | OUTPUTS |
| T | i | i |
| Z | j | j |
| Y | k | k |
| X | l | Wijkl |
| L | / | l |
Examples: The same diagonal metric
and with x = y = z = 1
1 ENTER^
2 ENTER^
1 ENTER^
2 XEQ "WIJKL" >>>> "R38=R?" and finally
W1212 = 0
---Execution time = 4m11s---
-Likewise
2 ENTER^
3 ENTER^
1 ENTER^
3 R/S
>>>> W2313 = 0
---Execution time = 3m08s---
-In fact, all the components of Weyl tensor = 0 if n <
4
-The number of independent components is n(n-3)(n+1)(n+2)/12
Notes:
-When you XEQ "WIJKL" the HP41 displays "R38=R?" to remind that R38
must contain the scalar curvature R before executing "WIJKL"
-This message is not displayed again if you press R/S to get another
component.
-When the program stops, register R08 contains Rijkl
-"WIJKL" uses 5 subroutine-levels, provided the "GIJ" do not call any
subroutine.
-So, this program cannot be called as more than a first level subroutine.
-With the same expression for gii if i = 3
and g44 = 1 + x2 y2 z2
t2 at x = y = z = t =1 (
n = 4 )
you will get for example - after executing "INIGC" -
R = -45/32 and W1414
= 0 W2424 = +3/16 W3434
= -3/16
h) Covariant Derivative of a Tensor Field
of order 2 or 1
-The usual partial derivative of a vector or a tensor is not a tensor,
but the covariant derivatives are.
-They involve the partial derivatives and the Christoffel symbols.
-The formulae to get a tensor are not the same if the vector or the
tensor is covariant or contravariant:
Covariant Vector - CF 01 CF 02
Dk Vi = ¶k Vi - Gmik Vm
Contravariant Vector - SF 01 CF 02
Dk Vi = ¶k Vi + Gikm Vm
Covariant Tensor of order 2 - CF 01 CF 02
Dk Tij = ¶k Tij - Gmik Tmj - Gmkj Tim
Contravariant Tensor of order 2 - SF 02
Dk Tij = ¶k Tij + Gikm Tmj + Gjkm Tim
Mixed Tensor of order 2 - SF 01 CF 02
Dk Tij = ¶k Tij + Gikm Tmj - Gmkj Tim
-"DKTIJ" will calculate the components of all these tensors, but not
those of order 3 , 4 , ...
Data Registers: R00: temp
R01 thru Rnn coordinates of the point ( n < 8 )
R09 = n R11 to R18 are used by "FD"
R39 = n(n+1)/2 R41 thru R40+n(n+)(n+2)/2 = gij
, 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"
-This listing also contains other programs below
| 01 LBL "DKTIJ"
02 STO 23 03 RDN 04 X<>Y 05 FRC 06 CF 10 07 X#0? 08 SF 10 09 X<> L 10 X<>Y 11 FS? 10 12 1 13 STO 22 14 RDN 15 STO 21 16 X<>Y 17 STO 19 18 CLX 19 STO 24 20 RCL 09 21 STO 25 22 FS? 10 23 GTO 10 24 LBL 09 25 RCL 23 26 RCL 22 27 RCL 25 28 FS? 02 29 X<>Y 30 CIJK 31 X=0? 32 GTO 09 33 STO 20 34 RCL 21 35 RCL 25 36 RCL 09 37 ST* Y 38 - 39 + 40 ENTER 41 SIGN 42 - 43 RCL 19 44 + 45 RCL IND X 46 XEQ IND X 47 RCL 20 48 * 49 FC? 02 50 CHS |
51 ST+ 24
52 LBL 09 53 DSE 25 54 GTO 09 55 RCL 09 56 STO 25 57 LBL 10 58 RCL 23 59 RCL 21 60 RCL 25 61 FC? 01 62 FS? 02 63 X<>Y 64 CIJK 65 X=0? 66 GTO 10 67 STO 20 68 RCL 25 69 RCL 22 70 RCL 09 71 ST* Y 72 - 73 + 74 RCL 19 75 + 76 ENTER 77 SIGN 78 - 79 RCL IND X 80 XEQ IND X 81 RCL 20 82 * 83 FC? 01 84 FS? 02 85 CHS 86 ST- 24 87 LBL 10 88 DSE 25 89 GTO 10 90 RCL 21 91 RCL 22 92 RCL 09 93 ST* Y 94 - 95 + 96 RCL 19 97 + 98 1 99 - 100 RCL IND X |
101 STO 00
102 RCL 10 103 RCL 23 104 XROM "FD" 105 ST+ 24 106 RCL 23 107 SIGN 108 RCL 19 109 RCL 21 110 FC? 10 111 RCL 22 112 RCL 24 113 RTN 114 LBL "RCURL" 115 XEQ 14 116 ST- 26 117 LBL 06 118 RCL 21 119 STO 22 120 ISG 22 121 LBL 07 122 RCL 26 123 RCL 22 124 + 125 RCL IND X 126 STO 00 127 RCL 10 128 RCL 21 129 XROM "FD" 130 STO 25 131 RCL 26 132 RCL 21 133 + 134 RCL IND X 135 STO 00 136 RCL 10 137 RCL 22 138 XROM "FD" 139 RCL 25 140 - 141 STO IND 23 142 CLX 143 SIGN 144 ST+ 22 145 ST+ 23 146 RCL 09 147 RCL 22 148 X<=Y? 149 GTO 07 150 SIGN |
151 ST+ 21
152 CLX 153 RCL 21 154 X<Y? 155 GTO 06 156 LBL 13 157 RCL 23 158 DSE X 159 E3 160 / 161 RCL 24 162 + 163 RCL 19 164 SIGN 165 RDN 166 RTN 167 LBL "RDIV" 168 STO 19 169 INT 170 RCL 09 171 STO 21 172 + 173 STO 25 174 CLX 175 STO 24 176 DSE 25 177 LBL 04 178 RCL 09 179 STO 22 180 CLX 181 STO 23 182 LBL 05 183 RCL 21 184 RCL 22 185 ENTER 186 CIJK 187 ST+ 23 188 DSE 22 189 GTO 05 190 RCL IND 25 191 STO 00 192 XEQ IND 00 193 ST* 23 194 RCL 10 195 RCL 21 196 XROM "FD" 197 RCL 23 198 + 199 ST+ 24 200 DSE 25 |
201 DSE 21
202 GTO 04 203 RCL 19 204 SIGN 205 RCL 24 206 RTN 207 LBL"RLAPGR" 208 CLX 209 STO 24 210 RCL 39 211 RCL 09 212 STO 22 213 STO 23 214 2 215 + 216 * 217 41 218 + 219 STO 26 220 ST+ 23 221 DSE 23 222 LBL 08 223 RCL 10 224 RCL 22 225 XROM "FD" 226 STO IND 23 227 DSE 23 228 DSE 22 229 GTO 08 230 RCL 09 231 STO 21 232 ST+ 23 233 LBL 01 234 RCL 21 235 STO 22 236 LBL 02 237 CLX 238 STO 20 239 RCL 09 240 STO 25 241 LBL 03 242 RCL 21 243 RCL 22 244 RCL 25 245 CIJK 246 X=0? 247 GTO 00 248 RCL IND 23 249 * 250 ST+ 20 |
251 LBL 00
252 DSE 23 253 DSE 25 254 GTO 03 255 RCL 09 256 ST+ 23 257 RCL 10 258 RCL 21 259 RCL 22 260 XROM "SD" 261 ST- 20 262 RCL 21 263 RCL 22 264 X=Y? 265 GTO 00 266 RCL 20 267 ST+ 20 268 RDN 269 LBL 00 270 ENTER 271 CLX 272 CIJK 273 RCL 20 274 * 275 ST- 24 276 DSE 22 277 GTO 02 278 DSE 21 279 GTO 01 280 RCL 23 281 E3 282 / 283 RCL 26 284 + 285 RCL 24 286 RTN 287 LBL "V<>V^" 288 LBL 16 289 XEQ 14 290 LBL 11 291 CLX 292 STO 14 293 SIGN 294 STO 12 295 LBL 12 296 RCL 21 297 RCL 12 298 1 299 CHS 300 FC? 01 |
301 CLX
302 CIJK 303 RCL IND 26 304 * 305 ST+ 14 306 CLX 307 SIGN 308 ST+ 26 309 ST+ 12 310 RCL 09 311 RCL 12 312 X<=Y? 313 GTO 12 314 RCL 14 315 STO IND 23 316 RCL 09 317 ST- 26 318 SIGN 319 ST+ 21 320 ST+ 23 321 LASTX 322 RCL 21 323 X<=Y? 324 GTO 11 325 XEQ 13 326 RTN 327 GTO 16 328 LBL 14 329 STO 19 330 INT 331 STO 26 332 RCL 39 333 RCL 09 334 ST+ Z 335 2 336 + 337 * 338 41 339 + 340 X<Y? 341 X<>Y 342 STO 23 343 STO 24 344 SIGN 345 STO 21 346 END |
( 589 bytes / SIZE
41+n(n+1)(n+2)/2 )
| STACK | T INPUTS | OUTPUTS | V INPUTS | OUTPUTS |
| T | bbb.eee | bbb.eee | / | / |
| Z | i | i | bbb.eee | bbb.eee |
| Y | j | j | i | i |
| X | k | Dk Tij or Dk Tij or Dk Tij | k | Dk Vi or Dk Vi |
| L | / | k | / | k |
Examples: With the same metric as above and
still at the same point ( x , y , z ) = ( 1 , 1 , 1 )
Covariant Tensor Field of order 2 - CF 01 CF 02
T11 = x2 + y z3
T12 = x2 y z3
T13 = x2 + y2 + z2
T21 = x2 y z
T22 = x y z T23
= x2 + y z
T31 = x + y z
T32 = x2 y2 z2
T33 = x + y + z
-Load the corresponding routines in main memory ( the names are arbitrary,
provided they have less than 7 characters - global LBLs )
| 01 LBL "T11"
02 RCL 01 03 X^2 04 RCL 03 05 ENTER^ 06 X^2 07 * 08 RCL 02 09 * 10 + 11 RTN |
12 LBL "T21"
13 RCL 01 14 X^2 15 RCL 02 16 * 17 RCL 03 18 * 19 RTN 20 LBL "T31" 21 RCL 01 22 RCL 02 |
23 RCL 03
24 * 25 + 26 RTN 27 LBL "T12" 28 RCL 01 29 X^2 30 RCL 02 31 * 32 RCL 03 33 ST* Y |
34 X^2
35 * 36 RTN 37 LBL "T22" 38 RCL 01 39 RCL 02 40 RCL 03 41 * 42 * 43 RTN 44 LBL "T32" |
45 RCL 01
46 RCL 02 47 RCL 03 48 * 49 * 50 X^2 51 RTN 52 LBL "T13" 53 RCL 01 54 X^2 55 RCL 02 |
56 X^2
57 RCL 03 58 X^2 59 + 60 + 61 RTN 62 LBL "T23" 63 RCL 01 64 X^2 65 RCL 02 66 RCL 03 |
67 *
68 + 69 RTN 70 LBL "T33" 71 RCL 01 72 RCL 02 73 RCL 03 74 + 75 + 76 RTN 77 END |
-Store the names of theses subroutines in 9 consecutive registers,
without disturbing R41 to R70 , or R39 , R09 , R10 , R21 to R25 which are
used by "DKTIJ"
-For instance, in R27 thru R35 -> control number = 27.035
T11 T12
T13
R27 R30 R33
T21 T22
T23 stored in R28
R31 R34
T31 T32
T33
R29 R32 R35
-If you want to get D3 T12
Covariant Tensor Field of order 2 - CF 01 CF 02
CF 01 CF 02 ( and still SF 00 since the metric is orthogonal )
27.035
ENTER^
1 ENTER^
2 ENTER^
3 XEQ "DKTIJ"
>>>> D3 T12 = 1/4
---Execution time = 16s---
Contravariant Tensor Field of order 2 - CF 01 SF 02
-Assuming we keep the same expressions but the tensor is Tij instead of Tij
CF 01 SF 02
27.035
ENTER^
1 ENTER^
2 ENTER^
3 XEQ "DKTIJ"
>>>> D3 T12 = 5/4
---Execution time = 16s---
Mixed Tensor Field of order 2 - SF 01 CF 02
-Assuming again the same expressions with the tensor Tij
SF 01 CF 02
27.035
ENTER^
1 ENTER^
2 ENTER^
3 XEQ "DKTIJ"
>>>> D3 T12 =
3/4
---Execution time = 16s---
Covariant Vector Field - CF 01 CF 02
-Let's keep the first 3 functions as the components of the covariant vector field
CF 01 CF 02
27.029
ENTER^
2 ENTER^
3 XEQ "DKTIJ"
>>>> D3 V2 = -1/4
---Execution time = 12s---
Contravariant Vector Field - SF 01 CF 02
-Let's take again the first 3 functions as the components of the contravariant vector field
SF 01 CF 02
27.029
ENTER^
2 ENTER^
3 XEQ "DKTIJ"
>>>> D3 V2 = +1/4
---Execution time = 12s---
i) Curl of a Covariant Vector Field,
Divergence of a Contravariant
Vector Field,
Laplacian & Gradient of a Scalar Field
-The partial derivatives of a vector is not a tensor and
Curlij V is defined as Di Vj
- Dj Vi
-But if we apply the formulae given in the paragraph above, the Christoffel
symbols vanish and finally it yields
Curlij V = ¶i Vj - ¶j Vi
-This tensor is obviously antisymmetric, so Curlii V = 0 and Curlji V = - Curlij V
-So, to save space, "RCURL" only calculates and stores in consecutive
registers Curlij V with i < j & i =
1 , 2 , ... , n-1 i-e n(n-1)/2 components
Data Registers:
R00: temp
• R01 thru Rnn coordinates of the point ( n < 8 )
• R09 = n R20 to R26 are used by "RCURL"
• R39 = n(n+1)/2 R41 thru R40+n(n+)(n+2)/2 =
gij , gij , Gkij
• R10 = h
Flags: /
Subroutine:
"FD"
| STACK | INPUTS | OUTPUTS |
| X | bbb.eee | bbb.eeeCURL |
| L | / | bbb.eee |
Example: Again the same metric at the same
point and the following covariant vector field:
X = x2 + y z3
Y = x2 y z2
Z = x + y z
| 01 LBL "X"
02 RCL 01 03 X^2 04 RCL 03 05 ENTER^ 06 X^2 07 * 08 RCL 02 09 * 10 + 11 RTN 12 LBL "Y" 13 RCL 01 14 X^2 15 RCL 02 16 * 17 RCL 03 18 X^2 19 * 20 RTN 21 LBL "Z" 22 RCL 01 23 RCL 02 24 RCL 03 25 * 26 + 27 RTN |
-If you store the function names in R31 R32 R33 -> control number =
31.033
31.033 XEQ "RCURL" >>>> 71.073 ---Execution time = 53s---
-And you get
R71 = -1
0 +1 -2
R72 = 2 so,
the complete result is
-1 0
-1
R73 = 1
2 1
0
Notes
-The HP41 stores the results in a block of consecutive registers starting
at RBB
where BB = Max ( bb+n , 41+n(n+1)(n+2)/2 )
-So, the data are not disturb.
-This definition will not give the same results as "CDGL" in cylindrical
and spherical coordinates.
because the usual formulae in 3D-Euclidean space involve what
is called "Vraies Grandeurs" in French ( "True Magnitudes" in English ?
)
which are neither the covariant nor the contravariant components,
except with an orthonormal basis
-In tensor calculus, this "vraie grandeur" is the product of the contravariant
component Vk by the modulus of the corresponding vector ek
of the basis
or the quotient of the covariant component Vk by
the same quantity:
Vk ek = Vk / ek
= (gkk)1/2 Vk = (gkk) -1/2
Vk ( No summation
here )
Divergence of a Contravariant Vector Field
-The divergence of a vector field V is given by
Div V = Di Vi = ¶i Vi + Giik Vk
-It is a scalar
Data Registers:
R00: temp
• R01 thru Rnn coordinates of the point ( n < 8 )
• R09 = n R19 to R25 are used
by "RDIV" • R39 = n(n+1)/2
R41 thru R40+n(n+)(n+2)/2 = gij , 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
Note:
-The HP41 stores the results in a block of consecutive registers starting
at RBB
where BB = Max ( bb+n , 41+n(n+1)(n+2)/2 )
8°) An Example of a non-diagonal metric Tensor
( 3D-manifolfd )
g11 = 1 + x2 y2
z2
g12 = x y
g22 = 1 + x2 + y2
+ z2
g13 = x z
g33 = 1 + x2 y2
+ x2 z2 + y2 z2
g23 = y z
-Load for instance the 6 routines in main memory:
| 01 LBL "G11"
02 CLX 03 SIGN 04 RCL 01 05 RCL 02 06 RCL 03 07 * 08 * 09 X^2 10 + 11 RTN 12 LBL "G22" |
13 CLX
14 SIGN 15 RCL 01 16 X^2 17 + 18 RCL 02 19 X^2 20 + 21 RCL 03 22 X^2 23 + 24 RTN |
25 LBL "G33"
26 CLX 27 SIGN 28 RCL 01 29 RCL 02 30 * 31 X^2 32 + 33 RCL 01 34 RCL 03 35 * 36 X^2 |
37 +
38 RCL 02 39 RCL 03 40 * 41 X^2 42 + 43 RTN 44 LBL"G12" 45 RCL 01 46 RCL 02 47 * 48 RTN |
49 LNL "G13"
50 RCL 01 51 RCL 03 52 * 53 RTN 54 LBL "G23" 55 RCL 02 56 RCL 03 57 * 58 RTN 59 END |
-At the point x = y = z = 1
Christoffel Symbols:
1 STO 01 STO 02 STO 03 ( SIZE 071 at least ) and if you choose h = 0.1
• SF 03
to avoid the test with alpha data - unuseful here since all the gijare
different from 0
• CF 04
to get the Gkij
0.1 ENTER^
3 XEQ "INIGC" >>>>
41.070
---Execution time = 21m24s---
-And you find in registers R41 thru R52
R41 = g11 = 2
R42 = g12 = 1 R44 = g13
= 1
R47 = g11 = +5/8 R48 = g12
= -1/8 R50
= g13 = -1/8
R43 = g22 = 4 R45 = g23 =
1
R49 = g22 = +7/24
R51 = g23 = -1/24
R46 = g33 = 4
R52 = g33 = +7/24
>>> R40 = g = det gij = 24
-And in registers R53 to R70
R53 = G111
= 5/8 R54 = G112
= 1/2 R56 = G113
= 3/8
R55 = G122 = -1/8
R57 = G123 = -3/8
R58 = G133 = -3/4
R59 = G211
= -1/8 R60 = G212
= 1/6 R62 = G213
= -5/24
R61 = G222 = 7/24
R63 = G223 = 5/24
R64 = G233 = -1/4
R65 = G311
= -1/8 R66 = G312
= -1/6 R68
= G313 = 11/24
R67 = G322 = -1/24
R69 = G323 = 13/24
R70 = G333 =
3/4
-We have also R39 = 6 i-e n(n+1)/2
Curvature Tensor:
2 ENTER^
3 ENTER^
1 ENTER^
3 XEQ "RIJKL" >>>> R2313 = -7/24
---Execution time = 86s---
2 ENTER^
1 ENTER^
2 ENTER^
1 XEQ "RIJK^L" >>>> R1212
= -13/96
---Execution time = 232s---
Ricci Tensor:
1 ENTER^
3 XEQ "RIJ" >>>> R13
= -17/96
---Execution time = 12m26s---
-All the components - with i <= j - are
R11 = +7/96
R22 = -79/288 R33 =
-5/144
R12 = -28/96
R23 = -277/288
R13 = -17/96
Scalar Curvature:
XEQ "RR" >>>> R = +11/72 = R38 ---Execution time = 59mn36s--- ( 40m55s with h < 0 )
>>> You can change the h-value in R10 to get faster results, for example, h = - 0.01 STO 10 )
Einstein Tensor: ( after executing "RR" )
2 ENTER^
3 XEQ "EIJ" >>>> "R38=R?"
and E23 = -299/288
---Execution time = 6m59s--- ( with h < 0 )
-All the components - with i <= j - are
E11 = -23/288
E22 = -167/288 E33
= -49/144
E12 = -53/144
E23 = -299/288
E13 = -73/288
Weyl Tensor: ( after executing "RR" )
1 ENTER^
2 ENTER^
1 ENTER^
2 XEQ "WIJKL" >>>> "R38=R?" and finally
W1212 = 0
-Not a surprise since W = 0 if n < 4.
Differential Operators:
-With the same fields as above ( FF X Y Z T11 ..... T33 ) you will get:
Lap(f) = 317/96 ---Execution time = 140s---
Curlij & Gradk are unchanged because they don't depend on the metric.
Div(X Y Z ) = 10.5 ---Execution time = 33s---
D3T12 = 31/24 ---Execution time = 18s---
-Likewise, D3T33 = 1 D1T23 = -5/12 D2T22 = 19/24 ... and so on ...
Notes:
-As you can see, though dim = 3 , several routines are very slow, especially
"RR" to get the scalar curvature.
-That's why a good emulator is recommended.
-Practically, there are often several gij = 0, even if the
metric is not diagonal and execution time could be much smaller.
-Remember also that the programs run faster if you choose h < 0.
9°) Schwarzchild Metric ( 4D-Spacetime )
-Expressed in spherical coordinates ( r , u , v , ct ), we have:
ds2 = ( 1 - a / r ) -1 dr2 + r2 ( du2 + Sin2 u dv2 ) - ( 1 - a / r ) c2 dt2 where a = 2GM/c2
-Let's store a in register R07 and load the following routines ( we
must add the 6 functions LBL "G12" ........... LBL "G34" )
| 01 LBL "G11"
02 CLX 03 SIGN 04 RCL 07 05 RCL 01 06 / 07 - 08 1/X |
09 RTN
10 LBL "G22" 11 RCL 01 12 X^2 13 RTN 14 LBL "G33" 15 RCL 01 16 RCL 02 |
17 SIN
18 * 19 X^2 20 RTN 21 LBL "G44" 22 RCL 07 23 RCL 01 24 / |
25 ENTER^
26 CLX 27 SIGN 28 - 29 RTN 30 LBL "G12" 31 LBL "G13" 32 LBL "G14" |
33 LBL "G23"
34 LBL "G24" 35 LBL "G34" 36 ASTO X 37 RTN 38 END |
-Let's take, for example r = 2 u = PI/6
v = 0 ct = 0 and a = 1
2 STO 01 PI 6 / STO 02 CLX STO 03 STO 04 1 STO 07
Initialization:
XEQ "RAD" CF 04 CF 03
-Dimension = n = 4 and if you choose h = 0.1
0.1 ENTER^
4 XEQ INIGC
>>> > 41.100
---Execution time = 8m18s---
And the Christoffel symbols are ( after rounding the results given by the HP41 )
R61 = G111
= -1/4 R62 = G112
= 0 R64 = G113
= 0 R67 = G144
= 0
R63 = G122 = -1
R65 = G123 =
0 R68 = G134
= 0
R66 = G133 = -1/4
R69 = G124 =
0
R70 = G114 = 1/16
R71 = G211
= 0 R72 = G212
= 1/2 R74 = G213
= 0
R77 = G244 = 0
R73 = G222 =
0 R75 = G223
= 0
R78 = G234 =
0
R76 = G233 = -sqrt(3)/4
R79 = G224 =
0
R80 = G214 =
0
R81 = G311
= 0 R82 = G312
= 0 R84 = G313
= 1/2 R87 = G314
= 0
R83 = G322 = 0
R85 = G323 = sqrt(3)
R88 = G324 = 0
R86 = G333 =
0 R89 = G334
= 0
R90 = G344 =
0
R91 = G411
= 0 R92 = G412
= 0 R94 = G413
= 0 R97 = G414
= 1/4
R93 = G422 =
0 R95 = G423
= 0 R98 = G424
= 0
R96 = G433 =
0 R99 = G434
= 0
R100 = G444 = 0
Curvature Tensor:
1 ENTER^
2 ENTER^
1 ENTER^
2 XEQ "RIJKL" >>>> R1212 = -1/2
( HP41 gives -0.500001321 )
---Execution time = 30s---
1 ENTER^
4 ENTER^
1 ENTER^
4 R/S >>>>
R1414 = -0.124999921 ( exact result
= -1/8 )
3 ENTER^
4 ENTER^
3 ENTER^
4 R/S >>>>
R3434 = 1/32
... and so on ...
Ricci Tensor:
2 ENTER^
3 XEQ "RIJ" >>>> R23
= 0
---Execution time = 3m14s---
-All the components are actually zero.
-However, roundoff-errors are unavoidable and the HP41 sometime returns
small numbers instead of 0, for example:
R11 = -0.000000820
Scalar Curvature:
XEQ "RR" >>>> R = -0.000000846 ( exact result = 0 )
Einstein Tensor:
2 ENTER^
3 XEQ "EIJ" >>>> E23
= 0
-Here again, all the components are 0
Weyl Tensor:
1 ENTER^
2 ENTER^
1 ENTER^
2 XEQ "WIJKL" >>>> "R38=R?" and finally
W1212 = -0.500000095 ( exact result = -1/2 )
2 ENTER^
3 ENTER^
2 ENTER^
3 R/S >>>>
W2323 = 0.500000097 ( exact result = +1/2 )
1 ENTER^
4 ENTER^
1 ENTER^
4 R/S >>>>
W1414 = -0.125000024 ( exact result = -1/8 )
3 ENTER^
4 ENTER^
3 ENTER^
4 R/S >>>>
W3434 = 0.031250006 ( exact result = +1/32 )
.... and so on ...
10°) Non-Riemannian Geometry
-In Riemannian manifolds, the metric tensor determines the connexion
by the Christoffel symbols.
-In more general differential manifolds, the linear connexions Gkij
are independant functions and may be non-symmetric.
>>> So, you have to load in main memory n3
functions called LBL "C111" LBL "C112" LBL "C121"
and so on...
-Registers R09 & R10 are also to be initialized before executing the programs below:
• R09 = n < 9
• R10 = h
-And of course,
• R01 = x1 ................. • Rnn = xn ( n < 9 )
>>> If the tensor depends on the metric tensor, execute "INIGC"
first - with SF 04 to avoid calculating the Christoffel symbols.
a) Curvature Tensor
(1-3)
-If CF 02, "2RIJK^L" calculates the curvature tensor defined as
Rlijk = ¶jGlik - ¶kGlij + GlmjGmik - GlmkGmij
-If SF 02 this routine calculates the curvature tensor defined by the transposed connexion Gkji
Rlijk = ¶jGlki - ¶kGlji + GljmGmki - GlkmGmji
-This listing also contains other programs below.
| 01 LBL "2RIJK^L"
02 LBL 10 03 STO 14 04 RDN 05 STO 13 06 RDN 07 STO 12 08 X<>Y 09 STO 11 10 CLX 11 STO 20 12 RCL 09 13 STO 15 14 LBL 01 15 RCL 15 16 RCL 13 17 FS? 02 18 X<>Y 19 RCL 14 20 -DIFF MANIF 21 ASTO X 22 XEQ IND X 23 SIGN 24 X#0? 25 LASTX 26 STO 26 27 RCL 11 28 RCL 12 29 FS? 02 30 X<>Y 31 RCL 15 32 -DIFF MANIF 33 ASTO X 34 XEQ IND X 35 SIGN 36 X#0? |
37 LASTX
38 ST* 26 39 RCL 15 40 RCL 12 41 FS? 02 42 X<>Y 43 RCL 14 44 -DIFF MANIF 45 ASTO X 46 XEQ IND X 47 SIGN 48 X#0? 49 LASTX 50 STO 19 51 RCL 11 52 RCL 13 53 FS? 02 54 X<>Y 55 RCL 15 56 -DIFF MANIF 57 ASTO X 58 XEQ IND X 59 SIGN 60 X#0? 61 LASTX 62 RCL 19 63 * 64 RCL 26 65 - 66 ST+ 20 67 DSE 15 68 GTO 01 69 RCL 11 70 STO 21 71 RCL 12 72 STO 22 |
73 FS? 02
74 X<>Y 75 RCL 14 76 STO 24 77 -DIFF MANIF 78 ASTO 00 79 RCL 10 80 RCL 13 81 STO 23 82 XROM "FD" 83 ST- 20 84 RCL 21 85 RCL 23 86 FS? 02 87 X<>Y 88 RCL 24 89 -DIFF MANIF 90 ASTO 00 91 RCL 10 92 RCL 22 93 XROM "FD" 94 ST+ 20 95 RCL 24 96 SIGN 97 RCL 21 98 RCL 22 99 RCL 23 100 RCL 20 101 RTN 102 GTO 10 103 LBL "RQIJ" 104 LBL 12 105 STO 23 106 X<>Y 107 STO 28 108 RCL 09 |
109 STO 24
110 CLX 111 STO 27 112 LBL 02 113 RCL 28 114 RCL 24 115 FS? 00 116 X<>Y 117 RCL 23 118 RCL 24 119 XROM "2RIJK^L" 120 ST+ 27 121 DSE 24 122 GTO 02 123 RCL 23 124 SIGN 125 RCL 28 126 RCL 27 127 RTN 128 GTO 12 129 LBL "RQ" 130 LBL 09 131 RCL 09 132 STO 28 133 CLX 134 STO 38 135 LBL 03 136 RCL 09 137 STO 29 138 LBL 04 139 RCL 28 140 RCL 29 141 ENTER 142 CLX 143 CIJK 144 X=0? |
145 GTO 00
146 STO 30 147 RCL 28 148 RCL 29 149 XROM "RQIJ" 150 RCL 30 151 * 152 ST+ 38 153 LBL 00 154 DSE 29 155 GTO 04 156 DSE 28 157 GTO 03 158 RCL 38 159 RTN 160 GTO 09 161 LBL "SIJK" 162 LBL 05 163 STO 13 164 RDN 165 STO 12 166 X<>Y 167 STO 11 168 R^ 169 -DIFF MANIF 170 ASTO X 171 XEQ IND X 172 SIGN 173 X#0? 174 LASTX 175 STO 14 176 RCL 11 177 RCL 12 178 RCL 13 179 -DIFF MANIF 180 ASTO X |
181 XEQ IND X
182 SIGN 183 X#0? 184 LASTX 185 RCL 14 186 - 187 STO 14 188 RCL 13 189 SIGN 190 RCL 11 191 RCL 12 192 RCL 14 193 RTN 194 GTO 05 195 LBL "SJ" 196 LBL 06 197 STO 11 198 RCL 09 199 STO 12 200 CLX 201 STO 15 202 LBL 07 203 RCL 11 204 RCL 12 205 ENTER 206 XROM "SIJ" 207 ST+ 15 208 DSE 12 209 GTO 07 210 RCL 11 211 SIGN 212 RCL 15 213 RTN 214 GTO 06 215 END |
( 349 bytes / SIZE var. )
| STACK | INPUTS | OUTPUTS |
| T | i | i |
| Z | j | j |
| Y | k | k |
| X | l | Rlijk |
| L | / | l |
Example: Let's take the 3D-manifold with
the same metric as above at the point x = y = z = 1
g11 = 1 + x2 y2
z2
g12 = 0
g22 = 1 + x2 + y2
+ z2
g13 = 0
g33 = 1 + x2 y2
+ x2 z2 + y2 z2
g23 = 0
-Load the 6 routines in main memory:
| 01 LBL "G11"
02 CLX 03 SIGN 04 RCL 01 05 RCL 02 06 RCL 03 07 * 08 * 09 X^2 10 + |
11 RTN
12 LBL "G22" 13 CLX 14 SIGN 15 RCL 01 16 X^2 17 + 18 RCL 02 19 X^2 20 + |
21 RCL 03
22 X^2 23 + 24 RTN 25 LBL "G33" 26 CLX 27 SIGN 28 RCL 01 29 RCL 02 30 * |
31 X^2
32 + 33 RCL 01 34 RCL 03 35 * 36 X^2 37 + 38 RCL 02 39 RCL 03 40 * |
41 X^2
42 + 43 RTN 44 LBL "G12" 45 LBL "G13" 46 LBL "G23" 47 ASTO X 48 RTN 49 END |
-But let's change the connexion as follows:
G111 =
x y2 z2 / g11
G122 = - x / g11
G112
= x2 y z2 / g11 = G121G123
= G132 = 0
G113
= x2 y2 z / g11 = G131G133
= - x ( y2 + z2 ) / g11
G211 =
- x2 y z2 / g22
G222 = y / g22
G212
= x / g22 = G221
G223 = G232
= z / g22
G213
= 0 = G231
G233 = - y ( x2
+ z2 ) / g22
G311 =
- x2 y2 z / g33
G322 = - z / g33
so, it is the same as the Riemannian connexion
G312
= 0 = G321
G323 = - G332
= y ( x2 + z2 ) / g33
except that G323
= - G332
G313
= x ( y2 + z2 ) / g33 = G331G333
= z ( x2 + y2 ) / g33
-We have to load these 27 functions in main memory:
| 01 LBL "C111"
02 XEQ "G11" 03 RCL 02 04 RCL 03 05 * 06 X^2 07 RCL 01 08 * 09 X<>Y 10 / 11 RTN 12 LBL "C121" 13 LBL "C211" 14 XEQ "G11" 15 RCL 01 16 RCL 03 17 * 18 X^2 19 RCL 02 20 * 21 X<>Y 22 / 23 RTN 24 LBL "C131" 25 LBL "C311" 26 XEQ "G11" 27 RCL 01 28 RCL 02 29 * 30 X^2 |
31 RCL 03
32 * 33 X<>Y 34 / 35 RTN 36 LBL "C221" 37 XEQ "G11" 38 RCL 01 39 CHS 40 X<>Y 41 / 42 RTN 43 LBL "C231" 44 LBL "C321" 45 CLX 46 RTN 47 LBL "C331" 48 XEQ "G11" 49 RCL 02 50 X^2 51 RCL 03 52 X^2 53 + 54 RCL 01 55 * 56 CHS 57 X<>Y 58 / 59 RTN 60 LBL "C112" |
61 XEQ "G22"
62 RCL 01 63 RCL 03 64 * 65 X^2 66 RCL 02 67 * 68 CHS 69 X<>Y 70 / 71 RTN 72 LBL "C122" 73 LBL "C212" 74 XEQ "G22" 75 RCL 01 76 X<>Y 77 / 78 RTN 79 LBL "C132" 80 LBL "C312" 81 CLX 82 RTN 83 LBL "C222" 84 XEQ "G22" 85 RCL 02 86 X<>Y 87 / 88 RTN 89 LBL "C232" 90 LBL "C322" |
91 XEQ "G22"
92 RCL 03 93 X<>Y 94 / 95 RTN 96 LBL "C332" 97 XEQ "G22" 98 RCL 01 99 X^2 100 RCL 03 101 X^2 102 + 103 RCL 02 104 * 105 CHS 106 X<>Y 107 / 108 RTN 109 LBL "C113" 110 XEQ "G33" 111 RCL 01 112 RCL 02 113 * 114 X^2 115 RCL 03 116 * 117 CHS 118 X<>Y 119 / 120 RTN |
121 LBL "C123"
122 LBL "C213" 123 CLX 124 RTN 125 LBL "C133" 126 LBL "C313" 127 XEQ "G33" 128 RCL 02 129 X^2 130 RCL 03 131 X^2 132 + 133 RCL 01 134 * 135 X<>Y 136 / 137 RTN 138 LBL "C223" 139 XEQ "G33" 140 RCL 03 141 CHS 142 X<>Y 143 / 144 RTN 145 LBL "C323" 146 SF 09 147 GTO 00 148 LBL "C233" 149 CF 09 150 LBL 00 |
151 XEQ "G33"
152 RCL 01 153 X^2 154 RCL 03 155 X^2 156 + 157 RCL 02 158 * 159 X<>Y 160 / 161 FC? 09 162 RTN 163 CHS 164 RTN 165 LBL "C333" 166 XEQ "G33" 167 RCL 01 168 X^2 169 RCL 02 170 X^2 171 + 172 RCL 03 173 * 174 X<>Y 175 / 176 RTN 177 END |
-Though this is not necessary to calculate Rlijk
let's execute "INIGC" with SF04 to calculate the metric tensor
( covariant & contravariant ) only
1 STO 01 STO 02 STO 03 CF 03 SF 04
0.1 ENTER^
3 XEQ "INIGC" >>>>
41.052
---Execution time = 33s---
-We get the same values as before:
R41 = g11 = 2
R42 = g12 = 0 R44 = g13
= 0
R47 = g11 = 1/2 R48 = g12
= 0 R50
= g13 = 0
R43 = g22 = 4 R45 = g23 =
0
R49 = g22 = 1/4
R51 = g23 = 0
R46 = g33 = 4
R52 = g33 = 1/4
-Then, with CF 02
3 ENTER^
2 ENTER^
3 ENTER^
XEQ "2RIJK^L" >>>> R3323
= 0
---Execution time = 62s---
-Likewise with SF 02 ( transposed connexion )
3 ENTER^
2 ENTER^
3 ENTER^
R/S
>>>> R3323 = +1/4
---Execution time = 59s---
b) Ricci Tensors & Segmental
Curvatures
-Ricci tensor is Rij = Rmimj
-There is also a segmental curvature tensor Qij = Rmmij which equals zero in Riemannian manifolds.
Qij = ¶iGmmj
- ¶jGmmi
is actually a curl.
If CF 00 "RQIJ" returns the components of
Rij
-If CF 02 with the connexion
Gkij
-If SF 02 with the connexion
Gkji
If SF 00 "RQIJ" returns the components of Qij
-If CF 02 with the connexion
Gkij
-If SF 02 with the connexion
Gkji
| STACK | INPUTS | OUTPUTS |
| Y | i | i |
| X | j | Rij or Qij |
| L | / | j |
Rij if CF 00 (
connexion Gkij if CF 02
, transpose connexion Gkji
if SF 02 )
Qij if SF 00
( connexion Gkij if CF
02 , transpose connexion Gkji
if SF 02 )
Examples:
CF 00 CF 02
2 ENTER^
3 XEQ "RQIJ" >>>>
R23 = - 3/8
---Execution time = 2m44s---
-Remark that R32 = - 7/8 Ricci tensor is no more symmetric
CF 00 SF 02
2 ENTER^
3 R/S
>>>> R23 = - 9/8
SF 00 CF 02
2 ENTER^
3 R/S
>>>> Q23 = 0
SF 00 SF 02
2 ENTER^
3 R/S
>>>> Q23 = 0
c) Scalar Curvatures
-The contracted tensors R = gij Rij and Q = gij Qij gives 2 scalars.
-To obtain correct results, "INIGC" must be executed before "RQ" to
store gij in the proper registers.
| STACK | INPUT | OUTPUT |
| X | / | R or Q |
R if CF 00 ( connexion Gkij
if CF 02 , transpose connexion Gkji
if SF 02 )
Q if SF 00 ( connexion Gkij
if CF 02 , transpose connexion Gkji
if SF 02 )
Examples:
CF 00 CF 02
2 ENTER^
3 XEQ "RQ" >>>> R
= - 17/32
---Execution time = 9m40s---
CF 00 SF 02
2 ENTER^
3 R/S
>>>> R = - 9/16
SF 00 CF 02
2 ENTER^
3 R/S
>>>> Q = 0
SF 00 SF 02
2 ENTER^
3 R/S
>>>> Q = 0
Notes:
-A more important modification of the connection would be necessary
to get Q # 0
-"RQ" uses 5 subroutine levels, so the functions "CIJK" must use less
than 2 subroutine levels...
d) Torsion Tensor
-Since the Christoffel symbols are symmetric in Riemannian manifolds,
Gkij
- Gkji = 0
-In the general case, we get a tensor of order 3 that is called the
torsion tensor
Skij = Gkij - Gkji ( some authors define Skij = (1/2) ( Gkij - Gkji ) others by Skij = ( Gkji - Gkij ) )
-Though Gkij is not
a tensor, Skij is a tensor !
| STACK | INPUTS | OUTPUTS |
| Z | i | i |
| Y | j | j |
| X | k | Skij |
| L | / | k |
Examples:
2 ENTER^
3 ENTER^
XEQ "SIJK" >>>>
S323 = +1
3 ENTER^
2 ENTER^
3 R/S
>>>> S332 = -1
e) Torsion Vector
-Contracting Skij we get a covariant
vector Sj = Skjk
| STACK | INPUT | OUTPUT |
| X | j | Sj |
| L | / | j |
Examples:
1 XEQ "SJ" >>>> S1 =
0
2 R/S
>>>> S2 = +1
3 R/S
>>>> S3 = 0
Note:
-I called this program "SJ" instead of "SI" to avoid a conflict with
the Sine Integral...
11°) Differential Geometry QRG
| Functions | DESCRIPTION | REMARKS |
| KHT | Curvature & Torsion 3D-parametric Curve X(t) Y(t) Z(t) h t -> Khi(t) X<>Y Tau(t) | / |
| KGS | Gaussian Curvature of a surface z = f(x,y) • R00 = f name h y x -> K(Gauss) | / |
| dF | 1st- 2nd- 3rd-derivative of a function of 1 var. • R00 = f name h x -> f'(x) RDN f"(x) RDN f'''(x) | / |
| dF2 | Derivatives of a function of 2 variables • R00 = f name h y x -> f 'x RDN f 'y RDN f ''xx RDN f ''yy | / |
| dF3 | Derivatives of a function of 3 variables • R00 = f name h z y x -> f 'x RDN f 'y RDN f 'z RDN Df | / |
| MDV | Mixed Derivative of a function of 2 var. • R00 = f name h y x -> f ''xy | / |
| BHRM | Biharmonic Operator of a function of 3 var. • R00 = f name h z y x -> D2f | / |
| CDGL | Curl-Diverg-Grad-Lapl
h z y
x -> Curl1 RDN Curl2 RDN Curl3
RDN Div E
E = [ X(x,y,z) Y(x,y,z) Z(x,y,z) ] L = 4.019 ( contr. numb. all results ) |
CF01 CF 02 = Rect Coord
SF 01 CF 02 = Cyl Coord CF 01 SF 02 = Sph. Coord. |
| FD | First Deriv. funct n var. ( n < 10 ) • R00 = f name • R01 = x1 , ..... , • Rnn = xn h i -> ¶f / ¶xi | SF03 -> X'=0 if X = alpha
h > 0 = order10 h < 0 = order 4 |
| SD | 2nd Deriv. funct n var. ( n < 10 ) • R00 = f name • R01 = x1 , ..... , • Rnn = xn h j i -> ¶2f / ¶xi¶xj | X'=0 if X = alpha
h > 0 = order10 h < 0 = order 4 |
| LAPL | Laplacian funct n var. ( n < 10 ) • R00 = f name • R01 = x1 , ..... , • Rnn = xn h n -> D f | X'=0 if X = alpha
h > 0 = order10 h < 0 = order 4 |
| LS | Solves linear system, p = small number for tiny elements, rr = nb rows bbb.eee,rr p -> det M | / |
| / | N-DIMENSIONAL RIEMANNIAN MANIFOLDS | / |
| INIGC | Initialization: • R01 = x1 , ..... , • Rnn = xn h n -> 41.eee = control number of gij , gij and Gkij | SF 04 -> gij , gij only |
| CIJK | Recalls Gkij assuming INIGC has been executed i j k -> Gkij | / |
| RIJKL | Covariant Curvature Tensor i
j k l
-> Rijkl = (1/2) ( ¶jk
gil + ¶il gjk
- ¶jl gik - ¶ik
gjl )
+ gmp ( GmjkGpil - GmjlGpik ) |
/ |
| RIJK^L | Curvature Tensor ( 1-3 ) i j k l -> Rlijk = glm Rmijk | / |
| RIJ | Ricci Tensor i j -> Rij = gkm Rikjm | / |
| RR | Scalar Riemannian Curvature no input -> R = gij Rij | / |
| EIJ | Einstein Tensor i j -> Eij = Rij - (1/2) R gij | XEQ "RR" first to get R38 = R |
| WIJKL | Weyl Tensor
i j k
l -> Wijkl = Rijkl + [ 1/(n-2) ] ( Ril
gjk - Rik gjl + Rjk gil
- Rjl gik )
+ [ R/(n-1)/(n-2) ] ( gik gjl - gil gjk ) |
XEQ "RR" first to get R38 = R |
| DKTIJ | Covariant Derivative of a Tensor bbb.eee
i j k
-> Dk Tij or Dk Tij or
Dk Tij
or a Vector bbb.eee i k -> Dk Vi or Dk Vi |
CF 01 CF 02 covariant
SF 01 1 contrav SF 02 2 contrav |
| RCURL | Curl of a Covariant Vector Field bbb.eee -> Curlij V = ¶i Vj - ¶j Vi | / |
| RDIV | Divergence of a Contrav. Vector Fied bbb.eee -> Div V = Di Vi = ¶i Vi + Giik Vk | / |
| RLAPGR | Laplacian & Gradient of a Scalar Field •
R00 = f name no input -> Lapl(f) = gjk ( ¶jk
f - Gijk ¶i
f )
X<>Y bbb.eee(Gradient) |
/ |
| V<>V^ | Exchanging Covariant & Contravariant Coord. of a Vector bbb.eee(Vi) -> bbb.eee(Vi) | SF 01 = inverse operation |
| / | GENERAL DIFFERENTIAL MANIFOLDS including NON-RIEMANNIAN MANIFOLDS | / |
| 2RIJK^L | 2 Curvature Tensors: with the connexion Gkij
& with the transpose connexion Gkji
i j k
l -> Rlijk
Rlijk = ¶jGlik - ¶kGlij + GlmjGmik - GlmkGmij OR Rlijk = ¶jGlki - ¶kGlji + GljmGmki - GlkmGmji |
CF 02 = first Tensor
SF 02 = second Tensor |
| RQIJ | Ricci Tensor(s)
i j -> Rij
= Rmimj
& Segmental Curvatue Tensor(s) i j -> Qij = Rmmij = ¶iGmmj - ¶jGmmi |
CF 00 = Ricci Tensor
SF 00 = Segmental Tensor CF 02 & SF 02 as above |
| RQ | Scalar Curvatures no input ->
R = gij Rij
no input -> Q = gij Qij |
CF 00 = R
SF 00 = Q CF 02 & SF 02 as above |
| SIJK | Torsion Tensor i j k -> Skij = Gkij - Gkji | / |
| SJ | Torsion Vector j -> Sj = Skjk | / |
References:
[1] Elie Cartan - "Leçons sur la géométrie
des espaces de Riemann" - Gauthier-Villars, Paris ( in French
)
[2] Denis-Papin & Kaufmann - "Cours de calcul Tensoriel Appliqué"
- Albin Michel ( in French )
[3] List
of Formulas in Riemannian Geometry
[4] Nikodem J. Poplawski - "Spacetime and fields"- Department
of Physics, Indiana University, Bloomington, IN 47405, USA
[5] Marie-Antoinette Tonnelat - Theorie unitaire affine du champ
physique. J. Phys. Radium, 1951, 12 (2),
pp.81-88. <10.1051/jphysrad:0195100120208100>.
<jpa-00234360>