Differential Geometry for the HP48/50

 Overview
 

-These programs use numerical differentiation to calculate:

   Several differential operators
   Usual problems in Euclidean differential geometry
   Tensors in Riemaniann geometry
   A few tensors in non-Riemaniann manifolds
 

Notes:

-Most of these routines may also be applied with complex arguments and complex manifolds ( Complex General Relativity )
-However, use (0,1) instead of i. Otherwise, meaningless results could be returned when the variable 'i' is used.

-Most of the formulas to estimate derivatives are of order 10 ( h = 0.1 is "often" a good choice )
-This produces a good precision for the 1st derivatives but a low precision for 6th derivatives.

-With formulas of order 10, the execution time may be prohibitive to compute all tensors..
-So, you can also choose 4th-order methods with"FD" and "SD" ( in the Riemaniann & Manifolds sub-directories ), to get faster - but less accurate - results.

-In this case, choose h < 0 to select 4th order formulae. ( h = - 0.01 is "often" a good choice )

-The functions f you need to derive must be programs that you place in the stack:

     •  With 1 variable, f must take x in level 1 and return f(x) in level 1
     •  With 2 variables, f takes x in level 2 & y in level 1 and returns f(x,y) in level 1
     •  With 3 variables, f takes x in level 3, y in level 2 & z in level 1 and returns f(x,y,z) in level 1

     •  It's different with n variables:  x₁ , ............. , x_n are stored in n variables 'X1' .............. 'XN'
        f must use the contents of  'X1' .............. 'XN'  and return  f( x₁ , ....... , x_n ) in level 1

-You can use local or global variables to program your functions, but avoid lower case variables because some of the programs also use them
  and it could yield meaningless results.

-The execution times have been measured with an HP48GX.
-Surprisingly, they are sometimes longer with an HP50G !
 

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,

            a_k b^k = a₁ b¹ + a₂ b² + a₃ b³

         g_ij aⁱ a^j = g₁₁ a¹ a¹ + g₁₂ a¹ a² + g₁₃ a¹ a³
                    + g₂₁ a² a¹ + g₂₂ a² a² + g₂₃ a² a³
                    + g₃₁ a³ a¹ + g₃₂ a³ a² + g₃₃ a³ a³

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

        g_ij aⁱ a^j = g₁₁ a¹ a¹ + g₂₂ a² a² + g₃₃ a³ a³  + 2 ( g₁₂ a¹ a² + g₁₃ a¹ a³ + g₂₃ a² a³ )
 

-The partial derivative     ¶_m g_ij   means   ¶g_ij / ¶x^m
-Likewise,          ¶_km g_ij   =  ¶g_ij / ¶x^k¶x^m
 
 
 

Function

Desciption

DIFFGEOM  d12F  d34F  d56F  KHT  KHTN  ROM     LNG     SRV     SKS     CRVL     CRVN     TOL     ROM     N  VAR2     MDV     MDV12     MDV21     MDV4     KGM  VAR3     CuRL     DVG     GRaD     LaPL     VLaP     CYL          CuRL          DVG          GRaD          LaPL          VLaP     SPH          CuRL          DVG          GRaD          LaPL          VLaP     MDV3     BHRM     THRM  NVAR     FD     SD     TD     CuRL     DVG     GRaD     HESS     LaPL     BHRM     KGM      h      N      X1      X2      X3      X4      X5      X6      X7      X8      X9  RIEMANNIAN      h      N      X1      X2      X3      X4      G11      G22      G33      G44      G12      G13      G14      G23      G24      G34      INIGC      GIJ      G^IJ      CIJK      GETCIJK      RIJKL      RIJK^L      RIJ      RSC      R      EIJ      WIJKL      CuRL      DVG      GRaD      LaPL      COV^      COV¯      FD      SD  MANIFOLDS      h      N      X1      X2      X3      X4      G11      G22      G33      G44      G12      G13      G14      G23      G24      G34      INIG      GIJ      G^IJ      RIJK^L      RIJ      QIJ      RRR      QQQ      SIJK      SI      C111      C121      ........      FD      SD

Directory  1st & 2nd Derivatives of a Function of 1 variable  4th & 5th Derivatives of a Function of 1 variable  5th & 6th Derivatives of a Function of 1 variable  Curvature & Torsion of a Curve ( 3 Dim )  Curvature & Torsion of a Curve ( n Dim )  Subdirectory = Romberg Method      Arc Length of a Curve defined by y = f(x)      Area of a Surface of Revolution      Area of a Surface defined by z = f(x,y)      Arc Length of a Parametric Curve ( 3-Dim )      Arc Length of a Parametric Curve ( N-Dim )      Tolerance      A Subroutine to use Romberg Method      N = 1 , 2 , 4 , 8 , ......  Subdirectory = Functions of 2 Variables      Mixed Derivative of a Function of 2 vatiables d2F/dxdy      d3F/dxdy2      d3F/dx2dy      d4F/dx2dy2      Gaussian Curvature & Mean Curvature of a Surface.  Subdirectory = Functions of 3 Variables      Curl of a 3D-Vector Field ( Rectangular Coordinates )      Divergence of a 3D-Vector Field      Gradient of a function of 3 variables      Laplacian of a function of 3 variables      3D-Vector Laplacian       Subdirectory = Cylindrical Coordinates          Curl of a 3D-Vector Field          Divergence of a 3D-Vector Field          Gradient of a function of 3 variables          Laplacian of a function of 3 variables         3D-Vector Laplacian       Subdirectory = Spherical Coordinates          Curl of a 3D-Vector Field          Divergence of a 3D-Vector Field          Gradient of a function of 3 variables          Laplacian of a function of 3 variables         3D-Vector Laplacian       d3F/dxdydz      Biharmonic Operator ( 3 Dim. )      Triharmonic Operator ( 3 Dim. )  Subdirectory = Functions of N Variables      First Derivatives of a Function of N variables      Second Derivatives of a Function of N variables      Third Derivatives of a Function of N variables      Curl of an N-Dim Vector Field ( antisymmetric tensor )      Divergence of an N-Dim Vector Field      Gradient of a Function of N variables      Hessian Matrix      Laplacian of a Function of N variables      Biharmonic Operator ( N Dim. )      Gaussian & Mean Curvatures of a Hypersurface.      h = a "small" number      N = Number of Variables      1st variable      2nd variable      3rd variable      you can add      4th variable      5th variable       other variables      6th variable      7th variable       'X10'  'X11'  ...  if need be.      8th variable      9th variable  Subdirectory = Riemannian Manifolds      h = a "small" number      N = Dimension of the Riemannian manifold      1st coordinate      2nd coordinate      3rd coordinate      4th coordinate      G11         Programs to compute      G22      G33         the covariant coordinates      G44      G12         of the metric tensor      G13      G14      G23      G24      G34      Initialization      Covariant Metric Tensor ( a symmetric Matrix )      Contravariant Metric Tensor ( Inverse of [ GIJ ] )      List of the Christoffel Symbols       Recalling Cijk       Curvature Tensor ( 0-4 )       Curvature Tensor ( 1-3 )      Ricci Tensor      Calculates the Scalar Riemannian Curvature      Scalar Riemannian Curvature      Einstein Tensor      Weyl Tensor      Curl of a Covariant Vector-Field      Divergence of a Contravariant Vector-Field      Gradient of a Scalar Field      Laplacian of a Scalar Field      Covariant Derivative of a Contravariant Vector-Field      Covariant Derivative of a Covariant Vector-Field      First Derivatives of a Function of N variables      Second Derivatives of a Function of N variables  Subdirectory = General Manifolds      h = a "small" number      N = Dimension of the Manifold      1st coordinate      2nd coordinate      3rd coordinate      4th coordinate      G11         Programs to compute      G22      G33         the covariant coordinates      G44      G12         of the metric tensor      G13      G14      G23      G24      G34      Initialization      Covariant Metric Tensor ( a symmetric Matrix )      Contravariant Metric Tensor ( Inverse of [ GIJ ] )      Curvature Tensor ( 1-3 )      Ricci Tensor      Segmental Curvature Tensor      Scalar Curvature      Segmental Curvature      Torsion Tensor      Torsion Vector      Affine Connexions ( Programs to compute the Cijk )      .............................      .............................      First Derivatives of a Function of N variables      Second Derivatives of a Function of N variables

 
 
 

1st & 2nd Derivatives of a Function of 1 variable
 

-"d12F"  calculates the 1st & 2nd derivatives of a function of 1 variable f(x)
 

Formulae:

       df/dx = (1/2520.h).[ 2100.( f₁ - f_-1 ) - 600.( f₂ - f_-2 ) + 150.( f₃ - f_-3 ) - 25.( f₄ - f_-4 ) + 2.( f₅ - f_-5 ) ] + O(h¹⁰)
     d²f/dx² = (1/25200.h²).[ -73766 f₀ + 42000.( f₁ + f_-1 ) - 6000.( f₂ + f_-2 ) + 1000.( f₃ + f_-3 ) - 125.( f₄ + f_-4 ) + 8.( f₅ + f_-5 ) ] + O(h¹⁰)

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

-They are exact for any polynomial of degree < 11
 
 
 

STACK	INPUTS	OUTPUTS
Level 3	<< f >>	/
Level 2	x	f ' (x)
Level 1	h	f '' (x)

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

-With  h = 0.1

   << EXP LASTARG LN + >>   ENTER
                       2                           ENTER                     f '(2) = 7.88905609897
                      0.1                           d12F    >>>>          f "(2) = 7.13905610484

Note:

-The program in level 3 must take x in Level 1 and return f(x) in level 1
 

3rd & 4th Derivatives of a Function of 1 variable
 

-"d34F"  calculates the 3rd & 4th derivatives of a function of 1 variable f(x)
 

Formulae:
 

-We have:

    d³f/dx³ ~ [ -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) ] / h³

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

   h⁴ f⁽⁴⁾(x) ~  ( 192654 F - 140196 G + 52428 H - 9738 I +1261 J - 82 K ) / 15120

       with   A = f(x+h) - f(x-h)  ,  B = f(x+2h) - f(x-2h)  ,  C = f(x+3h) - f(x-3h)  ,  D = f(x+4h) - f(x-4h)  , E = f(x+5h) - f(x-5h)

        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)
 

-They are exact for any polynomial of degree < 11
 
 

STACK	INPUTS	OUTPUTS
Level 3	<< f >>	/
Level 2	x	f (3) (x)
Level 1	h	f (4) (x)

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

-With  h = 0.1

   << EXP LASTARG LN + >>   ENTER
                       2                           ENTER                      f ⁽³⁾ (2) = 7.63905605159
                      0.1                           d34F    >>>>           f ⁽⁴⁾ (2) = 7.01405238095
 

5th & 6th Derivatives of a Function of 1 variable
 

-"d56F"  calculates the 5th & 6th derivatives of a function of 1 variable f(x)
 

Formulae:
 

-Let  A = f(x+h) - f(x-h)  ,  B = f(x+2h) - f(x-2h)  ,  C = f(x+3h) - f(x-3h)  ,  D = f(x+4h) - f(x-4h)  , E = f(x+5h) - f(x-5h)

        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)

-We have:

   h⁵ f⁽⁵⁾(x) ~  ( 1938 A - 1872 B + 783 C - 152 D + 13 E ) / 288
   h⁶ f⁽⁶⁾(x) ~  ( -233244 F + 184110 G - 88920 H + 24795 I -3610 J + 247 K ) / 4560
 

-They are exact for any polynomial of degree < 11
 
 
 

STACK	INPUTS	OUTPUTS
Level 3	<< f >>	/
Level 2	x	f (5) (x)
Level 1	h	f (6) (x)

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

-With  h = 0.1

   << EXP LASTARG LN + >>   ENTER
                       2                           ENTER                      f ⁽⁵⁾ (2) = 8.13909027778
                      0.1                           d56F    >>>>           f ⁽⁶⁾ (2) = 5.51559210526

Note:

-The 4th derivatives are already difficult to evaluate with a great precision, but here, the exact values are:

  f ⁽⁵⁾ (2) = 8.13905609893...
  f ⁽⁶⁾ (2) = 5.51405609893...

-Thus, hardly 4 digits are significant for the 6th derivative !
 

Curvature & Torsion of a Parametric Curve ( 3 Dimensions )
 

-We assume that the curve is parametrized by 3 functions  x(t) , y(t) , z(t)

-The curvature is calculated by   khi = | r' x r'' | / | r' |³                                      where  the vector  r = [ x(t) , y(t) , z(t) ]
                   and the torsion by   tau = ( r' x r'' ) .r''' / | r' x r'' |²                        x = cross-product  and  . = dot-product
 
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< x(t) >>	/
Level 4	<< y(t) >>	/
Level 3	<< z(t) >>	/
Level 2	t	khi(t)
Level 1	h	tau (t)

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

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

-If you choose  h = 0.1

    RAD

    << EXP  LASTARG  COS  *  >>  ENTER
    << EXP  LASTARG  SIN   *  >>   ENTER
                  <<  LN  >>                       ENTER
                         1.3                             ENTER               k:  0.194807823431
                         0.1                               KHT      >>>>   t:  0.123665238776

-The exact values are   khi (1.3) = 0.194807823  &  tau (1.3) = 0.123665529  ( rounded to 9 decimals )
 

Curvature & Torsion of a Parametric Curve ( N Dimensions )
 

-The method follows the Gram-Schmidt process ( cf references [1] & [2] )
 
 
 

STACK	INPUTS	OUTPUTS
Level n+3	<< x1(t) >>	/
Level ...	.................	/
Level 4	<< xn(t) >>	/
Level 3	t	/
Level 2	h	khi (t)
Level 1	n	tau (t)

 
Example:      The curve defined by

     x₁(t) = t⁴       x₃(t) = t³      at the point  t = 1
     x₂(t) = t²       x₄(t) = e^t

-If you choose  h = 0.1

    << 4 ^  >>   ENTER
    <<  SQ >>   ENTER
    << 3 ^  >>   ENTER
    << EXP >>  ENTER
            1           ENTER
           0.1         ENTER             k:  0.142277240599
            4           KHTN   >>>     t:   0.121469323539

Note:

-If n = 3 , the torsion given ny the Gram-Schmidt process is multiplied by SIGN ( DET ( r' , r'' , r''' ) )  to get the same sign as  KHT above.
 

Arc Length of a Curve defined by  y = f(x)
 

-The arc length of the curve  y = f(x)   ( a < x < b )  is given by   L = §_a^b  ( 1 + (y')² )^1/2 dx
  "LNG" doesn't use this formula and avoids the calculation of dy/dx . It simply applies Pythagoras' theorem.
 
 
 

STACK	INPUTS	OUTPUTS
Level 3	<< f >>	/
Level 2	a	/
Level 1	b	L

 
Example:    Calculate the arc length of the curve   y = ln x      1 < x < 3

-If you choose a tolerance of 10 ^-9  , store this value in the variable 'tol'  then:

     <<  LN  >>  ENTER
             1          ENTER
             3            LNG      >>>    2.30198753457

Note:

-The HP48 displays the successive approximations
 

Area of a Surface of Revolution
 

-The rotation of the curve  y = f(x)   ( a < x < b )  around x-axis generates a surface of revolution given by   S = 2.pi  §_a^b  y ( 1 + (y')² )^1/2 dx
  "SRV" avoids the calculation of dy/dx : the area of a truncated cone is used with Romberg method.
 
 
 

STACK	INPUTS	OUTPUTS
Level 3	<< f >>	/
Level 2	a	/
Level 1	b	S

 
Example:    Evaluate the area of the surface of revolution generated by the rotation of the curve   y = sin x  ( 0 < x < pi )  around the x-axis.

-If you choose the same tolerance of 10 ^-9  , store this value in the variable 'tol'  then:

     RAD

     << SIN >>  ENTER
             0          ENTER
            PI           SRV      >>>    14.4235994518

Note:

-The HP48 displays the successive approximations
 

Area of a Surface defined by z = f(x,y)
 

-"SKS" computes the area of a surface defined by:     z = f(x,y)       a < x < b  ,   c < y < d

-The result could be obtained by the double integral  §_a^b §_c^d   ( 1 + f_x² + f_y²  )^1/2 dx dy

   where  f_x = ¶f/¶x  and   f_y = ¶f/¶y  are the partial derivatives with respect to x and y respectively.

-But "SKS" avoids the calculus of the partial derivatives:
-The intervals [a,b] and [c,d] are divided into n parts, and the approximate area is the sum of the areas of triangles.
 
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< f >>	/
Level 4	a	/
Level 3	b	/
Level 2	c	/
Level 1	d	S

 
Example:   Evaluate the area of the surface defined by z = ( 25 - x² - y² )^1/2  ,  0 < x < 2 ,  0 < y < 3
 

-Let's choose 10 ^-6 for the tolerance ( store this positive value in 'tol' )

    <<  SQ  SWAP  SQ  +  25  SWAP  -  Ö   >>   ENTER                                 ( Ö  = SQRT )
                                         0                                   ENTER
                                         2                                   ENTER
                                         0                                   ENTER
                                         3                                     SKS        >>>   6.65439661003

-With  tol = 10 ^-8 , it yields  6.6543961212    ( more accurate, but slower )

Notes:

-The HP48 displays the successive approximations
-As usual with Romberg method, n is multiplied by 2 at each iteration but here,
  execution time is multiplied by 4: we approximate a double integral.
-You can stop the iterations: simply press any key. The program will stop at the next step.
 

Arc Length of a Parametric Curve ( 3-Dim )
 

 "CRVL" evaluates the integral  §_a^b  [ ( dX/dt )² + ( dY/dt )² + ( dZ/dt )² ] ^1/2  dt

-Romberg method is used with Pythagora's theorem.
 
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< x(t) >>	/
Level 4	<< y(t) >>	/
Level 3	<< z(t) >>	/
Level 2	a	/
Level 1	b	L

 
Example:     X(t) = t⁴ ,  Y(t) = t²  ,  Z(t) = exp t   ;  a = 0 ,  b = 1

-If you choose a tolerance of 10 ^-9  , store this value in the variable 'tol'  then:

    <<  4  ^  >>  ENTER
    <<   SQ  >>  ENTER
    <<  EXP >>  ENTER
             0          ENTER
             1           CRVL   >>>>   L = 2.34211645877

Note:

-The HP48 displays the successive approximations
 

Arc Length of a Parametric Curve ( N-Dim ) N < 10
 

 "CRVN" evaluates the integral  §_a^b  [ ( dX₁/dt )² + ( dX₂/dt )² + ............ + ( dX_N/dt )² ] ^1/2  dt
 
 
 
 

STACK	INPUTS	OUTPUTS
Level n+3	<< x1(t) >>	/
Level ...	.................	/
Level 4	<< xn(t) >>	/
Level 3	a	/
Level 2	b	/
Level 1	n	L

 
Example:     X₁(t) = t⁴ ,  X₂(t) = t²  ,  X₃(t) = t³  ,  X₄(t) = exp t   ;  a = 0 ,  b = 1

-If you choose a tolerance of 10 ^-9  , store this value in the variable 'tol'  then:

    <<  4  ^  >>  ENTER
    <<  SQ   >>  ENTER
    <<  3  ^  >>  ENTER
    <<  EXP >>  ENTER
             0          ENTER
             1          ENTER
             4           CRVN  >>>>      L = 2.57907923132

Notes:

-The HP48 displays the successive approximations
-This program employs DOLIST ,  SLIST ... so it will not work with an HP48S/SX
 

Romberg Method
 

-Suppose that a sequence  I_n  tends to I  as  n tends to infinity  and  that the  "errors"  I - I_n  are "nearly" proportional to 1/n²

  If you want to use Romberg method to estimate the limit I
  ROM must be called by a program with the following structure:

  <<  1  'n'  STO
         DO
         .......  calculates  I_n

         UNTIL   ROM
         END   OVER  2  +  DROPN
  >>

-Please EDIT for instance LNG or another program of this subdirectory to get examples.
 

Mixed Derivative of a Function of 2 variables
 

  MDV calculates  ¶²f / ¶x¶y   with the formula:
 

     f "_xy = ¶²f / ¶x¶y = (1/50400.h²).[ 73766 f₀₀ + 42000.( f₁₁ + f _-1-1 - f₁₀ - f _-10 - f₀₁ - f_0-1 ) + 6000.( - f ₂₂ - f _-2-2 + f₂₀ + f _-20 + f₀₂ + f_0-2 )
                                                         + 1000.( f₃₃ + f _-3-3 - f₃₀ - f _-30 - f₀₃ - f_0-3 ) + 125.( - f ₄₄ - f _-4-4 + f₄₀ + f _-40 + f₀₄ + f_0-4 )
                                                               + 8.( f₅₅ + f _-5-5 - f₅₀ - f _-50 - f₀₅ - f_0-5 ) ] + O(h¹⁰)

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

STACK	INPUTS	OUTPUTS
Level 4	<< f(x,y) >>	/
Level 3	x	/
Level 2	y	/
Level 1	h	f ''xy

   Where the program in level 4 takes x & y in level 2 & 1 respectively and returns f(x,y) in level 1

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

-If you choose h = 0.1

    <<  LN  SWAP  SQ  NEG  EXP  *  >>   ENTER
                                 1                                  ENTER
                                 2                                  ENTER
                                0.1                                 MDV      >>>>    f "_xy  = ¶²f / ¶x¶y = -0.367879442194

-The exact result is  - 1/e =  -0.367879441171...
 

Mixed Derivative d³F/dx²dy & d³F/dxdy²
 

 MD12 & MD21 employ a method of order 11 to estimate  f'''_xxy = ¶³f / ¶x²¶y  and  f'''_xyy = ¶³f / ¶y²¶x

    h³ f'''_xxy ~  SUM_k=1,2,...,5  A_k [ - 2 f_0k + 2 f_0-k  + ( f_kk - f_-k-k + f_-kk - f_k-k ) ]

 With    A₁ = 5/6 ,  A₂ = -5/84  ,  A₃ = +5/756  ,  A₄ = -5/8064  ,  A₅ = +1/31500
 
 
 

STACK	INPUTS	OUTPUTS
Level 4	<< f(x,y) >>	/
Level 3	x	/
Level 2	y	/
Level 1	h	f'''xxy or  f'''xyy

   Where the program in level 4 takes x & y in level 2 & 1 respectively and returns f(x,y) in level 1

Example:       f(x,y,z) = Ln ( 1 + x² y )    x = 1  ,  y = 2
 

-Again with h = 0.1

    << SWAP  SQ  *  1  +  LN  >>   ENTER
                              1                         ENTER
                              2                         ENTER
                             0.1                       MD21      >>>>   f'''_xxy = ¶³f / ¶x²¶y  =  - 0.370370410386

-The exact value is  -10/27 =  - 0.370370370.....

-With the same inputs,  MD12  returns   f'''_xyy = ¶³f / ¶x¶y²  =  - 0.148148108438

-The exact result is  -4/27 =  - 0.148148148...
 

.Mixed Derivative d⁴F/dx²dy²
 

 "MDV4" evaluates   f'''_xxyy = ¶⁴f / ¶x²¶y²  with a formula of order 12

  h⁴ f'''_xxyy ~  20881861 f₀₀ / 3240000  + SUM_k=1,....,5  A_k [ f_kk + f_-k-k + f_k-k + f_-kk - 2 ( f_k0 + f_-k0 + f_0k + f_0-k ) ]

  where   A₁ = 5/3 ,  A₂ = -5/84  ,  A₃ = +5/1134  ,  A₄ = -5/16128  ,  A₅ = +1/78750
 
 
 

STACK	INPUTS	OUTPUTS
Level 4	<< f(x,y) >>	/
Level 3	x	/
Level 2	y	/
Level 1	h	f'''xxyy

   Where the program in level 4 takes x & y in level 2 & 1 respectively and returns f(x,y) in level 1

Example:    f(x,y) = Ln ( x² + y³ )     x = 2 , y = 1
 

-Again with h = 0.1

    << 3  ^  SWAP  SQ +  LN   >>   ENTER
                              2                         ENTER
                              1                         ENTER
                             0.1                       MDV4       >>>>  f⁽⁴⁾_xxyy =  ¶⁴f / ¶x²¶y² =  - 0.0383994555318

-Exact result = - 0.0384
 

Gaussian Curvature & Mean Curvature of a Surface
 

  KGM calculates the Gaussian curvature K_G and the mean curvature K_m of a surface defined by a function  z = z(x,y)

-The Gaussian Curvature K_G is be computed by:

   K_G = [ ( ¶²z / ¶x² ) ( ¶²z / ¶y² ) - ( ¶²z / ¶x¶y )² ] / [ 1 + ( ¶z / ¶x )² + ( ¶z / ¶y )² ]²

-The mean curvature K_m is given by:

    K_m = (1/2) [ t ( 1+p² ) + r ( 1+q² ) - 2 p q s ] / ( 1+p²+q² )^3/2

   where  p = ¶z / ¶x , q = ¶z / ¶y,  r = ¶²z / ¶x²,  s = ¶²z / ¶x¶y,  t = ¶²z / ¶y²
 
 
 

STACK	INPUTS	OUTPUTS
Level 4	<< f(x,y) >>	/
Level 3	x	/
Level 2	y	KG
Level 1	h	Km

   Where the program in level 4 takes x & y in level 2 & 1 respectively and returns f(x,y) in level 1

Example:    A surface is defined by  z(x,y) = Ln ( 1 + x² y³ )    Calculate the Gaussian & mean curvatures at the point  (x,y) = (1,1)

-With h = 0.1

   <<  3  ^  SWAP  SQ  *  1  +  LN  >>     ENTER
                                  1                               ENTER
                                  1                               ENTER                KG :  -0.124572639507
                                 0.1                              KGM     >>>>    KM :  -0.171207007011
 

Notes:

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

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

Curl of a 3D-Vector Field - Rectangular Coordinates
 

Formula:                Given a 3dim-vector field:   E = ( f , g , h )
 

  Curl E = ( ¶h/¶y - ¶g/¶z , ¶f/¶z - ¶h/¶x , ¶g/¶x - ¶f/¶y )
 
 
 

STACK	INPUTS	OUTPUTS
Level 7	<< f(x,y,z) >>	/
Level 6	<< g(x,y,z) >>	/
Level 5	<< h(x,y,z) >>	/
Level 4	x	/
Level 3	y	/
Level 2	z	/
Level 1	h	Curl E

 
Example:     E = ( f , g , h )  = [ exp(-x²) Ln(y²+z) , x²y²z² , exp(x) y²z ]     With  x = 1 , y = 2 , z = 3

-If you choose  h = 0.1

  <<  SWAP  SQ  +  LN  SWAP  SQ  NEG  EXP  *  >>  ENTER
                            <<  *  *  SQ  >>                                    ENTER
                 <<  SWAP  SQ  *  SWAP  EXP  *  >>            ENTER
                                         1                                                ENTER
                                         2                                                ENTER
                                         3                                                ENTER
                                       0.1                                                CuRL           >>>>   [   8.6193819414
                                                                                                                               -32.5668277359
                                                                                                                                 71.7897831765  ]
 

Divergence of a 3D-Vector Field - Rectangular Coordinates
 

Formula:                Given a 3dim-vector field:   E = ( f , g , h )
 

  Div E  =  ¶f/¶x + ¶g/¶y + ¶h/¶z
 
 

STACK	INPUTS	OUTPUTS
Level 7	<< f(x,y,z) >>	/
Level 6	<< g(x,y,z) >>	/
Level 5	<< h(x,y,z) >>	/
Level 4	x	/
Level 3	y	/
Level 2	z	/
Level 1	h	Div E

 
Example:     E = ( f , g , h )  = [ exp(-x²) Ln(y²+z) , x²y²z² , exp(x) y²z ]     With  x = 1 , y = 2 , z = 3

-If you choose  h = 0.1

  <<  SWAP  SQ  +  LN  SWAP  SQ  NEG  EXP  *  >>  ENTER
                           <<  *  *  SQ  >>                                    ENTER
                 <<  SWAP  SQ  *  SWAP  EXP  *  >>            ENTER
                                         1                                                ENTER
                                         2                                                ENTER
                                         3                                                ENTER
                                       0.1                                                 DVG           >>>>    Div E = 45.4414066323
 

Gradient - Rectangular Coordinates
 

Formula:                Given a function of 3 variables f(x,y,z)
 

  Grad f = ( ¶f/¶x , ¶f/¶y , ¶f/¶z )
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< f(x,y,z) >>	/
Level 4	x	/
Level 3	y	/
Level 2	z	/
Level 1	h	Grad f

 
Example:     f(x,y,z) = exp(-x²) Ln(y²+z)      With  x = 1 , y = 2 , z = 3

-If you choose  h = 0.1

  <<  SWAP  SQ  +  LN  SWAP  SQ  NEG  EXP  *  >>   ENTER
                                         1                                                ENTER
                                         2                                                ENTER
                                         3                                                ENTER
                                       0.1                                                GRaD           >>>>  [ -1.43172068193
                                                                                                                               0.210216823546
                                                                                                                               5.25542059036 E-2 ]
 

Laplacian - Rectangular Coordinates
 

Formula:                Given a function of 3 variables f(x,y,z)
 

  Lapl f  = ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z²
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< f(x,y,z) >>	/
Level 4	x	/
Level 3	y	/
Level 2	z	/
Level 1	h	Lapl f

 
Example:     f(x,y,z) = exp(-x²) Ln(y²+z)      With  x = 1 , y = 2 , z = 3

-If you choose  h = 0.1

  <<  SWAP  SQ  +  LN  SWAP  SQ  NEG  EXP  *  >>   ENTER
                                         1                                                ENTER
                                         2                                                ENTER
                                         3                                                ENTER
                                       0.1                                                LaPL           >>>>  Lapl f = 1.40919744606

Note:

-The exact result is     1.40919744532.....
 

3D-Vector Laplacian - Rectangular Coordinates
 

-With rectangular coordinates, the coordinates of the vector Lplacian are simply the Laplacian of the coordinates.
-It's not so easy with cylindrical or spherical coordinates.
 
 

STACK	INPUTS	OUTPUTS
Level 7	<< f(x,y,z) >>	/
Level 6	<< g(x,y,z) >>	/
Level 5	<< h(x,y,z) >>	/
Level 4	x	/
Level 3	y	/
Level 2	z	/
Level 1	h	Lapl E

 
Example:     E = ( f , g , h )  = [ exp(-x²) Ln(y²+z) , x²y²z² , exp(x) y²z ]     With  x = 1 , y = 2 , z = 3

-If you choose  h = 0.1

  <<  SWAP  SQ  +  LN  SWAP  SQ  NEG  EXP  *  >>  ENTER
                            <<  *  *  SQ  >>                                    ENTER
                 <<  SWAP  SQ  *  SWAP  EXP  *  >>            ENTER
                                         1                                                ENTER
                                         2                                                ENTER
                                         3                                                ENTER
                                       0.1                                               VLaP           >>>>   [   1.40919744606
                                                                                                                               98
                                                                                                                               48.9290729151  ]
 

Curl of a 3D-Vector Field - Cylindrical Coordinates
 

Formula:                Given a 3dim-vector field:   E = ( f , g , h )
 

  Curl E = [ (1/r) ¶h/¶f - ¶g/¶z , ¶f/¶z - ¶h/¶r , (1/r) ( g + r ¶g/¶r - ¶f/¶f ) ]
 
 
 

STACK	INPUTS	OUTPUTS
Level 7	<< f(r,u,z) >>	/
Level 6	<< g(r,u,z) >>	/
Level 5	<< h(r,u,z) >>	/
Level 4	r	/
Level 3	u	/
Level 2	z	/
Level 1	h	Curl E

 
Example:         E = ( f , g , h )  = (  r z² sin²f , r² z , r³ z cos f )                   r = 2 , f = PI/5 , z = 1

-If you choose  h = 0.1

        <<  SWAP  SIN  *  SQ  *  >>           ENTER
  <<  SWAP  DROP  SWAP  SQ  *  >>    ENTER
 <<  SWAP  COS  *  SWAP  3  ^  *  >>   ENTER
                              2                                   ENTER
                          PI  5  /
                              1                                   ENTER
                            0.1                                   CuRL    >>>>   [ -6.35114100918
                                                                                                -8.32623792129
                                                                                                 5.04894348373 ]

Note:

 RAD mode is automatically set by this program
 

Divergence of a 3D-Vector Field - Cylindrical Coordinates
 

Formula:                Given a 3dim-vector field:   E = ( f , g , h )
 

  Div E  =  ¶f/¶r + (1/r) f + (1/r) ¶g/¶f + ¶h/¶z
 
 

STACK	INPUTS	OUTPUTS
Level 7	<< f(r,u,z) >>	/
Level 6	<< g(r,u,z) >>	/
Level 5	<< h(r,u,z) >>	/
Level 4	r	/
Level 3	u	/
Level 2	z	/
Level 1	h	Div E

 
Example:         E = ( f , g , h )  = (  r z² sin²f , r² z , r³ z cos f )                   r = 2 , f = PI/5 , z = 1

-If you choose  h = 0.1

        <<  SWAP  SIN  *  SQ  *  >>           ENTER
  <<  SWAP  DROP  SWAP  SQ  *  >>    ENTER
 <<  SWAP  COS  *  SWAP  3  ^  *  >>   ENTER
                              2                                   ENTER
                          PI  5  /
                              1                                   ENTER
                            0.1                                    DVG    >>>>   7.16311896062

Note:

 RAD mode is automatically set by this program
 

Gradient - Cylindrical Coordinates
 

Formula:                Given a function of 3 variables f(r,f,z)

 Grad f = ( ¶f/¶r , (1/r) ¶f/¶f , ¶f/¶z )
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< f(r,u,z) >>	/
Level 4	r	/
Level 3	u	/
Level 2	z	/
Level 1	h	Grad f

 
Example:             f(r,f,z) =  r³ z cos f                 r = 2 , f = PI/5 , z = 1

-If you choose  h = 0.1

 <<  SWAP  COS  *  SWAP  3  ^  *  >>   ENTER
                              2                                   ENTER
                          PI  5  /
                              1                                   ENTER
                            0.1                                   GRaD    >>>>   [ 9.70820393254
                                                                                               -2.35114100918
                                                                                                6.472135955 ]

Note:

 RAD mode is automatically set by this program
 

Laplacian - Cylindrical Coordinates
 

Formula:       Given a function of 3 variables f(r,f,z)

  Lapl f  = ¶²f / ¶r² + (1/r) ¶f/¶r + (1/r²) ¶²f / ¶f² + ¶²f / ¶z²
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< f(r,u,z) >>	/
Level 4	r	/
Level 3	u	/
Level 2	z	/
Level 1	h	Lapl f

 
Example:             f(r,f,z) =  r³ z cos f                 r = 2 , f = PI/5 , z = 1

-If you choose  h = 0.1

 <<  SWAP  COS  *  SWAP  3  ^  *  >>   ENTER
                              2                                   ENTER
                          PI  5  /
                              1                                   ENTER
                            0.1                                   LaPL    >>>>   12.9442719039

Note:

 RAD mode is automatically set by this program
 

3D-Vector Laplacian - Cylindrical Coordinates
 

-The formula is given here
 
 

STACK	INPUTS	OUTPUTS
Level 7	<< f(r,u,z) >>	/
Level 6	<< g(r,u,z) >>	/
Level 5	<< h(r,u,z) >>	/
Level 4	r	/
Level 3	u	/
Level 2	z	/
Level 1	h	Lapl E

 
Example:         E = ( f , g , h )  = (  r z² sin²f , r² z , r³ z cos f )                   r = 2 , f = PI/5 , z = 1

-If you choose  h = 0.1

        <<  SWAP  SIN  *  SQ  *  >>           ENTER
  <<  SWAP  DROP  SWAP  SQ  *  >>    ENTER
 <<  SWAP  COS  *  SWAP  3  ^  *  >>   ENTER
                              2                                   ENTER
                          PI  5  /
                              1                                   ENTER
                            0.1                                   VLaP    >>>>   [  1.69098300707
                                                                                                 3.95105651627
                                                                                                12.9442719039 ]

Note:

 RAD mode is automatically set by this program
 

Curl of a 3D-Vector Field - Spherical Coordinates
 

Formula:                Given a 3dim-vector field:   E = ( f , g , h )
 

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

STACK	INPUTS	OUTPUTS
Level 7	<< f(r,u,v) >>	/
Level 6	<< g(r,u,v) >>	/
Level 5	<< h(r,u,v) >>	/
Level 4	r	/
Level 3	u	/
Level 2	v	/
Level 1	h	Curl E

 
Example:         E = ( f , g , h )  = (  r sin² q cos² f , r²sin f  , r³ cos q cos² f )                   r = 2 , q = PI/3 , f = PI/5

-If you choose  h = 0.1

           <<  COS  SWAP  SIN  *  SQ  *  >>               ENTER
     <<  SIN  SWAP  DROP  SWAP  SQ  *  >>         ENTER
 <<  COS  SQ  SWAP  COS  *  SWAP  3  ^  *  >>   ENTER
                                   2                                                ENTER
                                PI  3  /
                                PI  5  /
                                  0.1                                                CuRL    >>>>   [ -3.37986734613
                                                                                                                  -6.05970708098
                                                                                                                    2.95989052817 ]

Note:

 RAD mode is automatically set by this program
 

Divergence of a 3D-Vector Field - Spherical Coordinates
 

Formula:                Given a 3dim-vector field:   E = ( f , g , h )
 

  Div E  = ¶f/¶r + (2/r) f + g (cos q)/(r sin q)  + (1/r) ¶g/¶q + (1/(r sin q)) ¶h/¶f
 
 
 

STACK	INPUTS	OUTPUTS
Level 7	<< f(r,u,v) >>	/
Level 6	<< g(r,u,v) >>	/
Level 5	<< h(r,u,v) >>	/
Level 4	r	/
Level 3	u	/
Level 2	v	/
Level 1	h	Div E

 
Example:         E = ( f , g , h )  = (  r sin² q cos² f , r²sin f  , r³ cos q cos² f )                   r = 2 , q = PI/3 , f = PI/5

-If you choose  h = 0.1

           <<  COS  SWAP  SIN  *  SQ  *  >>               ENTER
     <<  SIN  SWAP  DROP  SWAP  SQ  *  >>         ENTER
 <<  COS  SQ  SWAP  COS  *  SWAP  3  ^  *  >>   ENTER
                                   2                                                ENTER
                                PI  3  /
                                PI  5  /
                                  0.1                                                DVG    >>>>  -0.045010876787

Note:

 RAD mode is automatically set by this program
 

Gradient - Spherical Coordinates
 

Formula:                Given a function of 3 variables f(r,q,f)
 

 Grad f = ( ¶f/¶r , (1/r) ¶f/¶q , (1/(r sin q))  ¶f/¶f )
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< f(r,u,v) >>	/
Level 4	r	/
Level 3	u	/
Level 2	v	/
Level 1	h	Grad f

 
Example:      f(r,q,f)   =  r³ cos q cos²f                    r = 2 , q = PI/3 , f= PI/5

-If you choose  h = 0.1

 <<  COS  SQ  SWAP  COS  *  SWAP  3  ^  *  >>   ENTER
                                   2                                                ENTER
                                PI  3  /
                                PI  5  /
                                  0.1                                                GRaD    >>>>   [ 3.92705098312
                                                                                                                  -2.26728394224
                                                                                                                  -2.19637094272 ]

Note:

 RAD mode is automatically set by this program
 

Laplacian - Spherical Coordinates
 

Formula:                Given a function of 3 variables f(r,q,f)
 

  Lapl f  = ¶²f / ¶r² + (2/r) ¶f/¶r + (1/r²) ¶²f / ¶q² + (1/(r²tan q))  ¶f / ¶q + (1/(r²sin²q)) ¶²f / ¶f²
 
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< f(r,u,v) >>	/
Level 4	r	/
Level 3	u	/
Level 2	v	/
Level 1	h	Lapl f

 
Example:      f(r,q,f)   =  r³ cos q cos²f                    r = 2 , q = PI/3 , f = PI/5

-If you choose  h = 0.1

 <<  COS  SQ  SWAP  COS  *  SWAP  3  ^  *  >>   ENTER
                                   2                                                ENTER
                                PI  3  /
                                PI  5  /
                                  0.1                                                LaPL    >>>>   5.72103965058

Note:

 RAD mode is automatically set by this program
 

3D-Vector Laplacian - Spherical Coordinates
 

-The formula, even more complicated than the cylindrical case, is given here
 
 

STACK	INPUTS	OUTPUTS
Level 7	<< f(r,u,v) >>	/
Level 6	<< g(r,u,v) >>	/
Level 5	<< h(r,u,v) >>	/
Level 4	r	/
Level 3	u	/
Level 2	v	/
Level 1	h	Lapl E

 
Example:         E = ( f , g , h )  = (  r sin² q cos² f , r²sin f  , r³ cos q cos² f )                   r = 2 , q = PI/3 , f = PI/5

-If you choose  h = 0.1

           <<  COS  SWAP  SIN  *  SQ  *  >>               ENTER
     <<  SIN  SWAP  DROP  SWAP  SQ  *  >>         ENTER
 <<  COS  SQ  SWAP  COS  *  SWAP  3  ^  *  >>   ENTER
                                   2                                                ENTER
                                PI  3  /
                                PI  5  /
                                  0.1                                                VLaP    >>>>   [  1.04501087529
                                                                                                                   3.79418051536
                                                                                                                   5.10341187668 ]

Note:

 RAD mode is automatically set by this program
 

Mixed Derivative d³F/dxdydz
 

 "MDV3" employs a method of order 11 to estimate  f'''_xyz = ¶³f / ¶x¶y¶z

    h³ f'''_xyz ~  SUM_k=1,2,...,5  A_k [ f_k00 - f_-k00 + f_0k0 - f_0-k0 + f_00k - f_00-k + f_kkk - f_-k-k-k - ( f_kk0 - f_-k-k0 + f_k0k - f_-k0-k + f_0kk - f_0-k-k ) ]

 With    A₁ = 5/6 ,  A₂ = -5/84  ,  A₃ = +5/756  ,  A₄ = -5/8064  ,  A₅ = +1/31500
 
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< f(x,y,z) >>	/
Level 4	x	/
Level 3	y	/
Level 2	z	/
Level 1	h	f'''xyz

 
Example:     f(x,y,z) = exp(-x²) Ln(y²+z)      With  x = 1 , y = 2 , z = 3

-If you choose  h = 0.1

  <<  SWAP  SQ  +  LN  SWAP  SQ  NEG  EXP  *  >>   ENTER
                                         1                                                ENTER
                                         2                                                ENTER
                                         3                                                ENTER
                                       0.1                                               MDV3   >>>>   0.0600619522

-The exact result is   0.0600619495....
 

Biharmonic Operator ( 3-Dim )
 

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

       D² f = ¶⁴f / ¶x⁴ + ¶⁴f / ¶y⁴ + ¶⁴f / ¶z⁴ + 2 ( ¶⁴f / ¶x²¶y² + ¶⁴f / ¶x²¶z² + ¶⁴f / ¶y²¶z² )
 

-"BHRM" uses the following approximation
 

   h⁴ D² f  ~  A f₀₀₀
               + B₁ ( f₁₀₀ + f_-100 + f₀₁₀ + f_0-10 + f₀₀₁ + f_00-1 )  + C₁ ( f₁₁₀ + f_-1-10 + f_1-10 + f_-110 + f₁₀₁ + f_-10-1 + f_10-1 + f_-101 + f₀₁₁ + f_0-1-1 + f_01-1 + f_0-11 )
               + B₂ ( f₂₀₀ + f_-200 + f₀₂₀ + f_0-20 + f₀₀₂ + f_00-2 )  + C₂ ( f₂₂₀ + f_-2-20 + f_2-20 + f_-220 + f₂₀₂ + f_-20-2 + f_20-2 + f_-202 + f₀₂₂ + f_0-2-2 + f_02-2 + f_0-22 )
               + B₃ ( f₃₀₀ + f_-300 + f₀₃₀ + f_0-30 + f₀₀₃ + f_00-3 )  + C₃ ( f₃₃₀ + f_-3-30 + f_3-30 + f_-330 + f₃₀₃ + f_-30-3 + f_30-3 + f_-303 + f₀₃₃ + f_0-3-3 + f_03-3 + f_0-33 )
               + B₄ ( f₄₀₀ + f_-400 + f₀₄₀ + f_0-40 + f₀₀₄ + f_00-4 )  + C₄ ( f₄₄₀ + f_-4-40 + f_4-40 + f_-440 + f₄₀₄ + f_-40-4 + f_40-4 + f_-404 + f₀₄₄ + f_0-4-4 + f_04-4 + f_0-44 )
               + B₅ ( f₅₀₀ + f_-500 + f₀₅₀ + f_0-50 + f₀₀₅ + f_00-5 )  + C₅ ( f₅₅₀ + f_-5-50 + f_5-50 + f_-550 + f₅₀₅ + f_-50-5 + f_50-5 + f_-505 + f₀₅₅ + f_0-5-5 + f_05-5 + f_0-55 )

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

     A = 41523361 / 540000  ,   B₁ =  -4069 / 180             C₁ =  +10 / 3
                                                 B₂ =  +4969 / 1260          C₂ =  -10 / 84
                                                 B₃ =  -2201 / 3240           C₃ =  +10 / 1134
                                                 B₄ =  +371 / 4320            C₄ =  -10 / 16128
                                                 B₅ =  -5221 / 945000       C₅ =  +2 / 78750
 
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< f(x,y,z) >>	/
Level 4	x	/
Level 3	y	/
Level 2	z	/
Level 1	h	f'''xyz

 
Example:     f(x,y,z) = exp(-x²) Ln(y²+z)      With  x = 1 , y = 2 , z = 3

-If you choose  h = 0.1

  <<  SWAP  SQ  +  LN  SWAP  SQ  NEG  EXP  *  >>   ENTER
                                         1                                                ENTER
                                         2                                                ENTER
                                         3                                                ENTER
                                       0.1                                               BHRM   >>>>   -14.3426429833

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

Notes:

-91 evaluations of the function are performed (!).
-The formula is of order 10 but - due to cancellation of leading digits - there are seldom more than 7 or 8 significant digits.
 

Triharmonic Operator ( 3-Dim )
 

-The triharmonic operator is the trilaplacian operator = the Laplacian of the Laplacian of the Laplacian of a function f

-"THRM" employs a method of order 10  to evaluate
 

   D³ f = ¶⁶f / ¶x⁶ + ¶⁶f / ¶y⁶ + ¶⁶f / ¶z⁶ + 3 ( ¶⁶f / ¶x⁴¶y² + ¶⁶f / ¶x⁴¶z² + ¶⁶f / ¶y⁴¶z² + ¶⁶f / ¶x²¶y⁴ + ¶⁶f / ¶x²¶z⁴ + ¶⁶f / ¶y²¶z⁴ ) + 6 ¶⁶f / ¶x²¶y²¶z²
 

   h⁶ D³ f  ~  SUM_m=1....5  A_m ( f_m00 + f_-m00 + f_0m0 + f_0-m0 + f_00m + f_00-m - 6 f₀₀₀ )
                                     + B_m ( f_mm0 + f_-m-m0 + f_m-m0 + f_-mm0 + f_m0m + f_-m0-m + f_m0-m + f_-m0m + f_0mm + f_0-m-m + f_0m-m + f_0-mm - 12 f₀₀₀ )
                                     + C_m ( f_mmm + f_-m-m-m + f_mm-m + f_-m-mm + f_m-mm + f_-mm-m + f_m-m-m + f_-mmm - 8 f₀₀₀ )

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

     A₁ = +22997 / 120                  B₁ = - 2869 / 60                C₁ = + 10
     A₂ = - 12709 / 420                  B₂ = + 667 / 240               C₂ = - 5 / 56
     A₃ = + 858391 / 136080         B₃ = - 15007 / 68040        C₃ = + 5 / 1701
     A₄ = - 137843 / 161280          B₄ = + 5119 / 322560       C₄ = - 5 / 43008
     A₅ = + 894317 / 15750000     B₅ = - 739 / 1125000        C₅ = + 1 / 328125
 
 
 
 

STACK	INPUTS	OUTPUTS
Level 5	<< f(x,y,z) >>	/
Level 4	x	/
Level 3	y	/
Level 2	z	/
Level 1	h	f'''xyz

 
Example:     f(x,y,z) = exp(-x²) Ln(y²+z)      With  x = 1 , y = 2 , z = 3

-If you choose  h = 0.1

  <<  SWAP  SQ  +  LN  SWAP  SQ  NEG  EXP  *  >>   ENTER
                                         1                                                ENTER
                                         2                                                ENTER
                                         3                                                ENTER
                                       0.1                                               THRM   >>>>   133.533181
 

Notes:

-Fourth derivatives are already difficult to evaluate numerically but here,
  there will be seldom more than 5 or 6 significant digits, and sometimes only 3 or 4 !
-The function f is evaluated 131 times.
 

First Derivatives of a Function of N variables
 

  FD employs formulas of order 10,
  but in the subdirectories "Riemann" & "Manifolds",
   formulas of order 10 are used if Re(h) > 0
  and formulas of order 4 if Re(h) < 0  i-e

    df/dx   = (1/12h).[ f(x-2h) - 8.f(x-h) + 8.f(x+h) - f(x+2h) ]
 
 
 

STACK	INPUTS	OUTPUTS
Level 2	<< f(x1,....,xn) >>	/
Level 1	i	¶f / ¶xi

 
Example:      f(x) = exp(-x²).Ln(y²+z)    x = 1 , y = 1 , z = 1    with  h = 0.1
 

1°)  Store  n = 3  into  'N'   and     0.1  into  'h'

2°)  Store 1   1   1    into   'X1'  'X2'  'X3'    respectively

3°)    <<  X1  SQ  NEG  EXP  X2  SQ  X3  +  LN  *  >>   ENTER
                                             1                                                 FD           >>>>    f '_x1 = ¶f / ¶x₁ = -0.509989196833

-Likewise:

        <<  X1  SQ  NEG  EXP  X2  SQ  X3  +  LN  *  >>   ENTER
                                             2                                                 FD           >>>>    f '_x2 = ¶f / ¶x₂ = 0.367879439374

-Likewise:

        <<  X1  SQ  NEG  EXP  X2  SQ  X3  +  LN  *  >>   ENTER
                                             3                                                 FD           >>>>    f '_x3 = ¶f / ¶x₃ = 0.183939720617

Note:

-The contents of variables 'X1'  'X2'  ....  are temporarily modified during the calculations, but they are restored finally.
-So, if you interrupt the program before the end, check these variables before performing another calculation
 

Second Derivatives of a Function of N variables
 

  SD employs formulas of order 10,
  but in the subdirectories "Riemann" & "Manifolds",
  formulas of order 10 are used if Re(h) > 0
  and formulas of order 4 if Re(h) < 0  i-e

          d²f/dx² = (1/12h²).[ - f(x-2h) + 16.f(x-h) - 30.f(x) +16.f(x+h) - f(x+2h) ]

 and, for the crossed derivative

           f "_xy = ¶²f / ¶x¶y = (1/24.h²).[ 30 f₀₀ + 16.( f₁₁ + f _-1-1 - f₁₀ - f _-10 - f₀₁ - f_0-1 ) + ( - f ₂₂ - f _-2-2 + f₂₀ + f _-20 + f₀₂ + f_0-2 ) ]

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

STACK	INPUTS	OUTPUTS
Level 3	<< f(x1,....,xn) >>	/
Level 2	i	/
Level 1	j	¶2f / ¶xi¶xj

 
Example:       f(x) = exp(-x²).Ln(y²+z)    x = 1 , y = 1 , z = 1    with  h = 0.1
 

1°)  Store  n = 3  into  'N'   and     0.1  into  'h'

2°)  Store 1   1   1    into   'X1'  'X2'  'X3'    respectively

3°)    <<  X1  SQ  NEG  EXP  X2  SQ  X3  +  LN  *  >>   ENTER
                                             1                                             ENTER
                                             1                                                 SD           >>>>    f "_xx = ¶²f / ¶x² = 0.509989195329

-Likewise:

        <<  X1  SQ  NEG  EXP  X2  SQ  X3  +  LN  *  >>   ENTER
                                             1                                             ENTER
                                             2                                                 SD           >>>>    f "_xy = ¶²f / ¶x¶y = -0.735758876115

Note:

-The contents of variables 'X1'  'X2'  ....  are temporarily modified during the calculations, but they are restored finally.
-So, if you interrupt the program before the end, check these variables before performing another calculation
 

Third Derivatives of a Function of N variables
 

   TD employs the same formulas as d34F , MDV3 , MD12 ,
 
 
 

STACK	INPUTS	OUTPUTS
Level 4	<< f(x1,....,xn) >>	/
Level 3	i	/
Level 2	j	/
Level 1	k	¶3f / ¶xi¶xj¶xk

 
Example:           f(x,y,z,t) = exp(-x²) Ln( y² + z + t³ )     x = 1 , y = 2 , z = 3 , t = 1       h = 0.1
 

1°)  Store  n = 4  into  'N'   and     0.1  into  'h'

2°)  Store 1   2   3   1    into   'X1'  'X2'  'X3'   'X4'    respectively then, for instance:

  •  f'''_xyz = ¶³f / ¶x¶y¶z  = ?

        <<  X1  SQ  NEG  EXP  X2  SQ  X3  +  X4  3  ^  +  LN  *  >>   ENTER
                                                      1                                                       ENTER
                                                      2                                                       ENTER
                                                      3                                                          TD           >>>>     f'''_xyz = ¶³f / ¶x¶y¶z  = 0.045984928498

  •  f'''_xtt = ¶³f / ¶x¶t²  = ?

        <<  X1  SQ  NEG  EXP  X2  SQ  X3  +  X4  3  ^  +  LN  *  >>   ENTER
                                                      1                                                       ENTER
                                                      4                                                       ENTER
                                                      4                                                          TD           >>>>    f'''_xtt = ¶³f / ¶x¶t²  = - 0.448353075973

  •  f'''_yyy = ¶³f / ¶y³  = ?

        <<  X1  SQ  NEG  EXP  X2  SQ  X3  +  X4  3  ^  +  LN  *  >>   ENTER
                                                      2                                                       ENTER
                                                      2                                                       ENTER
                                                      2                                                          TD           >>>>     f'''_yyy = ¶³f / ¶y³  = - 0.0459849295635
 

Note:

-The contents of variables 'X1'  'X2'  ....  are temporarily modified during the calculations, but they are restored finally.
-So, if you interrupt the program before the end, check these variables before performing another calculation
 

Curl of an NDim-Vector Field - Rectangular Coordinates
 

-If E = ( f₁ , ......... , f_n ) ,  Curl E  is an antisymmetric tensor defined as  c_ij = ¶f _i/ ¶x_j  - ¶f _j/ ¶x_i   ( some authors use the opposite sign )

 CuRL  returns the antisymmetric matrix  [ c_ij ]
 
 
 

STACK	INPUTS	OUTPUTS
Level n	<< f1( x1 , ..... , xn ) >>	/
Level ...	.................	/
Level 2	<< fn-1( x1 , ..... , xn ) >>	/
Level 1	<< fn( x1 , ..... , xn ) >>	Curl  E

 
Example1:     E = ( f , g , h )  = [ exp(-x²) Ln(y²+z) , x²y²z² , exp(x) y²z ]     With  x = 1 , y = 2 , z = 3   and   h = 0.1

  Store  3  into  'N'   and  0.1  into  'h'
  Store  1  2  3  into  'X1'  'X2'  'X3'  respectively

  <<  X1  SQ  NEG  EXP  X2  SQ  X3  +  LN  *  >>   ENTER
               <<  X1  X2  X3  *  *  SQ  >>                        ENTER
         <<  X1  EXP  X2  SQ  X3  *  *  >>                     CuRL     gives the following matrix:

   [ [            0               71.7897831765   32.5668277359  ]
     [ -71.7897831765             0                8.6193819414   ]
     [ -32.5668277359  -8.6193819414              0              ] ]

-And the 3 independant components  [ c₂₃  c₃₁  c₁₂ ] =  [   8.6193819414     -32.5668277359      71.7897831765  ]
  form the usual Curl that we've found above - which is actually a pseudo-vector.

Example2:     E = ( u , v , w , p ) = ( x²y²z²t² ,  x y z³ t ,  x²y z t³ , xy + zt )     x = t = 1 , y = 2 , z = 3    with  h = 0.1

  Store 4 into 'N'  and  0.1  into  'h'
  Store 1 into  'X1'  &  'X4'  and  2  into 'X2'  &  3  into  'X3'

     <<  X1  X2  X3  X4  *  *  *  SQ  >>     ENTER
     <<  X1  X2  X3  3  ^  X4  *  *  *  >>     ENTER
  <<  X1  SQ  X2  X3  X4  3  ^ *  *  *  >>  ENTER
        << X1  X2  *  X3  X4  *  +  >>             CuRL    returns the matrix:

     [ [    0   18  -12  -70 ]
       [  -18   0   -51  -53 ]
       [    12  51    0   -17 ]
       [    70  53   17    0  ] ]
 

Divergence of an NDim-Vector Field - Rectangular Coordinates
 

-If E = ( f₁ , ......... , f_n ) ,  Div E = ¶f ₁/ ¶x₁  + ................. + ¶f _n/ ¶x_n
 
 

STACK	INPUTS	OUTPUTS
Level n	<< f1( x1 , ..... , xn ) >>	/
Level ...	.................	/
Level 2	<< fn-1( x1 , ..... , xn ) >>	/
Level 1	<< fn( x1 , ..... , xn ) >>	Div  E

 
Example:     E = ( u , v , w , p ) = ( x²y²z²t² ,  x y z³ t ,  x²y z t³ , xy + zt )     x = t = 1 , y = 2 , z = 3    with  h = 0.1

  Store 4 into 'N'  and  0.1  into  'h'
  Store 1 into  'X1'  &  'X4'  and  2  into 'X2'  &  3  into  'X3'

     <<  X1  X2  X3  X4  *  *  *  SQ  >>     ENTER
     <<  X1  X2  X3  3  ^  X4  *  *  *  >>     ENTER
  <<  X1  SQ  X2  X3  X4  3  ^ *  *  *  >>  ENTER
        << X1  X2  *  X3  X4  *  +  >>             DVG    >>>>   Div E = 104
 

Gradient of a Function of N Variables
 

    Grad f = ( ¶f ₁/ ¶x₁ , .............. , ¶f _n/ ¶x_n )
 
 

STACK	INPUTS	OUTPUTS
Level 1	<< f(x1,....,xn) >>	Grad f

 
Example:       exp(-x²) [ Ln(y²+z) + 1/t ]     x = t = 1 , y = 2 , z = 3    with  h = 0.1

  Store 4 into 'N'  and  0.1  into  'h'
  Store 1 into  'X1'  &  'X4'  and  2  into 'X2'  &  3  into  'X3'

  << X1  SQ  NEG  EXP  X2  SQ  X3  +  LN  X4  INV  +  *  >>    GRad   gives the 4-vector:

     [ -2.16747956713  0.210216823496  0.0525542059107  -0.367880413465 ]
 

Hessian Matrix
 

  HeSS calculates the coefficients of the Hessian matrix H defined by    H = [ h_i,j ]   with  h_i,j = ¶²f / ¶x_i¶x_j
  where f is a function of n variables
 
 
 

STACK	INPUTS	OUTPUTS
Level 1	<< f(x1,....,xn) >>	Hess f

 
Example:       exp(-x²) [ Ln(y²+z) + 1/t ]     x = t = 1 , y = 2 , z = 3    with  h = 0.1

  Store 4 into 'N'  and  0.1  into  'h'
  Store 1 into  'X1'  &  'X4'  and  2  into 'X2'  &  3  into  'X3'

  << X1  SQ  NEG  EXP  X2  SQ  X3  +  LN  X4  INV  +  *  >>    HeSS  returns in 25 seconds the 4x4 symmetric matrix:

                      [ [  2.16747956102    -0.420433648869    -0.105108411936    0.73576055319   ]
                        [ -0.420433648869  -0.0150154894179  -0.0300309737861     7.566.. E-10     ]
                        [ -0.105108411936  -0.0300309737861  -0.00750774532488   8.316...E-10     ]
                       [   0.73576055319       7.566.. E-10               8.316...E-10       0.735760824916 ] ]

-The small coefficients of the order of E-10 should be 0, but due to roundoff-errors...
 

Laplacian of a Function of N variables
 

 "LAPL" may be used to evaluate the Laplacian of a function of n variables.
 
 
 

STACK	INPUTS	OUTPUTS
Level 1	<< f(x1,....,xn) >>	Lapl f

 
Example:       exp(-x²) [ Ln(y²+z) + 1/t ]     x = t = 1 , y = 2 , z = 3    with  h = 0.1

  Store 4 into 'N'  and  0.1  into  'h'
  Store 1 into  'X1'  &  'X4'  and  2  into 'X2'  &  3  into  'X3'

  << X1  SQ  NEG  EXP  X2  SQ  X3  +  LN  X4  INV  +  *  >>    LaPL  >>>>  2.8807171512

 Exact result = 2.88071520999
 

Biharmonic Operator ( N-Dim )
 

-Formulas of order 10 are now used to estimate the n-dimensional biharmonic operator
 

         D² f = SUM_k=1,....,n  ¶⁴f / ¶x_k⁴ + 2  SUM_j<k  ¶⁴f / ¶x_j²¶y_k²
 
 
 

STACK	INPUTS	OUTPUTS
Level 1	<< f(x1,....,xn) >>	D2 f

 
Example:     f(x,y,z,t) = exp(-x²) Ln( y² + z + t³ )     x = 1 , y = 2 , z = 3 , t = 1         with  h = 0.1

  Store 4 into 'N'  and  0.1  into  'h'
  Store 1 into  'X1'  &  'X4'  and  2  into 'X2'  &  3  into  'X3'

  << X1  SQ  NEG  EXP  X2  SQ  X3  +  X4  3  ^  +  LN  *  >>  BHRM  >>>   D² f (1,2,3,1) = -14.93970202      ( in 19 seconds )

-The exact result is  -14.9397000647...
 

Note:

-The function f is evaluated  10 n²  + 1  times
 
 

N	EVAL
1	11
2	41
3	91
4	161
5	251
6	361
7	491
8	641
9	811

 
Example2:    9 variables,     f (x,y,z,t,u,v,w,r,s)  = [ exp(-x y z) + t u v w ] / Ln ( 1 + w r s )    x = y = ....... = s = 1    again with  h = 0.1

  Store 9 into 'N'  and  0.1  into  'h'
  Store 1 into  'X1'  'X2'   .......   'X9'

  << X1  X2  X3  *  *  NEG  EXP  X4  X5  X6  X7  *  *  *  +  X7  X8  X9  *  *  1  +  LN  /  >>  BHRM   returns in 1mn51s

            D² f (1,......,1) = 103.210626858

Notes:

-The exact value is  103.2389124...
-We only have 5 significant digits;
 

Gaussian Curvature & Mean Curvature of a HyperSurface
 

 -Here, we assume that the hypersurface is defined implicitly by   f( x₁ , ........... , x_n ) = 0

Formulae:

   K_G = - ( det W ) / ( - | Grad f | )ⁿ⁺¹      where  W is defined in paragraph 1°) c)

   K_m = [ Grad f x H x ^TGrad f  -  | Grad f |² Trace (H) ] / (n-1) / | Grad f |³     where  H = Hessian matrix
 
 

STACK	INPUTS	OUTPUTS
Level 2	/	KG
Level 1	<< f(x1,....,xn) >>	Km

 
Example:    The surface ( n = 3 )   x⁴ y³ z² - 1 = 0   with  x = y = z = 1

  Store 3 into 'N'  and  0.1  into  'h'   ( if you choose h = 0.1 )
  Store 1 into  'X1'  'X2'  'X3'

                                                                                                                K_G  =  0.256837098693
  <<  X1  4  ^  X2  3  ^  X3  SQ  *  *  1  -  >>   ENTER  KGM  >>>    K_m  =  0.518666289391                   ( in 16 seconds )
 

Warning:

-This program does not check that  f(x₁, ..... ,x_n) = 0
 

Initialization - Metric Tensor & Christoffel Symbols
 

-Always execute "INIGC" before calculating the components of the tensors below
 

-You have to store the functions g_ij in the variables 'G11'  'G22'  .......  'G12'  'G13'  .....
-Since the metric tensor is symmetric, 'G21' ...  are unuseful.

 INIGC  stores  the coefficients g_ij under the form of a symmetric matrix into 'GIJ',  the inverse matrix into 'G^IJ'
 and the Christoffel symbols  G^k_ij = (1/2) g^km (  ¶_i g_mj + ¶_j g_im -  ¶_m g_ij )   are then stored in the variable 'CIJK' under the form of a list,
 only  G^k_ij with  i <= j  because they are symmetric  G^k_ji = G^k_ij
 
 
 

STACK	INPUT	OUTPUT
Level 1	/	"DONE"

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

    g₁₁ = 1 + x² y² z²                            g₁₂ = x y
    g₂₂ = 1 + x² + y² + z²                     g₁₃ = x z
    g₃₃ = 1 + x² y² + x² z² + y² z²         g₂₃ = y z

-Load for instance these 6 routines into 'G11'  'G22'  'G33'  'G12'  'G13'  'G23'  respectively

    <<  X1  X2  X3  *  *  SQ  1  +  >>
    <<  X1  SQ  X2  SQ  X3  SQ  +  +  1  +  >>
    <<  X1  X2  *  SQ  X1  X3  *  SQ  X2  X3  *  SQ  +  +  1  +  >>
    <<  X1  X2  *  >>
    <<  X1  X3  *  >>
    <<  X2  X3  *  >>

 [ If , say g₁₂ = 0 , you can store 0 in 'G12' instead of  << 0 >> ]

-If we want to calculate the metric tensor, the Christoffel symbols and other tensors at the point  x = y = z = 1  with  h = -0.01

   Store  -0.01  into  'h'  and  3  into  'N'
   Store  1  into  'X1'  'X2'  'X3'

-Then press   INIGC   the program places  "DONE"  in level 1   after 95 seconds and you get

                                             [ [ 2  1  1 ]                                       [ [  5/8  -1/8   -1/8   ]
   The metric tensor  in  GIJ  =  [ 1  4  1 ]       its inverse in  G^IJ      [ -1/8  7/24  -1/24 ]
                                                [ 1  1  4 ] ]                                      [ -1/8 -1/24   7/24 ] ]

   The list of the Christoffel symbols  G^k_ij   in CIJK = { 0.625  0.5  0.375  ....................  0.75 }

   G¹₁₁ = 5/8     G¹₁₂ =  1/2     G¹₁₃ =  3/8
                      G¹₂₂ = -1/8   G¹₂₃ = -3/8
                                             G¹₃₃ = -3/4

   G²₁₁ = -1/8     G²₁₂ = 1/6       G²₁₃ = -5/24
                        G²₂₂ = 7/24   G²₂₃ = 5/24
                                               G²₃₃ = -1/4

   G³₁₁ = -1/8     G³₁₂ = -1/6     G³₁₃ =  11/24
                        G³₂₂ = -1/24  G³₂₃ =  13/24
                                               G³₃₃ =   3/4
 

-Use  GETCIJK  below to recall these coefficients
 

Recalling Cijk
 

-After executing INIGC
 
 

STACK	INPUTS	OUTPUTS
Level 3	i	/
Level 2	j	/
Level 1	k	Gkij

 
Example:

   1  ENTER
   2  ENTER
   3  GETCIJK  gives  G³₁₂ = -1/6
 

Curvature Tensor ( 0-4 )
 

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

     R_ijkl = (1/2) ( ¶_jk g_il + ¶_il g_jk - ¶_jl g_ik - ¶_ik g_jl ) + g_mp ( G^m_jkG^p_il - G^m_jlG^p_ik )
 
 

STACK	INPUTS	OUTPUTS
Level 4	i	/
Level 3	j	/
Level 2	k	/
Level 1	l	Rijkl

 
Examples:

 The same metric and with x = y = z = 1

  1  ENTER
  2  ENTER
  1  ENTER
  2  RIJKL  >>>>  R₁₂₁₂ = -5/24

-Likewise

  2  ENTER
  3  ENTER
  1  ENTER
  3   RIJKL    >>>>  R₂₃₁₃ = -7/24                             ---Execution time = 7s---

Notes:

-The curvature tensor has many components but it has also many symmetries
-So there are in fact N =  n² ( n² - 1 ) / 12  independant components i-e

   N = 1  if  n = 2
   N = 6  if  n = 3
   N =20 if  n = 4
 

Curvature Tensor ( 1-3 )
 

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

      R^l_ijk = g^lm R_mijk

-"RIJK^L" calculates these components
 
 
 

STACK	INPUTS	OUTPUTS
Level 4	i	/
Level 3	j	/
Level 2	k	/
Level 1	l	Rlijk

 
Examples:

  2  ENTER
  1  ENTER
  2  ENTER
  1  RIJK^L    >>>>  R¹₂₁₂ = -13/96

-Likewise

  3  ENTER
  1  ENTER
  3  ENTER
  2  RIJK^L  >>>>  R²₃₁₃ = -26/288                                      ---Execution time = 22s---
 

Ricci Tensor
 

-If we contract 2 indices of the curvature tensor ( the 2nd and the 4th ), we get the Ricci tensor which is symmetric

    R_ij = g^km  R_ikjm

-We can also contract other indices, but it would yield  -R_ij
 
 

STACK	INPUTS	OUTPUTS
Level 2	i	/
Level 1	j	Rij

 
Examples:

  1  ENTER
  1    RIJ     >>>>  R₁₁ = 7/96                                        ---Execution time = 60s---                  ( 86 seconds with  h > 0 )

-Likewise

  2  ENTER
  3    RIJ      >>>>  R₂₃ = -277/288                                      ---Execution time = 64s---

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

   R₁₁ = +7/96       R₂₂ =  -79/288        R₃₃ = -5/144
   R₁₂ = -28/96      R₂₃ = -277/288
   R₁₃ = -17/96
 

Scalar Riemannian Curvature
 

-Contracting the Ricci tensor, we get the scalar curvature    R = g^ij  R_ij
-RSC calculates and stores R in the variable 'R'
 
 

STACK	INPUT	OUTPUT
Level 1	/	R

 
Example:

   RSC  >>>>   R = +11/72                                             ---Execution time = 6m51s---        It would be much slower with h > 0
 

Notes:

-This program is probably the slowest one.
-With a 2D-surface, the scalar curvature is twice the Gaussian curvature given by KGM
 

>>> Always execute RSC before computing Einstein tensor or Weyl tensor. 'R' must contain R
 

Einstein Tensor
 

-Einstein tensor is a symmetric tensor defined as

   E_ij = R_ij - (1/2) R g_ij

-Only after executing RSC - which stores R into 'R' - execute EIJ to get the components of Einstein tensor
 
 
 

STACK	INPUTS	OUTPUTS
Level 2	i	/
Level 1	j	Eij

 
Examples:

  1  ENTER
  1     EIJ      >>>>  the HP48 displays  "XEQ RSC FIRST"  and finally     E₁₁ = -23/288                       ---Execution time = 61s---

-Likewise

  2  ENTER
  3     EIJ      >>>>  E₂₃ = -299/288                                                                                            ---Execution time = 65s---

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

   E₁₁ = -23/288         E₂₂ = -167/288      E₃₃ = -49/144
   E₁₂ = -53/144         E₂₃ = -299/288
   E₁₃ = -73/288

Notes:

-When you XEQ "EIJ" the HP48 displays "XEQ RSC FIRST" to remind that the variable 'R' must contain the scalar curvature R before executing "EIJ"

-The cosmolical constant  L  is omitted here: the complete tensor is     E_ij = R_ij - (1/2) R g_ij + L g_ij

-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

    D_i Tⁱ_j = 0  for all j    where    D_i is the covariant differentiation.

-From  E_ij come the equations of General Relativity !
 

Weyl Tensor
 

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

   W_ijkl = R_ijkl + [ 1/(n-2) ] ( R_il g_jk - R_ik g_jl + R_jk g_il - R_jl g_ik ) + [ R/(n-1)/(n-2) ] ( g_ik g_jl - g_il g_jk )

-Like Einstein tensor, execute "RSC" first before executing "WIJKL"
 
 
 

STACK	INPUTS	OUTPUTS
Level 4	i	/
Level 3	j	/
Level 2	k	/
Level 1	l	Wijkl

 
Example:    If n < 4 , Weyl tensor is always equal to 0.

  1  ENTER^
  2  ENTER^
  1  ENTER^
  2  WIJKL  >>>>      W₁₂₁₂ = 0
 

-The number of independent components is n(n-3)(n+1)(n+2)/12

Notes:
 

-When you XEQ "WIJKL" ( if n > 4 ) the HP48 displays "XEQ RSC FIRST"  to remind that 'R' must contain the scalar curvature R before executing WIJKL

>>>>  If you press LASTARG just after executing WIJKL, you'll get R_ijkl in level 1 and another partial sum in level 2

-The number of independent components is n(n-3)(n+1)(n+2)/12

-With the same g_ij  if  i & j < 4 ,  and  g₄₄ = x y z t   ,  g₁₄ = g₂₄ = g₃₄ = 0  ,  t = 1

  after storing 4 in 'N'  1  in  'X4'  exetuting INIGC  and  SRC ,  you will find   W₁₂₃₂ =   35 / 576
 

Curl of a CovariantVector Field
 

-The partial derivatives of a vector is not a tensor   and  Curl_ij V is defined as   D_i V_j - D_j V_i
-But if we apply the formulae given in one of the paragraph below, the Christoffel symbols vanish and finally it yields

     Curl_ij V = ¶_i V_j - ¶_j V_i

 CuRL  returns the antisymmetric matrix  [ c_ij ]
 
 

STACK	INPUTS	OUTPUTS
Level n	<< f1( x1 , ..... , xn ) >>	/
Level ...	.................	/
Level 2	<< fn-1( x1 , ..... , xn ) >>	/
Level 1	<< fn( x1 , ..... , xn ) >>	Curl  V

 
Example:      V = ( u , v , w ) = ( x²y²z² ,  x y z³  ,  x²y z )     x = y = z = 1   with  h = 0.1

  Store   0.1  into  'h'

     <<  X1  X2  X3   *  *   SQ  >>     ENTER
     <<  X1  X2  X3  3  ^  *  *  >>      ENTER
     <<  X1  SQ  X2  X3  *  *  >>       CuRL    returns the matrix:

              [ [ 0  -1   0  ]
                [ 1   0   -2 ]
                [ 0   2    0 ] ]
 

Divergence of a Contravariant Vector Field
 

-The divergence of a vector field V is given by

    Div V =     D_i Vⁱ=    ¶_i Vⁱ + Gⁱ_ik V^k

-It is a scalar
 
 

STACK	INPUTS	OUTPUTS
Level n	<< f1( x1 , ..... , xn ) >>	/
Level ...	.................	/
Level 2	<< fn-1( x1 , ..... , xn ) >>	/
Level 1	<< fn( x1 , ..... , xn ) >>	Div V

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

    X = x² + y z³
    Y = x² y z²
    Z = x + y z

   << X1  SQ  X2  X3  3  ^  *  +  >>   ENTER
      <<  X1  X3  *  SQ  X2  * >>         ENTER
          << X1  X2  X3  *  + >>             DVG       >>>>    Div V = 10.5
 

Gradient of a Scalar Field
 

-The gradient of a scalar field f is defined as  Grad_k f  =  ¶_k f    i-e  simply the usual partial derivatives: it is a covariant vector
 
 

STACK	INPUTS	OUTPUTS
Level 1	<< f(x1,....,xn) >>	Grad f

 
Example:      Again the same metric at the same point ( 1 , 1 , 1 ) and the following scalar field:    f(x,y,z) = x² + y² z³
 

    << X1 SQ  X2  SQ  X3  3  ^  *  +  >>   GRaD   >>>   Grad f =  [ 2  2  3 ]
 

Laplacian of a Scalar Field
 

-The Laplacian of the same scalar field f  is a scalar defined by

    Lapl(f) = g^jk ( ¶_jk f  -  Gⁱ_jk ¶_i f  )
 
 

STACK	INPUTS	OUTPUTS
Level 1	<< f(x1,....,xn) >>	Lapl f

 
Example:      Again the same metric at the same point ( 1 , 1 , 1 ) and the following scalar field:    f(x,y,z) = x² + y² z³
 

    << X1 SQ  X2  SQ  X3  3  ^  *  +  >>   LaPL   >>>   317/96
 

Covariant Derivative of a Contravariant Vector Field
 

-The usual partial derivative of a vector or a tensor is not a tensor, but the covariant derivatives are.
-They involve the partial derivatives and the Christoffel symbols.
-The formulae to get a tensor are not the same if the vector or the tensor is covariant or contravariant:
 

    D_i V^j =    ¶_i V^j + G^j_im V^m

 COV^   returns the matrix whose elements are defined by this formula
 
 

STACK	INPUTS	OUTPUTS
Level n	<< f1( x1 , ..... , xn ) >>	/
Level ...	.................	/
Level 2	<< fn-1( x1 , ..... , xn ) >>	/
Level 1	<< fn( x1 , ..... , xn ) >>	[ Di Vj ]

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

    X = x² + y z³
    Y = x² y z²
    Z = x + y z

   << X1  SQ  X2  X3  3  ^  *  +  >>   ENTER
      <<  X1  X3  *  SQ  X2  * >>         ENTER
          << X1  X2  X3  *  + >>             COV^     gives

    [ [   4.5     1.5     1.5    ]
      [   9/8   49/24  41/24 ]               for instance   D₂ V³  = 41 / 24
      [ 15/8   31/24  95/24 ] ]
 

Covariant Derivative of a Covariant Vector Field
 

-The usual partial derivative of a vector or a tensor is not a tensor, but the covariant derivatives are.
-They involve the partial derivatives and the Christoffel symbols.
-The formulae to get a tensor are not the same if the vector or the tensor is covariant or contravariant:
 

    D_i V_j =    ¶_i V_j - G^m_ij V_m
 

 COV¯   returns the matrix whose elements are defined by this formula
 
 

STACK	INPUTS	OUTPUTS
Level n	<< f1( x1 , ..... , xn ) >>	/
Level ...	.................	/
Level 2	<< fn-1( x1 , ..... , xn ) >>	/
Level 1	<< fn( x1 , ..... , xn ) >>	[ Di Vj ]

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

    X = x² + y z³
    Y = x² y z²
    Z = x + y z

   << X1  SQ  X2  X3  3  ^  *  +  >>   ENTER
      <<  X1  X3  *  SQ  X2  * >>         ENTER
          << X1  X2  X3  *  + >>             COV¯     gives

    [ [   9/8      7/6   -11/24 ]
      [   1/6    25/24   11/24 ]               for instance   D₃ V¹  = 37 / 24
      [ 37/24  35/24    5/4   ] ]
 

Differential Manifolds with Non-Symmetric Connexions
 

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

>>>  So, you have to store  n³  functions in the variables 'C111'   'C112'  'C121'   and so on...
         and the n(n+1) /2  functions g_ij in 'G11'  'G22'  'G12'  ......
 

-The variables 'h'  &  'N'  are also to be initialized before executing the programs below:
-And of course,  'X1'  'X2'  .......
 

>>>  If the tensor depends on the metric tensor, execute  'INIG' first to calculate and store the coefficients of [ g_ij ]  and its inverse.
 

Initialization
 

 INIG simply calculates the metric tensor and its inverse
 
 
 

STACK	INPUT	OUTPUT
Level 1	/	"DONE"

 
Example:   A 3-dimensional manifold, with the same metric as above and the same coordinates:  x = y = z = 1

    g₁₁ = 1 + x² y² z²                            g₁₂ = x y
    g₂₂ = 1 + x² + y² + z²                     g₁₃ = x z
    g₃₃ = 1 + x² y² + x² z² + y² z²         g₂₃ = y z

-Store the same programs in the variables 'G11'  'G22'  .....   3 in 'N'  and  -0.01  in  'h'

  INIG  stores the same matrix as above

                                             [ [ 2  1  1 ]                                       [ [  5/8  -1/8   -1/8   ]
   The metric tensor  in  GIJ  =  [ 1  4  1 ]       its inverse in  G^IJ      [ -1/8  7/24  -1/24 ]
                                                [ 1  1  4 ] ]                                      [ -1/8 -1/24   7/24 ] ]
 

-Now, suppose the connexion is defined by ( I don't know if it is realistic ):

   G¹₁₁ = x y² z²                 G¹₂₂ = 0
   G¹₁₂ = x² y z²   = - G¹₂₁ G¹₂₃ = - G¹₃₂ = 1
   G¹₁₃ = x² y² z   = - G¹₃₁  G¹₃₃ = 0

   G²₁₁ =  x² y z                   G²₂₂ =  1
   G²₁₂ =  x y² z   = - G²₂₁   G²₂₃ = G²₃₂ = 0
   G²₁₃ =  x y z²   = - G²₃₁   G²₃₃ = 1

   G³₁₁ =  x y                        G³₂₂ =  0
   G³₁₂ =  x z  =  - G³₂₁        G³₂₃ = - G³₃₂ = 1
   G³₁₃ =  y z  =  - G³₃₁        G³₃₃ =  2

-Store  <<  X1  X2  X3  *  SQ  *  >>  in  'C111'  ..................................  and  2  in  'C333'   ( 27 functions )
 

Curvature Tensor (1-3)
 

  RIJK^L calculates the curvature tensor defined as

    R^l_ijk = ¶_jG^l_ik - ¶_kG^l_ij + G^l_mjG^m_ik - G^l_mkG^m_ij
 
 
 

STACK	INPUTS	OUTPUTS
Level 4	i	/
Level 3	j	/
Level 2	k	/
Level 1	l	Rlijk

 
Examples:

  2  ENTER
  1  ENTER
  2  ENTER
  1  RIJK^L    >>>>  R¹₂₁₂ = 0

-Likewise

  3  ENTER
  1  ENTER
  3  ENTER
  2  RIJK^L  >>>>  R²₃₁₃ = 1
 

Ricci Tensor
 

-Ricci tensor is still  R_ij = R^m_imj
 
 
 

STACK	INPUTS	OUTPUTS
Level 2	i	/
Level 1	j	Rij

 
Examples:

  1  ENTER
  1    RIJ     >>>>  R₁₁ = 8

-Likewise

  2  ENTER
  3    RIJ      >>>>  R₂₃ =  5
 

Segmental Curvature
 

-There is also a segmental curvature tensor  Q_ij = R^m_mij   which equals zero in Riemannian manifolds.

      Q_ij = ¶_iG^m_mj - ¶_jG^m_mi  is actually a curl.
 
 

STACK	INPUTS	OUTPUTS
Level 2	i	/
Level 1	j	Qij

 
Examples:

  1  ENTER
  2    QIJ     >>>>  Q₁₂ = 3

-Likewise

  1  ENTER
  3    QIJ      >>>>  Q₁₃ =  2
 

Scalar Curvatures
 

-The contracted tensors  R = g^ij R_ij  and  Q = g^ij Q_ij  gives 2 scalars.

>>> To obtain correct results, "INIG" must be executed before RRR & QQQ  to store  g^ij  in the proper registers.
 
 

STACK	INPUT	OUTPUT
Level 1	/	R or Q

 
Examples:   With the same metric & connexion and at the same point, we get:

   RRR  >>>   R = 181/24

   QQQ >>>  Q = 0
 

Torsion Tensor
 

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

     S^k_ij = G^k_ij -  G^k_ji                ( some authors define  S^k_ij = (1/2) ( G^k_ij -  G^k_ji )    others by   S^k_ij = ( G^k_ji -  G^k_ij )  )

-Though G^k_ij is not a tensor, S^k_ij is a tensor !
 
 

STACK	INPUTS	OUTPUTS
Level 3	i	/
Level 2	j	/
Level 1	k	Skij

 
Examples:

   1  ENTER
   2  ENTER
   3   SIJK    >>>>   S³₁₂ = +2

   3  ENTER
   2  ENTER
   1   SIJK      >>>>  S¹₃₂ = -2
 

Torsion Vector
 

-Contracting  S^k_ij  we get a covariant vector  S_j = S^k_jk
 
 

STACK	INPUTS	OUTPUTS
Level 1	i	Si

 
Examples:

   1   SI    >>>>    S₁ =  4
   2   SI    >>>>    S₂ =  0
   3   SI    >>>>    S₃ = -2
 
 

References:

 [1]  https://en.wikipedia.org/w/index.php?title=Differential_geometry_of_curves&oldid=740874421
 [2]  https://en.wikipedia.org/w/index.php?title=Frenet%E2%80%93Serret_formulas&oldid=731811089
 [3]  Ron Goldman - "Curvature formulas for implicit curves and surfaces"
 [4]  Elie Cartan - "Leçons sur la géométrie des espaces de Riemann" - Gauthier-Villars, Paris   ( in French )
 [5]  Denis-Papin & Kaufmann - "Cours de calcul Tensoriel Appliqué" - Albin Michel    ( in French )
 [6]  List of Formulas in Riemannian Geometry
 [7]  Nikodem J. Poplawski - "Spacetime and fields"- Department of Physics, Indiana University, Bloomington, IN 47405, USA
 [8]  Marie-Antoinette Tonnelat - Theorie unitaire affine du champ physique. J. Phys. Radium, 1951, 12 (2),
        pp.81-88. <10.1051/jphysrad:0195100120208100>. <jpa-00234360>
 [9]  Gerald Kaiser - "Quantum Physics, Relativity, and Complex Spacetime: Towards a New Synthesis"
[10] Giampiero Esposito - "From Spinor Geometry to Complex General Relativity"
[11] Jean E. Charon - "Complex General Relativity"

-The last reference is controversial, but there are many interesting ideas in Charon's theory...