DERIVITGR Module


 Overview
 

-Many programs in this module are also in DERIVE41 module.

-But several have been improved, especially CDGL+
  which also gives now the vectorial Laplacian in cylindrical & spherical coordinates.

-2 mixed derivatives are computed by new formulae which require less function-evaluations
  without beeing less accurate ( perhaps even with a better precision ).

-There are also a few routines to evaluate integrals with Gaussian formulas and to solve differential equations.
 
 

Notes:

-Most of the formulas for numerical differentiation are of order 10.
-This produces a good precision for the 1st derivatives but a low precision for 6th derivative.
 
 
 
XROM  Function  Desciption
 31,00
 31,01
 31,02
 31,03
 31,04
 31,05
 31,06
 31,07
 31,08
 31,09
 31,10
 31,11
 31,12
 31,13
 31,14
 31,15
 31,16
 31,17
 31,18
 31,19
 31,20
 31,21
 31,22
 31,23
 31,24
 31,25
 31,26
 31,27
 31,28
 31,29
 31,30
 31,31
 31,32
 31,33
 31,34
 31,35
 31,36
 31,37
 31,38
-DERIVITGR
 "DET"
 "dF"
 "d4F"
 "dF2"
 "dF3"
 "MDV"
 "MDV12"
 "MDV21"
 "MDV3"
 "MDV4"
 "FD"
 "SD"
 "TD"
 "BHRM"
 "CDGL+"
 "DAL"
 "HEAT"
 "LAPN"
 "THRM"
 "HESSGR"
 "CRVL"
 "KGM"
 "KGMN"
 "KHT"
 "KHTN"
 "RK7"
 "2RK7"
 "RKN6"
 "2RKN6"
 "GRK"
 "GL16"
 "DTI"
 "MIT"
 "MIT1"
 "ISSPH"
 NORM
 X^3
 X/3
 Section Header
 Determinant of a Square Matrix of order n ( n < 18 )
 1st 2nd & 3rd Derivatives of a Function of 1 var.
 4th 5th 6th Derivatives of a Function of 1 variable
 Derivatives of a Function of 2 variables
 Derivatives of a Function of 3 variables
 Mixed Derivative of a Function of 2 vatiables
 d3F/dxdy2
 d3F/dx2dy
 d3F/dxdydz
 d4F/dx2dy2
 First Derivatives of a Function of N variables
 Second Derivatives of a Function of N variables
 Third Derivatives of a Function of N variables
 Biharmonic Operator ( 3 Dim. )
 Curl, Divergence, Gradient, Laplacian, Vector-Lapl.
 D'Alembertian Operator
 Heat Operator
 Laplacian of a Function of N variables
 Triharmonic Operator ( 3 Dim. )
 Hessian Matrix + Gradient
 Arc Length of a Parametric 3D-Curve
 Gaussian Curvature & Mean Curvature of a Surface.
 Gaussian Curvature & Mean Curv. of a Hypersurface
 Curvature & Torsion of a Curve ( 3 Dim )
 Curvature & Torsion of a Curve ( n Dim )
 Differential Equation with Runge-Kutta of order 7
 Systems of 2 ODE with Runge-Kutta of order 7
 Conservative ODE, Runge-Kutta-Nystrom of order 6
 Systems of 2 conservative ODE with RKN6
 Systems of n ODE with Gill-Runge-Kutta method
 Gauss-Legendre Integration ( 16-point Formula )
 Double & Tripple Integrals with 3-point Gauss-Leg.
 Multiple Integrals with GL3
 Multiple Integrals from -1 to +1 with GL3
 Integral on the Surface of a Sphere
 Norm of a 3D-vector
 An M-Code Routine to Calculate X^3
 An M-Code Routine to Calculate X/3

 

Determinant of a Square Matrix
 

-"DET" is very similar to "LS1" ( cf "Linear and non-linear Systems for the HP-41" )
-You can store the coefficients of the matrix from Rbb to Ree provided  bb > 00
 
 

Data Registers:         R00 = det A at the end                   ( Registers Rbb thru Rbb-1+n2 are to be initialized before executing "DET" )

                                  Rbb   = a11  ....................     Ree   = ann

Flags: /
Subroutines: /
 
 
      STACK        INPUT      OUTPUT
           X     bbb.eeenn     Deteminant
           L             /     bbb.eeenn

   where n is the order of the square matrix and  bbb > 000

Example:    Calculate the following determinant assuming the coefficients are stored in R11 to R19

            |  4  9  2  |                  R11  R14  R17
    D =  |  3  5  7  |      into      R12  R15  R18      respectively
            |  8  1  6  |                  R13  R16  R19

    11.01903   XEQ "DET"  >>>>   Det = 360    ( in 10 seconds )

Notes:

-This program does not check that the control number is valid.
-"DET" is called by "KGMN"
 

1st 2nd & 3rd 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 = 13s---
                             RDN   f "(2) = 7.139056635 = R04
                             RDN  f '''(2) = 7.639051753 = R05
 

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

-"d4F"  calculates the 4th, 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:

   h4 f(4)(x) ~  ( 192654 F - 140196 G + 52428 H - 9738 I +1261 J - 82 K ) / 15120
   h5 f(5)(x) ~  ( 1938 A - 1872 B + 783 C - 152 D + 13 E ) / 288
   h6 f(6)(x) ~  ( -233244 F + 184110 G - 88920 H + 24795 I -3610 J + 247 K ) / 4560
 

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

                                         R01 = x     R03 = f(x)         R05 = d5f/dx5    R07 to R09: temp
                                         R02 = h     R04 = d4f/dx4   R06 = d6f/dx6

Flags:  /
Subroutine:      A program which computes f(x) assuming x is in X-register upon entry
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             /         f(6)(x)
           Y             h         f(5)(x)
           X             x         f(4)(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 "d4F"  >>>>   f '(2) = 7.013760979 = R04                                                     ---Execution time = 14s---
                               RDN   f "(2) = 8.140451389 = R05
                               RDN  f '''(2) = 5.637850877 = R06
 

Derivatives of a Function of 2 variables
 

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

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

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

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

Flags: /
Subroutine:      A program which computes f(x,y) assuming x is in X-register & y is in Y-register upon entry
 
 
 
      STACK        INPUTS      OUTPUTS
           T             /   f "yy = 2f / y2
           Z            h   f "xx = 2f / x2
           Y             y    f 'y = f / y
           X             x    f 'x = f / x

 
Example:     f(x) = exp(-x2).Ln(y)    x = 1 , y = 2
 
 
 01  LBL "T"
 02  X^2
 03  CHS
 04  E^X
 05  X<>Y
 06  LN
 07  *
 08  RTN 

 
    T  ASTO 00

-If we choose h = 0.1

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

-You have also:

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

Derivatives of a Function of 3 variables
 

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

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

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

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

Flags: /
Subroutines:   "dF" and a program which computes f(x,y,z) assuming x is in X-register, y is in Y-register and z is in Z-register upon entry.
 
 
 
      STACK        INPUTS      OUTPUTS
           T             h            Df
           Z             z    f 'z = f / z
           Y             y    f 'y = f / y
           X             x    f 'x = f / x

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

Example:     f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 2 , z = 3
 
 
 01  LBL "T"
 02  X^2
 03  CHS
 04  E^X
 05  RDN
 06  X^2
 07  +
 08  LN
 09  R^
 10  *
 11  RTN 

 
      T  ASTO 00

-If we choose h = 0.1

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

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

Mixed Derivative of a Function of 2 variables
 

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

     f "xy = 2f / xy = (1/h2) SUMj=1,2,3,4,5  Aj ( fjj + f-j-j - fj-j - f-jj )

   where  A1 = 5/12     A2 = -5/84     A3 = 5/504     A4 = -5/4032     A5 = 1/12600
 

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

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

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

 
Example:     f(x) = exp(-x2).Ln(y)    x = 1 , y = 2
 
 
  01  LBL "T"
  02  X^2
  03  CHS
  04  E^X
  05  X<>Y
  06  LN
  07  *
  08  RTN 

 
    T  ASTO 00

-If we choose h = 0.1

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

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

Notes:

-The function is evaluated 20 times instead of 31 with the previous version.
-So, the program is faster and even seems a little more accurate.
 

Mixed Derivative d3F/dx2dy & d3F/dxdy2
 

 "MDV21" & "MDV12" employ a method of order 11 to estimate  f'''xxy = 3f / x2y  and  f'''xyy = 3f / y2x

    h3 f'''xxy ~  SUMk=1,2,...,5  Ak [ - 2 f0k + 2 f0-k  + ( fkk - f-k-k + f-kk - fk-k ) ]

 With    A1 = 5/6 ,  A2 = -5/84  ,  A3 = +5/756  ,  A4 = -5/8064  ,  A5 = +1/31500
 

Data Registers:             R00 = f name                                 ( Register R00 is to be initialized before executing "MDV21" or "MDV12 )

                                         R01 = x     R03 = h
                                         R02 = y     R04 = f'''xxy = 3f / x2y  or  f'''xyy = 3f / y2x     R05-R06  temp

Flag:  F02            CF 02 ->  f'''xxy = 3f / x2y
                             SF 02  ->  f'''xyy = 3f / y2x

Subroutine:  A program that takes x , y  in registers X , Y respectively and returns  f(x,y)  in X-register.
 
 
      STACK        INPUTS      OUTPUTS
           T             h             /
           Z             z             /
           Y             y             /
           X             x    f'''xxy or f'''xyy

     CF 02 -> f'''xxy
     SF 02 -> f'''xyy

Example:       f(x,y,z) = Ln ( 1 + x2 y )    x = 1  ,  y = 2
 
 
 01  LBL "T"
 02  X^2
 03  *
 04  LN1+X
 05  RTN

 
    "T"  ASTO 00    and if you choose  h = 0.1

    0.1   ENTER^
     2     ENTER^
     1     XEQ "MDV21"  >>>>   f'''xxy = 3f / x2y  =  - 0.370369497                ---Execution time = 25s---

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

-With the same inputs,  XEQ "MDV12"  returns   f'''xyy = 3f / y2x  =  - 0.148148268

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

Mixed Derivative d3F/dxdydz
 

 "MDV3" employs a new method of order 11 to estimate  f'''xyz = 3f / xyz

    h3 f'''xyz ~  SUMk=1,2,...,5  Ak [  fkkk - f-k-k-k - fkk-k + f-k-kk - fk-kk + f-kk-k - f-kkk + fk-k-k ) ]

 With    A1 = 5/24 ,  A2 = -5/336  ,  A3 = +5/3024  ,  A4 = -5/32256  ,  A5 = +1/126000
 

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

                                         R01 = x     R04 = h
                                         R02 = y     R05 = f'''xyz = 3f / xyz
                                         R03 = z     R06-R07:  temp
Flags: /
Subroutine:  A program that takes x , y , z in registers X , Y , Z respectively and returns  f(x,y,z)  in X-register.
 
 
      STACK        INPUTS      OUTPUTS
           T             h             /
           Z             z             /
           Y             y             /
           X             x          f'''xyz

 
Example:     f(x,y,z) = Ln ( 1 + x2 + 2.y + z3 )    x = y = z = 1
 
 
 01  LBL "T"
 02  X^2
 03  X<>Y
 04  ST+ X
 05  +
 06  X<>Y
 07  ENTER^
 08  X^2
 09  *
 10  +
 11  LN1+X
 12  RTN

 
    "T"  ASTO 00    and if you choose  h = 0.1

    0.1   ENTER^
     1     ENTER^  ENTER^   XEQ "MDV3"  >>>>   f'''xyz = 3f / xyz  =  0.192000090                ---Execution time = 39s---

-The exact value is  0.192

Notes:

-With this new version, the function is evaluated 40 times instead of 70.
-The program becomes faster ( 39s instead of 62s in the example above ) and the precision better.
 

Mixed Derivative d4F/dx2dy2
 

 "MDV4" evaluates   f'''xxyy = 4f / x2y2  with a formula of order 12

  h4 f'''xxyy ~  20881861 f00 / 3240000  + SUMk=1,....,5  Ak [ fkk + f-k-k + fk-k + f-kk - 2 ( fk0 + f-k0 + f0k + f0-k ) ]

  where   A1 = 5/3 ,  A2 = -5/84  ,  A3 = +5/1134  ,  A4 = -5/16128  ,  A5 = +1/78750
 

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

                                         R01 = x     R03 = h
                                         R02 = y     R04 = f""xxyy = 4f / x2y2     R05-R06-R07:  temp

Flags:  /
Subroutine:  A program that takes x , y  in registers X , Y respectively and returns  f(x,y)  in X-register.
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             h             /
           Y             y             /
           X             x         f(4)xxyy

 
Example:    f(x,y) = Ln ( x2 + y3 )     x = 2 , y = 1
 
 
 01  LBL "T"
 02  X^2
 03  X<>Y
 04  ENTER^
 05  X^2
 06  *
 07  +
 08  LN
 09  RTN

 
    "T"  ASTO 00    and if you choose  h = 0.1

    0.1  ENTER^
     1    ENTER^
     2    XEQ "MDV4"   >>>>  f(4)xxyy4f / x2y2 =  - 0.038395989                                             ---Execution time = 33s---

-Exact result = - 0.0384
 

First Derivatives of a Function of N variables ( N < 11 )
 

 "FD" employs the same formula as above
 

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

                                R01 = x1 ,    R02 = x2 , .......... ,    Rnn = xn      n < 11

                                           R10 = h          R11 temp         R13 = i        R15 = xi
                                           R12 =  f / xi       R14 = 0 , h , 2h , ....     R16: temp

Flags:  /
Subroutine:   A program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn
 
 
      STACK        INPUTS      OUTPUTS
           Y             h             /
           X             i    f 'xi = f / xi

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

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

-With  h = 0.1

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

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

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

Second Derivatives of a Function of N variables ( N < 11 )
 
 

 "SD" employs the same formula as above
 
 

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

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

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

 
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.509989206   ( exact derivative = 0.509989195 )                   ---Execution time = 12s---

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

Third Derivatives of a Function of N variables ( N < 11 )
 

 "TD" employs the same formula as above
 
 

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

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

                                           R10 = h              R12-R18:  temp
                                           R11 =  3f / xixjxk
Flags:  F09-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 / xixj

 
Example:             f(x,y,z,t) = exp(-x2) Ln( y2 + z + t3 )     x = 1 , y = 2 , z = 3 , t = 1
 
 
 01  LBL "T"
 02  RCL 01
 03  X^2
 04  CHS
 05  E^X
 06  RCL 02
 07  X^2
 08  RCL 03
 09  +
 10  RCL 04
 11  ENTER^
 12  X^2
 13  *
 14  +
 15  LN
 16  *
 17  RTN
 18  END

 
    "T"  ASTO 00

     1   STO 01
     2   STO 02
     3   STO 03
     1   STO 04

    f'''xyz = 3f / xyz  = ?

-With   h = 0.1

     0.1  ENTER^
      1    ENTER^
      2    ENTER^
      3    XEQ "TD"   >>>>    f'''xyz = 3f / xyz  = 0.045984935                                          ---Execution time = 57s intead of 88s---

-Exact result = 0.045984930....

    f'''xtt = 3f / xt2  = ?

-With   h = 0.1

     0.1  ENTER^
      1    ENTER^
      4    ENTER^    XEQ "TD"   >>>>    f'''xtt = 3f / xt2  = - 0.448353739                                    ---Execution time = 42s---

-Exact result = - 0.448353069....

    f'''yyy = 3f / y3  = ?

-With   h = 0.1

     0.1  ENTER^
      2    ENTER^   ENTER^    XEQ "TD"   >>>>    f'''yyy = 3f / y3  = - 0.045984326                               ---Execution time = 15s---

-Exact result = - 0.0459849301...
 

Hessian Matrix + Gradient
 

 "HESSGR" calculates and sores the coefficients of the Hessian matrix H and the coefficients of the gradient G into registers R26  R27 ..... in the following way:

       W =   H  Gt      where  Gt  = transposed G              So, we get a square matrix W of order (n+1)
                 G  0

  with        H = [ hi,j ]   with  hi,j = 2f / xixj
                G = [ gk ]    with  gk = f / xi

-This form is useful to calculate the Gaussian curvature of a hypersurface.
 
 

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

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

                                 R11 = h       R12 to R25 :  temp    R26 ...... = W

Flag:  F10
Subroutines:   "FD"  "SD"  and a program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn
 
 
      STACK        INPUTS      OUTPUTS
           Y             h             /
           X          n < 11      26.eeemm

   Where  m = n+1

Example:    f(x,y,z) = x4 y3 z2 - 1    x = y = z = 1
 
 
 01  LBL "T"
 02  RCL 01
 03  X^2
 04  RCL 03
 05  *
 06  X^2
 07  RCL 02
 08  ST* Y
 09  X^2
 10  *
 11  1
 12  -
 13  RTN

 
  "T"  ASTO 00

   1  STO 01  STO 02  STO 03   and  if you choose  h = 0.1

   0.1   ENTER^
    3     XEQ "HESSGR"  >>>>    26.04104                                    ---Execution time = 111s---

-And we find the matrix W in registers  R26 to R41

      12   12   8   4
      12    6    6   3
       8     6    2   2
       4     3    2   0

                                    12  12  8
-The Hessian matrix is   12   6   6  and the gradient is  [ 4  3  2 ]
                                     8    6   2
 

Laplacian of a Function of N variables ( N < 11 )
 

 "LAPN" may be used to evaluate the Laplacian of a function of n variables, provided n < 11
 
 
 
      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 "LAPN"  >>>>    Laplacian(f) =  2f / x2 + 2f / y2 + 2f / z2 + 2f / t2 = 0.673013929                  ---Execution time = 54s---

 ( the last 2 decimals should be 32 )
 

D'Alembertian Operator
 

-The d'Alembertian of a function of 4 variables f(x,y,z,t) is defined as:

   (1/c2) 2f / t2 - ( 2f / x2 + 2f / y2 + 2f / z2 )       where  c = speed of light

-Here, we assume the units are choosen so that c = 1
 
 
      STACK        INPUT      OUTPUT
           X             h  d'Alembertian(f)

 
Example:    f(x,y,z,t) = exp(-t) ( x2 + y2 + z2 )       x = y = z = t = 1
 
 
 01 LBL "T"
 02 RCL 01
 03 X^2
 04 RCL 02
 05 X^2
 06 RCL 03
 07 X^2
 08 +
 09 +
 10 RCL 04
 11 CHS
 12 E^X
 13 *
 14 RTN

 
   T  ASTO 00

   1  STO 01  STO 02  STO 03  STO 04

  0.1  XEQ "DAL"  >>>>   d'Alembertian(f) = -1.103638055                                   ---Execution time = 46s---

Note:

-Several authors define the d'Alembertian as  2f / x2 + 2f / y2 + 2f / z2 - (1/c2) 2f / t2
-Simply change the sign of the result.
 

Heat Operator
 

-It is defined by   f / t - D ( 2f / x2 + 2f / y2 + 2f / z2 )      where D is a positive constant
 
 
      STACK        INPUTS      OUTPUTS
           Y             h             /
           X             D         Heat(f)

 
Example:      f(x,y,z,t) = exp(-t) ( x2 + y2 + z2 )       x = y = z = t = 1       D = 0.7
 
 
 01 LBL "T"
 02 RCL 01
 03 X^2
 04 RCL 02
 05 X^2
 06 RCL 03
 07 X^2
 08 +
 09 +
 10 RCL 04
 11 CHS
 12 E^X
 13 *
 14 RTN

 
   T  ASTO 00

   1  STO 01  STO 02  STO 03  STO 04

  0.1  ENTER^
  0.7  XEQ "HEAT"  >>>>   Heat(f) = -2.648731654                                   ---Execution time = 44s---
 

Biharmonic Operator ( 3-Dim )
 

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

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

-"BHRM3" uses 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 "BHRM3" )

                                         R01 = x        R04 = h
                                         R02 = y        R05 = D2 f
                                         R03 = z        R06-R07-R08-R09-R10:  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.
 
 
      STACK        INPUTS      OUTPUTS
           T            h             /
           Z             z             /
           Y             y             /
           X             x          D2 f

 
Example:     f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 2 , z = 3
 
 
 01  LBL "T"
 02  X^2
 03  CHS
 04  E^X
 05  RDN
 06  X^2
 07  +
 08  LN
 09  R^
 10  *
 11  RTN 

 
      T  ASTO 00

-If you choose h = 0.1

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

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

-With h = 0.2  it yields   -14.34286490 and with h = 0.15  we get  -14.34270102

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 5 or 6 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
 

   D3 f = 6f / x6 + 6f / y6 + 6f / z6 + 3 ( 6f / x4y2 + 6f / x4z2 + 6f / y4z2 + 6f / x2y4 + 6f / x2z4 + 6f / y2z4 ) + 6 6f / x2y2z2
 

   h6 D3 f  ~  SUMm=1....5  Am ( fm00 + f-m00 + f0m0 + f0-m0 + f00m + f00-m - 6 f000 )
                                     + Bm ( fmm0 + f-m-m0 + fm-m0 + f-mm0 + fm0m + f-m0-m + fm0-m + f-m0m + f0mm + f0-m-m + f0m-m + f0-mm - 12 f000 )
                                     + Cm ( fmmm + f-m-m-m + fmm-m + f-m-mm + fm-mm + f-mm-m + fm-m-m + f-mmm - 8 f000 )

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

     A1 = +22997 / 120                  B1 = - 2869 / 60                C1 = + 10
     A2 = - 12709 / 420                  B2 = + 667 / 240               C2 = - 5 / 56
     A3 = + 858391 / 136080         B3 = - 15007 / 68040        C3 = + 5 / 1701
     A4 = - 137843 / 161280          B4 = + 5119 / 322560       C4 = - 5 / 43008
     A5 = + 894317 / 15750000     B5 = - 739 / 1125000        C5 = + 1 / 328125
 
 

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

                                         R01 = x        R04 = h
                                         R02 = y        R05 = D3 f
                                         R03 = z        R06 to R12:  temp   ( R07 = 6 f(x,y,z) )
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.
 
 
 
      STACK        INPUTS      OUTPUTS
           T          h > 0             /
           Z             z             /
           Y             y             /
           X             x     D3 f(x,y,z)

 
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 "THRM"  >>>> D3 f (1,2,3) = 133.675100 = R05                                                 ---Execution time = 2m42s---

-With a similar program, the HP48 - which works with 12 digits - gives  133.533179
 

Notes:

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

Arc Length of a Parametric Curve ( 3-Dim )
 

 "CRVL" evaluates the integral  §ab  [ ( dX/dt )2 + ( dY/dt )2 + ( dZ/dt )2 ] 1/2  dt

-Romberg method is used with Pythagora's theorem.
 
 

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

                                        R01 = a      R03 = L      R04 to R09: temp    R20, R21, .... are used by "ROM"
                                        R02 = b
Flags: /
Subroutines:  3 programs LBL "X"  LBL "Y"  LBL "Z"  that take t in X-register and return  X(t) , Y(t) , Z(t)  in register X  respectively
 
 
      STACK        INPUTS      OUTPUTS
           Y             a             /
           X             b        L(a,b)

 
Example:    X(t) = t4 ,  Y(t) = t2  ,  Z(t) = exp t   ;  a = 0 ,  b = 1
 
 
 01  LBL "X"
 02  X^2
 03  X^2
 04  RTN 
 05  LBL "Y"
 06  X^2
 07  RTN
 08  LBL "Z"
 09  E^X"
 10  RTN

 
     0     ENTER^
     1     XEQ "CRVL"  >>>>   L = 2.342116456                        ---Execution time = 2m49s---

Notes:

-The HP41 displays the successive approximations
-Exact result  = 2.342116459....
 

Gaussian Curvature & Mean Curvature of a Surface
 

-"KGM" calculates the Gaussian curvature KG and the mean curvature Km of a surface defined by a function  z = z(x,y)

-The Gaussian Curvature KG is be computed by:

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

-The mean curvature Km is given by:

    Km = (1/2) [ t ( 1+p2 ) + r ( 1+q2 ) - 2 p q s ] / ( 1+p2+q2 )3/2

   where  p = z / x , q = z / y,  r = 2z / x2,  s = 2z / xy,  t = 2z / y2
 
 

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

                                         R01 to R15:  temp
Flags: /
Subroutines:   "dF2" and a program that takes x in X-register, y in Y-register and returns z(x,y) in X-register
 
 
      STACK        INPUTS      OUTPUTS
           Z             h             /
           Y             y            Km
           X             x            KG

 
Example:   A surface is defined by  z(x,y) = Ln ( 1 + x2 y3 )    Calculate the Gaussian & mean curvatures at the point  (x,y) = (1,1)
 
 
 01  LBL "T"
 02  X^2
 03  X<>Y
 04  ST* Y
 05  X^2
 06  *
 07  LN1+X
 08  RTN

 
  "T"  ASTO 00   and if you choose h = 0.1

     0.1  ENTER^
      1    ENTER^    XEQ "KGM"  >>>>   KG = -0.124570010                             ---Execution time = 49s instead of 58s---
                                                     X<>Y  Km = -0.171204283

Notes:

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

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

Gaussian Curvature & Mean Curvature of a HyperSurface
 

 -Here, we assume that the hypersurface is defined implicitly by   f( x1 , ........... , xn ) = 0            ( n < 11 )

Formulae:

   KG = - ( det W ) / ( - | Grad f | )n+1      where  W is defined in paragraph 1°) c)

   Km = [ Grad f x H x TGrad f  -  | Grad f |2 Trace (H) ] / (n-1) / | Grad f |3     where  H = Hessian matrix
 
 

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

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

                                 R11 = h       R12 to R25 :  temp    R26 ...... = W

Flag:  F10
Subroutines:  "HESSGR" and a program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn
 
 
      STACK        INPUTS      OUTPUTS
           Z             /        det W
           Y             h            Km
           X         n < 11            KG

 
Example:   The surface ( n = 3 )   x4 y3 z2 - 1 = 0   with  x = y = z = 1
 
 
 01  LBL "T"
 02  RCL 01
 03  X^2
 04  RCL 03
 05  *
 06  X^2
 07  RCL 02
 08  ST* Y
 09  X^2
 10  *
 11  1
 12  -
 13  RTN

 
  "T"  ASTO 00

   1  STO 01  STO 02  STO 03   and  if you choose  h = 0.1

   0.1   ENTER^
    3     XEQ "KGMN"  >>>>      KG  =  0.256837099                             ---Execution time = 137s---
                                     RDN      Km  =  0.518666290
                                    X<>Y  det W = - 216

Warning:

-This program does not check that  f(x1, ..... ,xn) = 0
-The matrix W is not saved.
 

Curvature & Torsion of a Parametric Curve ( 3 Dimensions )
 

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

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

Data Registers:              R00 & R05 to R15:  temp
                                         R01 = t  R02 = h   R03 = khi   R04 = tau
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
 
 
      STACK        INPUTS      OUTPUTS
           Y             h          tau(t)
           X             t          khi(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 khi (t) and the torsion tau (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"  >>>>  khi (1.3)  = 0.194807823                       ---Execution time = 46s---
                                 X<>Y  tau (1.3)  = 0.123665258

-The exact values are   khi (1.3) = 0.194807823  &  tau (1.3) = 0.123665529
 

Curvature & Torsion of a Parametric Curve ( N Dimensions )
 

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

Data Registers:              R00 & R05 to R09 & R11 to R5n+12:  temp

                                         R01 = t    R03 = khi
                                         R02 = h   R04 = tau
Flags: /
Subroutines:    "dF" and  n functions named  "X1"  "X2" ........... "XN"   that compute x1(t) , x2(t) , ........ , xn(t)  assuming t is in X-register upon entry
 
 
 
      STACK        INPUTS      OUTPUTS
           Z          n < 62             /
           Y             h          tau(t)
           X             t          khi(t)

 
Example:   The curve defined by

     x1(t) = t4       x3(t) = t3      at the point  t = 1
     x2(t) = t2       x4(t) = et
 
 
 01  LBL "X1"
 02  X^2
 03  X^2
 04  RTN
 05  LBL "X2"
 06  X^2
 07  RTN
 08  LBL "X3"
 09  ENTER^
 10  X^2
 11  *
 12  RTN
 13  LBL "X4"
 14  E^X
 15  RTN

 
-If you choose  h = 0.1

     4   ENTER^
   0.1  ENTER^
     1   XEQ "KHTN"   >>>>   khi (1)  = 0.142277239                       ---Execution time = 71s---
                                   X<>Y  tau (1)  = 0.121469324

Note:

-If n = 3, the sign of tau may be different from the result given by "KHT" because Gram-Schmidt process doesn't always give a direct basis.
-So use KHT in this case, which is also faster.
 

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.
-Flag F01 is cleared if flag F02 is set ( lines 02-03 ).
 

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

-The coordinates of the vector Laplacian are simply the Laplacians of the coordinates
  ( It's not so simple in cylindrical & spherical coordinates... )
 

      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

-For the vector-laplacian, the formulas are given  here
 

      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

-For the vector-laplacian, the formula, even more complicated than the cylindrical case, is given here
 
 

Data Registers:              R00 = h        ( at the end )

                                        R01 = x         R04-R05-R06 = Curl E       R08-R09-R10 = Grad (X)    R17 = Lap (X)    R20 to R22 = Vector Laplacian
                                        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.
 
 
 
      STACK       INPUTS      OUTPUTS
           T            h         Div E
           Z            x3        Curl3 E
           Y            x2        Curl2 E
           X            x1        Curl1 E
           L             /         4.022

    4.022 = 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 = 90s---
                                   RDN    -32.56682780  = R05
                                   RDN     71.78978318  = R06
                                   RDN     45.44140669  = R07

-So,  Curl E = ( 8.619381890 ,  -32.56682780 ,  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.61938200 , 32.61938189 , 10.87312737 )

    Lapl f = 1.409197689
    Lapl g = 98                               which are also in R20 thru R22
    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.351141010  = R04                                                                     ---Execution time = 94s---
                     RDN    -8.326237923  = R05
                     RDN     5.048943482  = R06
                     RDN     7.163118954  = R07

-So,  Curl E = ( -6.351141010 ,  -8.326237923 ,  5.048943482 )

 and    Div E = 7.163118954

-You also find in registers R08 to R19:

    Grad f = ( 0.345491503 , 0.951056518 , 1.381966010 )
    Grad g = ( 4 , 0 , 4 )
    Grad h = ( 9.708203933 , -2.351141010 , 6.472135948 )

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

-And the coordinates of the vector-laplacian in R20 to R22

  VLap =  ( 1.690982985 , 3.951056518 , 12.94427209 )

-In cylindrical coordinates, the 3rd coordinate of the vector-laplacian = the Laplacian of the 3rd coordinate.
-This is not true in spherical coordinates:
 

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.379867348  = R04                                                                     ---Execution time = 145s---
                     RDN    -6.059707088  = R05
                     RDN     2.959890530  = R06
                     RDN    -0.045010876  = R08

-So,  Curl E = ( -3.379867348 ,  -6.059707088 ,  2.959890530 )

 and    Div E = -0.045010876

-You also find in registers R08 to R19:

    Grad f = ( 0.490881375 , 0.566820985 , -0.823639103 )
    Grad g = ( 2.351141010 , 0 , 1.868344718 )
    Grad h = ( 3.927050990 , -2.267283945 , -2.196370945 )

    Lapl f = 0.018237234
    Lapl g = 2.742997604
    Lapl h = 5.721039300

-And the coordinates of the vector-laplacian in R20 to R22

  VLap =  ( 1.045010859 , 3.794180276 , 5.103411526 )
 

Differential Equation with a Runge-Kutta Formula of order 7
 

  "RK7" solves the differential equation

       y' = f(x,y)   with initial value   y(x0) = y0
 

-This method requires 9 evaluations of the function per step
 
 

Data Registers:         R00 =  f name                                                ( Registers R00 to R04 are to be initialized before executing "RK7" )

                                    R01 = x0         R03 =  h = the stepsize                      R05  to R12: temp
                                    R02 = y0         R04 = N = the number of steps         R20 to R54: Coefficients

Flags: /
Subroutine:  a program calculating f(x;y) assuming x is in X-register and y is in Y-register upon entry.
 
 
 
      STACK        INPUTS      OUTPUTS
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:     y(0) = 1 and  dy/dx = -2 x.y     Evaluate  y(1)
 
 
 01  LBL "T"
 02  *
 03  ST+ X
 04  CHS
 05  RTN

 
  "T" ASTO 00

   0    STO 01
   1    STO 02    and if we choose  h = 0.1
 0.1   STO 03
 10    STO 04    XEQ "RK7"   >>>>        x  =  1                                                       ---Execution time = 97s---
                              X<>Y                     y(1) = 0.3678794415

-The exact result is 0.3678794412...

Notes:

-Press R/S to continue with the same h and N.
-The next step will be faster because the constants do not have to be stored again ( 86s )
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 27 = 128 )
 

System of 2 Differential Equations with a Runge-Kutta Formula of order 7
 

  "2RK7" solves the system

       y' = f(x,y,z)     with initial values     y(x0) = y0
       z' = g(x,y,z)                                   z(x0) = z0
 

Data Registers:        R00 =  f name                                           ( Registers R00 to R05 are to be initialized before executing "2RK7" )

                                   R01 = x0         R04 = h = stepsize                       R06 to R19: temp    R55: counter
                                   R02 = y0         R05 = N = number of steps          R20 to R54: Coefficients
                                   R03 = z0
Flags: /
Subroutine:  a program calculating f(x;y;z) in Y-register and g(x;y;z) in X-register
                       assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Z             /       z(x0+N.h) 
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:     y(0) = 1 ;  z(0) = 0  ;  dy/dx = z , dz/dx = -2 x.z -2 y   ;    calculate  y(1) and z(1)
 
 
 01  LBL "U"
 02  RCL Z
 03  *
 04  +
 05  ST+ X
 06  CHS
 07  RTN

 
"U" ASTO 00

 0    STO 01
 1    STO 02
 0    STO 03   and if we choose  h = 0.1
0.1  STO 04
10   STO 05    XEQ "2RK7"   >>>>        x = 1                                                       ---Execution time = 157s---
                              RDN                      y(1) =  0.3678794408
                              RDN                      z(1) = -0.7357588818

-Exact results are  y = 0.3678794412...  &  z = -0.7357588823...

Notes:

-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when h is divided by 2 , errors are approximately divided by 27 = 128 )
 

Differential Equation with a Runge-Kutta-Nystrom Formula of order 6
 

 "RKN6" solves the conservative differential equation

     y" = f(x,y)     with the initial values:     y(x0) = y0     y'(x0) = y'0

-Five stages are sufficient to get a 6-th order formula.
 
 

Data Registers:       R00 = f name                        ( Registers R00 thru R05  are to be initialized before executing "RKN6" )

                                  R01 = x0         R04 = h  = stepsize                   R06 to R10: temp
                                  R02 = y0         R05 = N = number of steps
                                  R03 = y'0
Flags: /
Subroutine:  a program which calculates f(x,y) in X-register, assuming  x and y  are in registers X and Y upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Z             /       y'(x0+N.h) 
           Y             /       y(x0+N.h) 
           X             /          x0+N.h

 
Example:     y(0) = 1  ,   y'(0) = 1  ,   y" = - y ( x2 + y2 ) 1/2

 >>> Calculate  y(1) & y'(1)

-Load this short routine:
 
 
 01  LBL "T"
 02  X^2
 03  RCL Y
 04  X^2
 05  +
 06  SQRT
 07  *
 08  CHS
 09  RTN

 
  "T"  ASTO 00

   0  STO 01  STO 03  1  STO 02

-With h = 0.1  and  N = 10

   0.1  STO 04  10  STO 05

  XEQ "RKN6"   >>>>     x   = 1                                                                       ---Execution time = 72s---
                           RDN   y(1) =  0.536630617
                           RDN   y'(1) = -0.860171927

-The exact results, rounded to 10D are

    y(1) =  0.5366306164
    y'(1) = -0.8601719268

Note:

-Press R/S  to continue the calculations ( after changing h & N in registers R04 & R05 if need be )
 

System of 2 Differential Equations with Runge-Kutta-Nystrom of order 6
 

  "2RKN6" solves the system:

     y" = f(x,y,z)         y(x0) = y0     y'(x0) = y'0
     z" = g(x,y,z)         z(x0) = z0     z'(x0) = z'0
 

Data Registers:       R00:  subroutine name                                  ( Registers R00 thru R07  are to be initialized before executing "2RKN6" )

                                  R01 = x0         R04 = h =  stepsize                           R06 = y'0          R08 to R16: temp
                                  R02 = y0         R05 = N = number of steps               R07 = z'0
                                  R03 = z0
Flags: /
Subroutine:  a program which calculates f(x;y;z) in Y-register and g(x;y;z) in X-register
                       assuming x is in X-register , y is in Y-register and z is in Z-register upon entry.
 
 
      STACK        INPUTS      OUTPUTS
           Z             /     z(x0+N.h)
           Y             /     y(x0+N.h)
           X             /       x0+N.h

 
Example:

     d2y/dx2 =  -y.z                  y(0) = 2    y'(0) = 1         Evaluate  y(1) , z(1) , y'(1) , z'(1)    with  h = 0.1  &  N = 10
     d2z/dx2 =  x.( y + z )         z(0) = 1    z'(0) = 1
 
 
 01  LBL "T"
 02  RDN
 03  STO Z
 04  X<>Y
 05  CHS
 06  ST* Z
 07  -
 08  R^
 09  *
 10  RTN

 
  "T"  ASTO 00

    0     STO 01              1  STO 06
    2     STO 02              1  STO 07
    1     STO 03
  0.1    STO 04
   10    STO 05

 XEQ "2RKN6"  >>>>     x   =  1                                                                                            ---Execution time = 109s---
                            RDN   y(1) = 1.531356647     and   R06 = y'(1) = -2.312840139
                            RDN   z(1) = 2.620254282              R07 = z'(1) =  2.941748401

-The precision is almost perfect, since the exact results - rounded to 9D - are

      y(1) =  1.531356646      y'(1) = -2.312840137
      z(1) =  2.620254281      z'(1) =  2.941748399

-To continue the calculations, simply press R/S  ( after changing h & N in registers R04 & R05 if need be )
 

System of N Differential Equations with Gill-Runge-Kutta ( order 4 )
 

-The Gill's method is a 4-stage 4th order Runge-Kutta method to solve differential equations,
  but the coefficients are choosen so that the formulae may be written
  to use only 3 registers per variable instead of 4 registers in the "classical" RK4.
-So, we can solve larger systems of ODE.

-We have to solve the system:

     dy1/dx = f1(x,y1, ... ,yn)
     dy2/dx = f2(x,y1, ... ,yn)
             ..................                         with the initial values:     yi(x0)         i = 1 , ... , n

     dyn/dx = fn(x,y1, ... ,yn)
 

Data Registers:             R00 = subroutine name         ( Registers R00 thru R03 and R20 thru R20+n are to be initialized before executing "GRK" )

                                        R01 = n  = number of equations = number of functions              R04 thru R17: temp        R18 & R19 are unused
                                        R02 = h  = stepsize
                                        R03 = N = number of steps

                                        R20 = x0
                                        R21 = y1(x0)                                                                           R21+n thru R20+3n: temp
                                        R22 =  y2(x0)
                                          ...............         ( initial values )

                                        R20+n = yn(x0)
Flags: /
Subroutine:   A program which calculates and stores  f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in registers  R21+n, ... , R20+2n  respectively
                        with  x, y1, ... , yn in regiters R20, R21, ... , R20+n
 
 
      STACK        INPUTS      OUTPUTS
           T             /      y3(x0+N.h)
           Z             /      y2(x0+N.h)
           Y             /      y1(x0+N.h)
           X             /         x0+N.h

 
Example:      Solve the system:

      dy1/dx = - y1y2y3
      dy2/dx =  x ( y1+y2-y3 )
      dy3/dx =  xy1- y2y3

   with the initial values:   y1(0) = 1 ; y2(0) = 1 ; y3(0) = 2
 
 
 01  LBL "T"
 02  RCL 21
 03  RCL 22
 04  RCL 23
 05  *
 06  *
 07  CHS
 08  STO 24
 09  RCL 21
 10  RCL 22
 11  +
 12  RCL 23
 13  -
 14  RCL 20
 15  *
 16  STO 25
 17  LASTX
 18  RCL 21
 19  *
 20  RCL 22
 21  RCL 23
 22  *
 23  -
 24  STO 26
 25  RTN

 
 "T"  ASTO 00

 Initial values:

           0   STO 20
           1   STO 21    STO 22
           2   STO 23

 n  =  3     STO 01    ( 3 functions )
 h  = 0.1   STO 02
 N = 10   STO 03    XEQ "GRK"  >>>>      x    = 1                     = R20               ---Execution time = 145s---
                                                     RDN    y1(1) = 0.258210908  = R21                                                                         y1(1) = 0.258207906
                                                     RDN    y2(1) = 1.157620520  = R22         the exact results, rounded to 9D are          y2(1) = 1.157623981
                                                     RDN    y3(1) = 0.842179307  = R23                                                                         y3(1) = 0.842178312

-To continue the calculations, simply press R/S  ( after changing h & N in registers R02 & R03 if need be )

-There is no built-in error estimate, so we must do the calculations again with a smaller stepsize h.
-With h = 0.05 & N = 20, it yields:

    y1(1) = 0.258208088
    y2(1) = 1.157623789
    y3(1) = 0.842178380
 

Notes:

-Since the SIZE = 3n+21 instead of 4n+11 with RK4, we save registers only for n > 10
-Theoretically, "GRK" could solve a system of 99 differential equations.
-The method seems slightly more accurate than RK4 to solve the gravitational n-body problem...

Gauss-Legendre Integration ( 16-point Formula )
 

Data Registers:         R00 = Function name                       ( Registers R00 thru R03 are to be initialized before executing "GL16" )

                                    R01 = a
                                    R02 = b
                                    R03 = n = number of subintervals over which the formula is applied.  R04 = Integral         R05 thru R26: temp

Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.
 
 
 
      STACK       INPUTS     OUTPUTS
           X            /     §ab f(x).dx

 
Example:     Evaluate  §03  e-x^4 dx
 
 
 01  LBL "T"
 02  X^2
 03  X^2
 04  CHS
 05  E^X
 06  RTN

 
 "T"  ASTO 00

  0    STO 01
  3    STO 02

      n = 1    1  STO 03  XEQ "GL16"  >>>>  0.906402825                                    ---Execution time = 19s---
      n = 2    2  STO 03        R/S           >>>>  0.906402476                                    ---Execution time = 27s---
 

Double & Tripple Integrals with 3-point Gauss-Legendre Formula
 

 "DTI"  evaluates:

      §ab §u(x)v(x)  f(x,y) dx dy                                 if flag F03 is clear  and

      §ab §u(x)v(x)  §w(x,y)t(x,y)   f (x,y,z) dx dy dz    if flagF03 is set

-Gauss-Legendre 3-point formula is used
-You can divide each interval of integration by n  to get a better precision

-To save data registers and execution time,

    one subroutine must calculate u(x) in X-register & v(x) in Y-register
    and - for tripple integrals- another one must return w(x,y) in X-register and t(x,y) in Y-register
 
 

Data Registers:             R00 = f name                                   ( Registers R00 thru R04 or R05 are to be initialized before executing "DTI" )

                                        R01 = a           R04 = uv-name           R06 = Integral          R07 to R20:temp
                                        R02 = b           R05 = wt-name
                                        R03 = n

Flag:  CF 03 for double integrals
           SF 03 for tripple integrals

Subroutines:  2 or 3 subroutines:

      -The first one must take x, y and z from the X-register,the Y-register and the Z-register respectively and calculate f(x;y;z).
      -The second one takes x from the X-register and calculates u(x) and v(x) in registers X & Y respectively
      -The third one takes x and y from X and Y registers respectively and returns w(x,y) & t(x,y) in registers X & Y respectively
 
 
      STACK        INPUTS      OUTPUTS
           X             /             I

 
Example1:   Evaluate  I = §12 §xx^2 sqrt(1+x4 y4) dx dy.
 
 
 
 01  LBL "T"
 02  *
 03  X^2
 04  X^2
 05  SIGN
 06  LAST X
 07  +
 08  SQRT
 09  RTN
 10  LBL "U"
 11  X^2
 12  LASTX
 13  RTN

 
   "T" ASTO 00

    1     STO 01
    2     STO 02

  "U"  ASTO 04

         CF 03

    n = 1     1  STO 03  XEQ "DTI"  >>>   15.45937082 = R06                            ---Execution time = 10s---
    n = 2     2  STO 03      R/S          >>>   15.46673275
    n = 4     4  STO 03      R/S          >>>   15.46686031
    n = 8     8  STO 03      R/S          >>>   15.46686245  ... all the digits are correct!

Note:

-When n is multiplied by 2, execution time is roughly multiplied by 4)
 

Example2:      Evaluate   I = §12 §xx^2 §x+yxy  xyz / sqrt ( x2 + y2 + z2 ) dx dy dz
 
 
 01  LBL "T"
 02  ENTER^
 03  X^2
 04  R^
 05  ST* Z
 06  X^2
 07  +
 08  R^
 09  ST* Z
 10  X^2
 11  +
 12  SQRT
 13  /
 14  RTN
 15  LBL "U"
 16  X^2
 17  LASTX
 18  RTN
 19  LBL "W"
 21  STO Z
 22  X<>Y
 23  ST* Z
 24  +
 25  RTN

 
    " T"   ASTO 00

      1      STO 01
      2      STO 02

    "U"    ASTO 04
    "W"   ASTO 05

        SF 03

         n = 1      1  STO 03    XEQ "DTI"  >>>>    I1 =  0.765014888                                    ---Execution time = 38s---
         n = 2      2  STO 03         R/S         >>>>    I2 =  0.770640690                                   ---Execution time = 4m06s---
         n = 4      4  STO 03         R/S         >>>>    I4 =  0.770731245
         n = 8      8  STO 03         R/S         >>>>    I8 =  0.770732669

Notes:

-The exact value is 0.7707326901 rounded to 10D
-Unfortunately, with triple integrals, when n is multiplied by 2, execution time is roughly multiplied by 8 !
 

Multiple Integrals with 3-point Gauss-Legendre Formula
 

 "MIT" also uses Gauss-Legendre 3-point formula to estimate:

    §ab §u1(x1)v1(x1)  §u2(x1,x2)v2(x1,x2) .........  §u(n-1)(x1,...,x(n-1))v(n-1)(x1,...,x(n-1))   f(x1,x2,....;xn)  dx1dx2....dxn   ( n < 7 )

    i-e integrals , double integrals , ...... , up to sextuple integrals, provided the subroutines do not contain any XEQ.

-This limitation is only due to the fact that the return stack can accomodate up to 6 pending return addresses.
 

>>> The same trick as above is used to save execution time: only one subroutine for each pair of endpoints
 
 

Data Registers:       R00 = m =  number of  §

                                 R01 = x1 ;  R02 = x2 ; .......... ;  R06 = x6

                                 R07 = 7 ; R08 = 4 ; R09 = 2 ; R10 = 1.6 ; R11 = 3.6 ; R12 = 0.151/2 ; R13 = 0.61/2   ;    R14 thru R17 = control numbers

                     R18 = f name

                     R19 = n = number of sub-intervals       R25 = u1v1 name         R31 = u2v2 name         .....................     R49 = u5v5 name
                     R20 = a                                               R26 = u1(x1)               R32 = u2(x1;x2)           .....................     R50 = u5(x1;x2;...;x5)
                     R21 = b                                               R27 = v1(x1)               R33 = u2(x1;x2)           .....................     R51 = v5(x1;x2;...;x5)
                     R22 = (b-a)/n                                      R28 = (v1-u1)/n           R34 = (v2-u2)/n           .....................     R52 = (v5-u5)/n
                     R23 = index                                         R29 = index                R35 = index                .....................     R53 =  index
                     R24 = partial sum,                                R30 = partial sum        R36 = partial sum       .....................      R54 = partial sum
               and, finally, the integral
 

Flags: /
Subroutines:  The functions  f ;  uv1 ; uv2 ; ....... are to be computed assuming  x1 is in R01 ; x2 is in R02 ; ........

Example:    Evaluate   I  =  §13  §xx^2  §x+yx.y  §zx+z   ln(x2+y/z+t) dx.dy.dz.dt
 
 
 01 LBL "T"
 02 RCL 01
 03 X^2
 04 RCL 02
 05 RCL 03
 06 /
 07 +
 08 RCL 04
 09 +
 10 LN
 11 RTN
 12 LBL "U"
 13 RCL 01
 14 X^2
 15 LASTX
 16 RTN
 17 LBL "V"
 18 RCL 01
 19 RCL 02
 20 ST* Y
 21 RCL 01
 22 +
 23 RTN
 24 LBL "W"
 25 RCL 01
 26 RCL 03
 27 ST+ Y
 28 RTN
 
 
      XEQ "MIT"            >>>>    " A^B"
  1 ENTER 3   R/S       >>>>     " FNAME?"

         T   R/S      >>>>      " UV1 NAME?"
         U   R/S     >>>>      " UV2 NAME?"
         V   R/S     >>>>      " UV3 NAME?"
         W   R/S     >>>>     " UV4NAME?"      press  R/S   without any entry when all function names have been keyed in

                R/S     >>>>     " N=?"                    ( number of sub-intervals )  let's try     n = 1

           1   R/S     >>>>       I = 160.4523151 = R24                                        ---Execution time = 2m59s---

-To recalculate I with another n-value, key in this value and press  R/S or XEQ 10

  for example,

    with  n = 2     2   R/S   >>>>   160.6314953
    and   n = 4     4   R/S   >>>>   160.6342731

Notes:

-If n is multiplied by 2, execution time is approximately multiplied by 16 for quadruple integrals , by 64 for sextuple integrals.
-We can use Lagrange interpolation formula to improve our results by extrapolation to the limit.
-In this example, we find  I = 160.6343171
 

Multiple Integrals with all borns of integration -1+1
 

-Again Gauss-Legendre 3-point formula to estimate

    §-1+1 §-1+1 .......................... §-1+1   f(x1,x2,....;xn)  dx1dx2....dxn
 

Data Registers:         R00 = n =  number of  §                                   ( Register R00 is to be initialized before executing "MIT1" )

                                     R01 = x1 ,  R02 = x2 , .......... ,  Rnn = xn      Rn+1 ..... R2n:  counters    R2n+1 to R2n+6: temp

Flags: /
Subroutine:  A program  LBL "T"  to compute f(x1,x2,....;xn) assuming  x1 is in R01 , x2 is in R02 , ........
 
 
 
      STACK        INPUT      OUTPUT
           X             /             I

 
Example:     Evaluate     I  =  §-11  §-11  §-11  §-11   Ln ( PI^2 + x + y + z + t )  dx dy dz dt
 
 
 01 LBL "T"
 02 PI
 03 X^2
 04 RCL 01
 05 RCL 02
 06 RCL 03
 07 +
 08 +
 09 +
 10 RCL 04
 11 +
 12 LN
 13 RTN

 
   4  STO 00   XEQ "MIT1"  >>>>  I = 36.51975046                             ---Execution time = 4m26s---
 

Notes:

-If the endpoints are different from -1 +1, make a linear change of argument.

-The intervals of integration cannot be divided, so the precision will remain low in most cases.

-The number of § is only limited by the memory
-Practically, due to execution time, it would be difficult to evaluate centuple integrals...
 

Integration on the Surface of a Sphere
 

 "ISSPH" employs one of the formulae given in reference [7] to estimate

      I = §§x^2+y^2+z^2=R^2    f(x,y,z)  ds    where ds is the element of surface on a sphere
 

Data Registers:             R00 = function name                                      ( Register R00 is to be initialized before executing "ISSPH" )

                                         R01 = Integral   R02 = R   R03 to R08: temp
Flags: /
Subroutine:  A program that takes x , y , z in registers  X , Y , Z  respectively and returns  f(x,y,z) in X-register.
 
 
 
      STACK        INPUT      OUTPUT
           X             R       Integral

    where R = radius of the sphere

Example:    I = §§x^2+y^2+z^2=R^2    Ln ( 16 + x3 + y2 + z )  dx dy dz     with  R = 2
 
 
 01 LBL "T"
 02 X^3
 03 X<>Y
 04 X^2
 05 +
 06 +
 07 16
 08 +
 09 LN
 10 RTN

 
  "T"  ASTO 00

   2   XEQ "ISSPH"  >>>>  I = 142.2311076                                         ---Execution time = 79s---

Notes:

-Exact result = 142.231086....
-The formula is exact if f is a polynomial and deg f < 14
 

Norm of a 3D-Vector
 

-Just a short M-Code routine to calculate the norm of a 3D-vector in the stack
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             z             z
           Y             y             y
           X             x   sqrt(x2+y2+z2)
           L             /            x

 
Example:     V(3,4,5)

  5  ENTER^
  4  ENTER^
  3  XEQ "NORM"  >>>  7.071067812

Note:

-There is no check for alpha data.
 

2 M-code Routines
 

 X^3  takes x and returns x^3 in X-register. x is saved in L-register
 X/3   replaces x by X/3 in X-register. L-register is unchanged.

 No check for alpha data.
 
 
 

References:

[1]  Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  https://en.wikipedia.org/w/index.php?title=Differential_geometry_of_curves&oldid=740874421
[3]  https://en.wikipedia.org/w/index.php?title=Frenet%E2%80%93Serret_formulas&oldid=731811089
[4]  Ron Goldman - "Curvature formulas for implicit curves and surfaces"
[5]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5
[6]  F. Scheid -  "Numerical Analysis" - Mc Graw Hill   ISBN 07-055197-9
[7]  Z. P. Bazant - "Efficient Numerical Integration on the Surface of a Sphere"