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+n² are to be initialized before executing "DET" )

                                •  Rbb   = a₁₁  ....................   •  Ree   = a_nn

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

    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³ + O(h¹⁰)

       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 = d³f/dx³    R06 to R08: temp
                                         R02 = h     R04 = d²f/dx²

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:

   h⁴ f⁽⁴⁾(x) ~  ( 192654 F - 140196 G + 52428 H - 9738 I +1261 J - 82 K ) / 15120
   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
 

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

                                         R01 = x     R03 = f(x)         R05 = d⁵f/dx⁵    R07 to R09: temp
                                         R02 = h     R04 = d⁴f/dx⁴   R06 = d⁶f/dx⁶

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 = ¶²f / ¶x² ; f "_yy = ¶²f / ¶y²  and   f "_xy  = ¶²f / ¶x¶y
 
 

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

                    R01 = x                     R03 = f '_x    = ¶f / ¶x         R06 =  f '_y   = ¶f / ¶y           R09 =  f "_xy  = ¶²f / ¶x¶y
                    R02 = y                     R04 = f "_xx   = ¶²f / ¶x²      R07 = f "_yy  = ¶²f / ¶y²
                                                     R05 = f "'_xxx = ¶³f / ¶x³      R08 = f "'_yyy = ¶³f / ¶y³

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(-x²).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  = ¶²f / ¶x² =  0.509989188      the exact result is     0.509989195
                               RDN  f "_yy = ¶²f / ¶y² = -0.091969868      the exact result is    -0.091969860
                         LASTX  f "_xy  = ¶²f / ¶x¶y = -0.367879439    the exact result is    -0.367879441

-You have also:

       R05 = f "'_xxx = ¶³f / ¶x³   = 1.019980483  exact result = 1.019978390
       R08 = f "'_yyy = ¶³f / ¶y³ =  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 = ¶²f / ¶x² ;  f "_yy = ¶²f / ¶y² ;  f "_zz = ¶²f / ¶z²

   and  f "'_xxx = ¶³f / ¶x³      f "'_yyy = ¶³f / ¶y³      f "'_zzz = ¶³f / ¶z³
 

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 = ¶²f / ¶x²       R07 = f "_yy = ¶²f / ¶y²       R10 = f "_zz = ¶²f / ¶z²      R13 = x  R14 = y  R15 = z
                                         R05 = f "'_xxx = ¶³f / ¶x³     R08 = f "'_yyy = ¶³f / ¶y³     R11 = f "'_zzz = ¶³f / ¶z³

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 ) = ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z²

Example:     f(x) = exp(-x²).Ln(y²+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 = ¶²f / ¶x² =  1.413720682             ---Execution time = 52s---
                                RDN   f '_y  = ¶f / ¶y  =    0.210216822              R07 = f "_yy = ¶²f / ¶y² = -0.015015424
                                RDN   f '_z  = ¶f / ¶z   =    0.052554204              R10 = f "_zz =  ¶²f / ¶z² = -0.007507569
    RDN    ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z² =     1.409197689      and

     R05 = f "'_xxx  = ¶³f / ¶x³  =  2.863447421
     R08 = f "'_yyy = ¶³f / ¶y³  = -0.042900947
     R11 = f "'_zzz  = ¶³f / ¶z³ =   0.002145569
 

Mixed Derivative of a Function of 2 variables
 

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

     f "_xy = ¶²f / ¶x¶y = (1/h²) SUM_j=1,2,3,4,5  A_j ( f_jj + f_-j-j - f_j-j - f_-jj )

   where  A₁ = 5/12     A₂ = -5/84     A₃ = 5/504     A₄ = -5/4032     A₅ = 1/12600
 

-This formula is exact for any polynomial of degree < 11               (   f_km =  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 = ¶²f / ¶x¶y     R06: temp
                                         R03 = h
Flags: /
Subroutine:         A program which computes f(x,y) assuming x is in X-register and y is in Y-register upon entry.
 
 

STACK	INPUTS	OUTPUTS
Z	h	/
Y	y	/
X	x	f "xy = ¶2f / ¶x¶y

 
Example:     f(x) = exp(-x²).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  = ¶²f / ¶x¶y = -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 d³F/dx²dy & d³F/dxdy²
 

 "MDV21" & "MDV12" 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
 

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 = ¶³f / ¶x²¶y  or  f'''_xyy = ¶³f / ¶y²¶x     R05-R06  temp

Flag:  F02            CF 02 ->  f'''_xxy = ¶³f / ¶x²¶y
                             SF 02  ->  f'''_xyy = ¶³f / ¶y²¶x

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 + x² 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 = ¶³f / ¶x²¶y  =  - 0.370369497                ---Execution time = 25s---

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

-With the same inputs,  XEQ "MDV12"  returns   f'''_xyy = ¶³f / ¶y²¶x  =  - 0.148148268

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

Mixed Derivative d³F/dxdydz
 

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

    h³ f'''_xyz ~  SUM_k=1,2,...,5  A_k [  f_kkk - f_-k-k-k - f_kk-k + f_-k-kk - f_k-kk + f_-kk-k - f_-kkk + f_k-k-k ) ]

 With    A₁ = 5/24 ,  A₂ = -5/336  ,  A₃ = +5/3024  ,  A₄ = -5/32256  ,  A₅ = +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 = ¶³f / ¶x¶y¶z
                                         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 + x² + 2.y + z³ )    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 = ¶³f / ¶x¶y¶z  =  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 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
 

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

                                         R01 = x     R03 = h
                                         R02 = y     R04 = f""_xxyy = ¶⁴f / ¶x²¶y²     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 ( x² + y³ )     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⁽⁴⁾_xxyy =  ¶⁴f / ¶x²¶y² =  - 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 = x₁ ,  •  R02 = x₂ , .......... ,  •  Rnn = x_n      n < 11

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

Flags:  /
Subroutine:   A program which computes  f(x₁, ..... ,x_n)  assuming  x₁ is in R01 , ............ , x_n is in Rnn
 
 

STACK	INPUTS	OUTPUTS
Y	h	/
X	i	f 'xi = ¶f / ¶xi

 
Example:     f(x) = exp(-x²).Ln(y²+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 / ¶x₁ = -0.509989198   ( The last decimal should be a 5 )                    ---Execution time = 10s---

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

  0.1   ENTER^
    3       R/S      >>>>   f '_x3 = ¶f / ¶x₃ = 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 = x₁ ,  •  R02 = x₂ , .......... ,  •  Rnn = x_n

                                           R10 = h      R11 temp             R13 = x_i           R15 = 0 , h , 2h , ....
                                           R12 =  ¶²f / ¶x_i¶x_j        R14 = x_j            R16 =  i    R17 = j   R18: temp
Flag:  F10
Subroutine:   A program which computes  f(x₁, ..... ,x_n)  assuming  x₁ is in R01 , ............ , x_n is in Rnn
 
 
 

STACK	INPUTS	OUTPUTS
Z	h	/
Y	j	/
X	i	f "ij = ¶2f / ¶xi¶xj

 
Example:       f(x) = exp(-x²).Ln(y²+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 = ¶²f / ¶x² = 0.509989206   ( exact derivative = 0.509989195 )                   ---Execution time = 12s---

 0.1   ENTER^
   1    ENTER^
   2       R/S     >>>>   f "_xy = ¶²f / ¶x¶y = -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 = x₁ ,  •  R02 = x₂ , .......... ,  •  Rnn = x_n

                                           R10 = h              R12-R18:  temp
                                           R11 =  ¶³f / ¶x_i¶x_j¶x_k
Flags:  F09-F10
Subroutine:   A program which computes  f(x₁, ..... ,x_n)  assuming  x₁ is in R01 , ............ , x_n is in Rnn
 
 
 

STACK	INPUTS	OUTPUTS
Z	h	/
Y	j	/
X	i	f "ij = ¶2f / ¶xi¶xj

 
Example:             f(x,y,z,t) = exp(-x²) Ln( y² + z + t³ )     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 = ¶³f / ¶x¶y¶z  = ?

-With   h = 0.1

     0.1  ENTER^
      1    ENTER^
      2    ENTER^
      3    XEQ "TD"   >>>>    f'''_xyz = ¶³f / ¶x¶y¶z  = 0.045984935                                          ---Execution time = 57s intead of 88s---

-Exact result = 0.045984930....

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

-With   h = 0.1

     0.1  ENTER^
      1    ENTER^
      4    ENTER^    XEQ "TD"   >>>>    f'''_xtt = ¶³f / ¶x¶t²  = - 0.448353739                                    ---Execution time = 42s---

-Exact result = - 0.448353069....

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

-With   h = 0.1

     0.1  ENTER^
      2    ENTER^   ENTER^    XEQ "TD"   >>>>    f'''_yyy = ¶³f / ¶y³  = - 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  G^t      where  G^t  = transposed G              So, we get a square matrix W of order (n+1)
                 G  0

  with        H = [ h_i,j ]   with  h_i,j = ¶²f / ¶x_i¶x_j
                G = [ g_k ]    with  g_k = ¶f / ¶x_i

-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 = x₁ ,  •  R02 = x₂ , .......... ,  •  Rnn = x_n

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

Flag:  F10
Subroutines:   "FD"  "SD"  and a program which computes  f(x₁, ..... ,x_n)  assuming  x₁ is in R01 , ............ , x_n is in Rnn
 
 

STACK	INPUTS	OUTPUTS
Y	h	/
X	n < 11	26.eeemm

   Where  m = n+1

Example:    f(x,y,z) = x⁴ y³ z² - 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(-x².t).Ln(y²+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) =  ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z² + ¶²f / ¶t² = 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/c²) ¶²f / ¶t² - ( ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z² )       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) ( x² + y² + z² )       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  ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z² - (1/c²) ¶²f / ¶t²
-Simply change the sign of the result.
 

Heat Operator
 

-It is defined by   ¶f / ¶t - D ( ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z² )      where D is a positive constant
 
 

STACK	INPUTS	OUTPUTS
Y	h	/
X	D	Heat(f)

 
Example:      f(x,y,z,t) = exp(-t) ( x² + y² + z² )       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

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

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

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

                                         R01 = x        R04 = h
                                         R02 = y        R05 = D² 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(-x²).Ln(y²+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"  >>>> D² 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
 

   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
 
 

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

                                         R01 = x        R04 = h
                                         R02 = y        R05 = D³ 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(-x²).Ln(y²+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"  >>>> D³ 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  §_a^b  [ ( dX/dt )² + ( dY/dt )² + ( dZ/dt )² ] ^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) = t⁴ ,  Y(t) = t²  ,  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 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²
 
 

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 + x² y³ )    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"  >>>>   K_G = -0.124570010                             ---Execution time = 49s instead of 58s---
                                                     X<>Y  K_m = -0.171204283

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.
 

Gaussian Curvature & Mean Curvature of a HyperSurface
 

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

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
 
 

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

                              •  R01 = x₁ ,  •  R02 = x₂ , .......... ,  •  Rnn = x_n

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

Flag:  F10
Subroutines:  "HESSGR" and a program which computes  f(x₁, ..... ,x_n)  assuming  x₁ is in R01 , ............ , x_n is in Rnn
 
 

STACK	INPUTS	OUTPUTS
Z	/	det W
Y	h	Km
X	n < 11	KG

 
Example:   The surface ( n = 3 )   x⁴ y³ z² - 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"  >>>>      K_G  =  0.256837099                             ---Execution time = 137s---
                                     RDN      K_m  =  0.518666290
                                    X<>Y  det W = - 216

Warning:

-This program does not check that  f(x₁, ..... ,x_n) = 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' |³                                      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
 
 

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

-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 x₁(t) , x₂(t) , ........ , x_n(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

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

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  = ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z²      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  = ¶²f / ¶r² + (1/r) ¶f/¶r + (1/r²) ¶²f / ¶f² + ¶²f / ¶z²      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  = ¶²f / ¶r² + (2/r) ¶f/¶r + (1/r²) ¶²f / ¶q² + (1/(r²tan q))  ¶f / ¶q + (1/(r²sin²q)) ¶²f / ¶f²      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(-x²) Ln(y²+z) , x²y²z² , exp(x) y²z ]     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 z² sin²f , r² z , r³ 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 sin² q cos² f , r²sin f  , r³ cos q cos² 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(x₀) = y₀
 

-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 = x₀      •   R03 =  h = the stepsize                      R05  to R12: temp
                                  •  R02 = y₀      •   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 2⁷ = 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(x₀) = y₀
       z' = g(x,y,z)                                   z(x₀) = z₀
 

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

                                 •  R01 = x₀       •  R04 = h = stepsize                       R06 to R19: temp    R55: counter
                                 •  R02 = y₀       •  R05 = N = number of steps          R20 to R54: Coefficients
                                 •  R03 = z₀
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 2⁷ = 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(x₀) = y₀     y'(x₀) = y'₀

-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 = x₀      •   R04 = h  = stepsize                   R06 to R10: temp
                               •   R02 = y₀      •   R05 = N = number of steps
                               •   R03 = y'₀
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 ( x² + y² ) ^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(x₀) = y₀     y'(x₀) = y'₀
     z" = g(x,y,z)         z(x₀) = z₀     z'(x₀) = z'₀
 

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

                               •   R01 = x₀      •   R04 = h =  stepsize                        •   R06 = y'₀          R08 to R16: temp
                               •   R02 = y₀      •   R05 = N = number of steps            •   R07 = z'₀
                               •   R03 = z₀
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:

     d²y/dx² =  -y.z                  y(0) = 2    y'(0) = 1         Evaluate  y(1) , z(1) , y'(1) , z'(1)    with  h = 0.1  &  N = 10
     d²z/dx² =  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:

     dy₁/dx = f₁(x,y₁, ... ,y_n)
     dy₂/dx = f₂(x,y₁, ... ,y_n)
             ..................                         with the initial values:     y_i(x₀)         i = 1 , ... , n

     dy_n/dx = f_n(x,y₁, ... ,y_n)
 

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 = x₀
                                      •  R21 = y₁(x₀)                                                                           R21+n thru R20+3n: temp
                                      •  R22 =  y₂(x₀)
                                          ...............         ( initial values )

                                      •  R20+n = y_n(x₀)
Flags: /
Subroutine:   A program which calculates and stores  f₁ ( x, y₁ , ... , y_n) , ... , f_n ( x, y₁ , ... , y_n)  in registers  R21+n, ... , R20+2n  respectively
                        with  x, y₁, ... , y_n 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:

      dy₁/dx = - y₁y₂y₃
      dy₂/dx =  x ( y₁+y₂-y₃ )
      dy₃/dx =  xy₁- y₂y₃

   with the initial values:   y₁(0) = 1 ; y₂(0) = 1 ; y₃(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    y₁(1) = 0.258210908  = R21                                                                         y₁(1) = 0.258207906
                                                     RDN    y₂(1) = 1.157620520  = R22         the exact results, rounded to 9D are          y₂(1) = 1.157623981
                                                     RDN    y₃(1) = 0.842179307  = R23                                                                         y₃(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:

    y₁(1) = 0.258208088
    y₂(1) = 1.157623789
    y₃(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  §₀³  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:

      §_a^b §_u(x)^v(x)  f(x,y) dx dy                                 if flag F03 is clear  and

      §_a^b §_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 = §₁² §_x^{x^2} sqrt(1+x⁴ y⁴) 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 = §₁² §_x^{x^2} §_x+y^xy  xyz / sqrt ( x² + y² + z² ) 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"  >>>>    I₁ =  0.765014888                                    ---Execution time = 38s---
         n = 2      2  STO 03         R/S         >>>>    I₂ =  0.770640690                                   ---Execution time = 4m06s---
         n = 4      4  STO 03         R/S         >>>>    I₄ =  0.770731245
         n = 8      8  STO 03         R/S         >>>>    I₈ =  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:

    §_a^b §_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(x₁,x₂,....;x_n)  dx₁dx₂....dx_n   ( 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 = x₁ ;  R02 = x₂ ; .......... ;  R06 = x₆

                                 R07 = 7 ; R08 = 4 ; R09 = 2 ; R10 = 1.6 ; R11 = 3.6 ; R12 = 0.15^1/2 ; R13 = 0.6^1/2   ;    R14 thru R17 = control numbers

                     R18 = f name

                     R19 = n = number of sub-intervals       R25 = u₁v₁ name         R31 = u₂v₂ name         .....................     R49 = u₅v₅ name
                     R20 = a                                               R26 = u₁(x₁)               R32 = u₂(x₁;x₂)           .....................     R50 = u₅(x₁;x₂;...;x₅)
                     R21 = b                                               R27 = v₁(x₁)               R33 = u₂(x₁;x₂)           .....................     R51 = v₅(x₁;x₂;...;x₅)
                     R22 = (b-a)/n                                      R28 = (v₁-u₁)/n           R34 = (v₂-u₂)/n           .....................     R52 = (v₅-u₅)/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 ;  uv₁ ; uv₂ ; ....... are to be computed assuming  x₁ is in R01 ; x₂ is in R02 ; ........

Example:    Evaluate   I  =  §₁³  §_x^{x^2}  §_x+y^x.y  §_z^x+z   ln(x²+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⁺¹   f(x₁,x₂,....;x_n)  dx₁dx₂....dx_n
 

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

                                     R01 = x₁ ,  R02 = x₂ , .......... ,  Rnn = x_n      Rn+1 ..... R2n:  counters    R2n+1 to R2n+6: temp

Flags: /
Subroutine:  A program  LBL "T"  to compute f(x₁,x₂,....;x_n) assuming  x₁ is in R01 , x₂ is in R02 , ........
 
 
 

STACK	INPUT	OUTPUT
X	/	I

 
Example:     Evaluate     I  =  §_-1¹  §_-1¹  §_-1¹  §_-1¹   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 + x³ + y² + 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"