hp41programs

Derive

Numerical Differentiation for the HP-41


Overview
 

 1°) Functions of 1 variable

  a)  First & Second Derivatives
  b)  First & Second Derivatives - More accurate formulae

 2°) Functions of 2 variables

  a)  First & Second Derivatives
  b)  Mixed Derivative

    b-1) Program#1
    b-2) Program#2 - more accurate formula

 3°) Functions of 3 variables

 4°) An Iterative Method

 5°) Extrapolation to the limit

 6°) Functions of n variables ( n < 11 )

  a) First Derivatives
  b) Second Derivatives
  c) Iterative Method
  d) Extrapolation to the limit
 
 

1°) Functions of 1 variable
 

     a) First & Second Derivatives
 

Formulae:        df/dx   = (1/12h).[ f(x-2h) - 8.f(x-h) + 8.f(x+h) - f(x+2h) ] + O(h4)
                        d2f/dx2 = (1/12h2).[ - f(x-2h) + 16.f(x-h) - 30.f(x) +16.f(x+h) - f(x+2h) ] + O(h4)

-These formulas ( without + O(h4) ) are exact for any polynomial of degree < 5
 

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

                                         R01 = x        R03 = f(x)
                                         R02 = h        R04 = f '(x)      R05 = f "(x)
Flags: /
Subroutine:      A program which computes f(x) assuming x is in X-register upon entry
 
 

01  LBL "dF"
02  STO 01
03  X<>Y
04  STO 02
05  +
06  XEQ IND 00
07  STO 04
08  STO 05
09  RCL 01
10  RCL 02
11  -
12  XEQ IND 00
13  ST- 04
14  ST+ 05
15  8
16  ST* 04
17  ST+ X
18  ST* 05
19  RCL 01
20  RCL 02
21  ST+ X
22  STO 03
23  +
24  XEQ IND 00
25  ST- 04
26  ST- 05
27  RCL 01
28  RCL 03
29  -
30  XEQ IND 00
31  ST+ 04
32  ST- 05
33  RCL 01
34  XEQ IND 00
35  STO 03
36  30
37  *
38  ST- 05
39  RCL 02
40  ST/ 05
41  12
42  *
43  ST/ 04
44  ST/ 05
45  RCL 05        
46  RCL 04
47  END

 
   ( 75 bytes / SIZE 007 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h          f "(x)
           X             x          f '(x)

 
Example:     f(x) = exp(-x2)   Calculate  f '(1) & f "(1)

  LBL "T"          any global LBL , 6 characters or less
  X^2
  CHS
  E^X
  RTN               T  ASTO 00

-If we choose h = 0.03

  0.03  ENTER^
    1     XEQ "dF"  >>>>   f '(1) = -0.735758961     the exact values are  f '(1) = -0.735758882
                            X<>Y   f "(1) = 0.735757408                              and  f "(1) = 0.735758882

-Choosing the best h-value is not easy but h ~ 0.03 "often" produces good results.
 

     b) First & Second Derivatives - More Accurate Formulae
 

-The "Handbook of Mathematical Functions" gives many formulas ( table 25.2 ) ,  up to order 6
-But we can also build formulas of higher order and the following program uses formulas of order 10:

       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)

-These are exact for any polynomial of degree < 11     (  f(x+k.h) is denoted fk to simplify these expressions )
 

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

                                         R01 = x        R03 = f(x)
                                         R02 = h        R04 = f '(x)      R05 = f "(x)
Flags: /
Subroutine:      A program which computes f(x) assuming x is in X-register upon entry
 
 

  01  LBL "dF+"
  02  STO 01
  03  X<>Y
  04  STO 02
  05  +
  06  XEQ IND 00
  07  STO 04
  08  STO 05
  09  RCL 01
  10  RCL 02
  11  -
  12  XEQ IND 00
  13  ST- 04
  14  ST+ 05
  15  42
  16  ST* 05
  17  ST+ X
  18  ST* 04
  19  RCL 01
  20  RCL 02
  21  ST+ X
  22  +
  23  XEQ IND 00
  24  STO 03
  25  RCL 01
  26  RCL 02
  27  ST+ X
  28  -
  29  XEQ IND 00
  30  STO Y
  31  RCL 03
  32  -
  33  24
  34  *
  35  ST+ 04
  36  X<> 03
  37  +
  38  6
  39  *
  40  ST- 05
  41  RCL 01
  42  RCL 02
  43  3
  44  *
  45  +
  46  XEQ IND 00
  47  STO 03
  48  ST+ 05
  49  RCL 01
  50  RCL 02
  51  3
  52  *
  53  -
  54  XEQ IND 00
  55  ST+ 05
  56  RCL 03
  57  -
  58  6
  59  *
  60  ST- 04
  61  RCL 01
  62  RCL 02
  63  4
  64  *
  65  +
  66  XEQ IND 00
  67  STO 03
  68  ST- 04
  69  RCL 01
  70  RCL 02
  71  4
  72  *
  73  -
  74  XEQ IND 00
  75  ST+ 04
  76  RCL 03
  77  +
  78  8
  79  /
  80  ST- 05
  81  25
  82  ST* 04
  83   E3
  84  ST* 05
  85  RCL 01
  86  RCL 02
  87  5
  88  *
  89  +
  90  XEQ IND 00
  91  STO 03
  92  ST+ X
  93  ST+ 04
  94  RCL 01
  95  RCL 02
  96  5
  97  *
  98  -
  99  XEQ IND 00
100  ST- 04
101  ST- 04
102  RCL 03
103  +
104  8
105  *
106  ST+ 05
107  RCL 01
108  XEQ IND 00
109  STO 03
110  73766
111  *
112  ST- 05
113  RCL 02
114  ST/ 05
115  2520
116  *
117  ST/ 04
118  10
119  *
120  ST/ 05
121  RCL 05
122  RCL 04
123  END

 
   ( 182 bytes / SIZE 006 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h          f "(x)
           X             x          f '(x)

 
Example:     f(x) = exp(-x2)   Calculate  f '(1) & f "(1)
 

  LBL "T"
  X^2
  CHS
  E^X
  RTN               T  ASTO 00

-If we choose h = 0.1

   0.1   ENTER^
    1     XEQ "dF+"  >>>>   f '(1) = -0.735758886   (  the exact values are  f '(1) = -0.735758882
                              X<>Y   f "(1) = 0.735758849                                and  f "(1) = 0.735758882  )

    h ~ 0.1   "often" produces very good results for the first derivative ( 8 or 9 digit accuracy ), a little less good for the second derivative.
 

Note:   "dF" or "dF+" can be used to compute the radius of curvature of a curve defined by  y = f(x). For instance
 
 

 01  LBL "RHO"
 02  XEQ "dF+"
 03  X^2
 04  1
 05  +
 06  ENTER^
 07  SQRT
 08  *
 09  X<>Y
 10  /
 11  ABS
 12  END

 
-With y =  exp(-x2)   and x = 1 ,  it yields:   rho = 2.600834159   ( exact rho = 2.600834029 )
 

2°) Functions of 2 variables
 

     a) First & Second Derivatives
 

-"dF2" calls "dF+" or "dF" to compute the 4 partial derivatives  f 'x = f / x ; f 'y = f / y ;  f "xx = 2f / x2 ; f "yy = 2f / y2
 

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

                                         R01 = x        R03 = h            R05 = f 'x = f / x        R07 = f "xx = 2f / x2
                                         R02 = y        R04 = f(x,y)      R06 = f 'y = f / y        R08 = f "yy = 2f / y2          R09 & R10: temp

Flags: /
Subroutines:   "dF+" or "dF" and a program which computes f(x,y) assuming x is in X-register and y is in Y-register upon entry.
 
 

01  LBL "dF2"   
02  STO 07
03  RDN
04  STO 09
05  RCL 00
06  STO 10
07  "Y"
08  ASTO 00
09  RDN
10  XEQ "dF+"  
11  STO 06
12  X<>Y
13  STO 08
14  "X"
15  ASTO 00
16  RCL 02
17  RCL 07
18  XEQ "dF+"
19  STO 05
20  X<>Y
21  STO 07       
22  RCL 10
23  STO 00
24  RCL 09
25  X<> 02
26  X<> 03
27  STO 04
28  RCL 08
29  RCL 07
30  RCL 06       
31  RCL 05
32  CLA
33  RTN
34  LBL "X"
35  RCL 09
36  X<>Y
37  GTO IND 10
38  LBL "Y"
39  RCL 07
40  GTO IND 10
41  END

 
   ( 73 bytes / SIZE 011 )
 
 

      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

  LBL "T"
  X^2
  CHS
  E^X
  X<>Y
  LN
  *
  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
                               RDN   f 'y  = f / y  =    0.183939720      the exact result is     0.183939721
                               RDN  f "xx  = 2f / x2 =  0.509989206      the exact result is     0.509989195
                               RDN  f "yy = 2f / y2 = -0.091969921      the exact result is    -0.091969860

-Whence  Laplacian ( f ) = 2f / x2 + 2f / y2 =  0.418019286   ( exact value = 0.418019335 )
 

-In this example, execution time = 29 seconds
-Of course, the results are not so accurate if you use "dF" instead of "dF+"
 
 

     b) Mixed Derivative
 

        b-1) Program#1
 

-The following routine employs a method of order 4

      f "xy = 2f / xy = (1/24.h2).[ 30 f00 + 16.( f11 + f -1-1 - f10 - f -10 - f01 - f0-1 ) + ( - f 22 - f -2-2 + f20 + f -20 + f02 + f0-2 ) ] + O(h4)

-This formula is exact for any polynomial of degree < 5               (   fkm denotes  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.
 
 

01  LBL "MDV" 
02  STO 01
03  RDN
04  STO 02
05  X<>Y
06  STO 03
07  CLX
08  STO 06
09  XEQ 01
10  16
11  *
12  STO 05
13  XEQ 01
14  ST- 05
15  RCL 02
16  RCL 01
17  XEQ IND 00
18  STO 04
19  30
20  *
21  RCL 05
22  +
23  24
24  RCL 03
25  X^2
26  *
27  /
28  STO 05
29  RTN
30  LBL 01
31  RCL 03
32  ST+ 06
33  RCL 02
34  RCL 01
35  RCL 06
36  +
37  XEQ IND 00
38  STO 04
39  RCL 02
40  RCL 01
41  RCL 06
42  -
43  XEQ IND 00
44  ST+ 04
45  RCL 02
46  RCL 06
47  +
48  RCL 01
49  XEQ IND 00
50  ST+ 04
51  RCL 02
52  RCL 06
53  -
54  RCL 01
55  XEQ IND 00
56  ST+ 04
57  RCL 02
58  RCL 01
59  RCL 06
60  ST+ Z
61  +
62  XEQ IND 00
63  ST- 04
64  RCL 02
65  RCL 01
66  RCL 06
67  ST- Z
68  -
69  XEQ IND 00
70  RCL 04
71  -
72  END

 
   ( 102 bytes / SIZE 007 )
 
 

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

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

  LBL "T"
  X^2
  CHS
  E^X
  X<>Y
  LN
  *
  RTN               T  ASTO 00

-If we choose h = 0.03

  0.03  ENTER^
    2     ENTER^
    1     XEQ "MDV"  >>>>   f "xy = 2f / xy = -0.367879722      ( exact value = -1/e = -0.367879441 )
 
 

        b-2) Program#2 - more accurate formula
 

-Like "dF+" , "MDV+" uses a method of order 10:

     f "xy = 2f / xy = (1/50400.h2).[ 73766 f00 + 42000.( f11 + f -1-1 - f10 - f -10 - f01 - f0-1 ) + 6000.( - f 22 - f -2-2 + f20 + f -20 + f02 + f0-2 )
                                                         + 1000.( f33 + f -3-3 - f30 - f -30 - f03 - f0-3 ) + 125.( - f 44 - f -4-4 + f40 + f -40 + f04 + f0-4 )
                                                               + 8.( f55 + f -5-5 - f50 - f -50 - f05 - f0-5 ) ] + O(h10)

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

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

                                         R01 = x        R04 = f(x,y)
                                         R02 = y        R05 = f "xy = 2f / 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.
 
 

01  LBL "MDV+"
02  STO 01
03  RDN
04  STO 02
05  X<>Y
06  STO 03
07  CLX
08  STO 06
09  XEQ 01
10  42
11  *
12  STO 05
13  XEQ 01
14  6
15  *
16  ST- 05
17  XEQ 01
18  ST+ 05
19  XEQ 01
20  8
21  /
22  ST- 05
23   E3
24  ST* 05
25  XEQ 01
26  8
27  *
28  ST+ 05
29  RCL 02
30  RCL 01
31  XEQ IND 00
32  STO 04
33  73766
34  *
35  RCL 05
36  +
37  50400
38  RCL 03
39  X^2
40  *
41  /
42  STO 05
43  RTN
44  LBL 01
45  RCL 03
46  ST+ 06
47  RCL 02
48  RCL 01
49  RCL 06
50  +
51  XEQ IND 00
52  STO 04
53  RCL 02
54  RCL 01
55  RCL 06
56  -
57  XEQ IND 00
58  ST+ 04
59  RCL 02
60  RCL 06
61  +
62  RCL 01
63  XEQ IND 00
64  ST+ 04
65  RCL 02
66  RCL 06
67  -
68  RCL 01
69  XEQ IND 00
70  ST+ 04
71  RCL 02
72  RCL 01
73  RCL 06
74  ST+ Z
75  +
76  XEQ IND 00
77  ST- 04
78  RCL 02
79  RCL 01
80  RCL 06
81  ST- Z
82  -
83  XEQ IND 00
84  RCL 04
85  -
86  END

 
    ( 134 bytes / SIZE 007 )
 
 

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

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

  LBL "T"
  X^2
  CHS
  E^X
  X<>Y
  LN
  *
  RTN               T  ASTO 00

-If we choose h = 0.1

   0.1   ENTER^
    2     ENTER^
    1     XEQ "MDV+"  >>>>   f "xy = 2f / xy = -0.367879484      ( exact value = -1/e = -0.367879441 )

-In this example, execution time = 29 seconds
-The formula requires 31 evaluations of the function.
 
 

3°) Functions of 3 variables
 

-"dF3" calls "dF+" or "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
 

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

                                         R01 = x        R04 = h              R06 = f 'x = f / x         R09 = f "xx = 2f / x2
                                         R02 = y        R05 = f(x,y,z)      R07 = f 'y = f / y        R10 = f "yy = 2f / y2          R12 & R13: temp
                                         R03 = z                                   R08 = f 'z = f / z         R11 = f "zz = 2f / z2

Flags: /
Subroutines:   "dF+" or "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.
 
 

01  LBL "dF3"   
02  STO 06
03  RDN
04  STO 09
05  RDN
06  STO 12
07  RCL 00
08  STO 13
09  "Z"
10  ASTO 00
11  RDN
12  XEQ "dF+"
13  STO 08
14  X<>Y
15  STO 11
16  "Y"
17  ASTO 00
18  RCL 02
19  RCL 09
20  XEQ "dF+"  
21  STO 07
22  X<>Y
23  STO 10
24  "X"
25  ASTO 00
26  RCL 02
27  RCL 06
28  XEQ "dF+"  
29  STO 06
30  X<>Y
31  X<> 09
32  X<> 02
33  STO 04
34  RCL 12
35  X<> 03
36  STO 05
37  RCL 13
38  STO 00
39  RCL 09
40  RCL 10
41  RCL 11        
42  +
43  +
44  RCL 08
45  RCL 07
46  RCL 06
47  CLA
48  RTN
49  LBL "X"
50  RCL 09
51  RCL 12
52  X<> Z
53  GTO IND 13
54  LBL "Y"
55  RCL 12
56  X<>Y
57  RCL 06
58  GTO IND 13
59  LBL "Z"
60  RCL 09
61  RCL 06
62  GTO IND 13
63  END

 
    ( 108 bytes / SIZE 014 )
 
 

      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

  LBL "T"
  X^2
  CHS
  E^X
  RDN
  X^2
  +
  LN
  R^
  *
  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         and       R09 = f "xx = 2f / x2 =  1.413720635
                                RDN   f 'y  = f / y  =    0.210216822                      R10 = f "yy = 2f / y2 = -0.015015516
                                RDN   f 'z  = f / z   =    0.052554204                      R11 = f "zz2f / z2 = -0.007507659
    RDN    2f / x2 + 2f / y2 + 2f / z2 =     1.409197460
 

-In this example, execution time = 50 seconds
 

4°) An Iterative Method
 

-The following program computes a derivative d with an initial h0-value and one of our formulas above
-Then, the calculation is performed again with  h1 = h0/sqrt(2)  ( lines 12-13 )
-The iteration continues ( with hi+1 = hi/sqrt(2) )  until the difference | di+1 - di |  is no more decreasing or (di) is no more monotonic.   [  di = d(hi) ]

-In other words, the iteration stops when the relation  0 <= ( di+1 - di ) / ( di - di-1 ) < 1  is no more satisfied.
-This relation is given in the PPC ROM user's manual.

-Actually, I've remarked that this test could stop the routine prematurely - for instance in example 1
- so I've replaced this "< 1" above by "< 1.1"  ( line 34 )
-Another value might be better.

-The result is the next to last di and   | di+1 - di |   is used as an error estimate
 

Data Registers:           •  R00 = function name    •  R20 = derivative name     ( Registers R00 & R20 are to be initialized before executing "DRV"  )

    function of 1 variable:     R01 = x    R02 = h                                                                    R19 = 2 , 3 or 4
    function of 2 variables:   R01 = x    R02 = y    R03 = h                                                    R21 = di    R22 = di+1
    function of 3 variables:   R01 = x    R02 = y    R03 = z    R04 = h

Flags:   F10    -      Set flag F02 for a function of 2 variables              ( Clear flags F02 & F03 for a function of 1 variable )
                               Set flag F03 for a function of 3 variables

Subroutines:   A program which computes f(x) [ or f(x,y)  or f(x,y,z) ]  assuming  x is in X-register  [ y in Y-register , z in Z-register ] upon entry
                         and a routine like "dF" , "dF+" , "dF2" , "dF3" , "MDV" ... etc ... which computes the required derivative with the same inputs as "DRV"

-Line 34 may be replaced by another number slighly greater than 1
 
 

01  LBL "DRV"
02  XEQ IND 20
03  STO 21
04  2
05  FS? 02
06  3
07  FS? 03
08  4
09  STO 19
10  SF 10
11  LBL 01
12  2
13  SQRT
14  ST/ IND 19
15  RCL 04
16  RCL 03
17  RCL 02
18  RCL 01
19  XEQ IND 20
20  ENTER^
21  X<> 22
22  FS?C 10       
23  GTO 01
24  VIEW X
25  ST- Y
26  ENTER^
27  X<> 21
28  -
29  X=0?
30  GTO 02
31  /
32  X<0?
33  GTO 02
34  1.1
35  X>Y?
36  GTO 01       
37  LBL 02
38  2
39  SQRT
40  ST* IND 19
41  RCL 22
42  RCL 21
43  -
44  ABS
45  RCL 21       
46  END

 
   ( 77 bytes / SIZE 023 )
 
 

      STACK       INPUTS1       INPUTS2       INPUTS3      OUTPUTS
           T             /             /             h0            /
           Z             /             h0             z            /
           Y             h0             y             y   error estimate
           X             x             x             x     derivative

   INPUTS1  for a function of 1 variable
   INPUTS2  for a function of 2 variables
   INPUTS3  for a function of 3 variables

Example1:     f(x) = exp(-x2)    evaluate  f "(1)

  LBL "T"
  X^2
  CHS
  E^X
  RTN               "T"  ASTO 00

  LBL "D2"
  XEQ "dF+"
  X<>Y
  RTN              "D2"  ASTO 20

-With h0 = 0.3

   CF 02  CF 03
   0.3  ENTER^
    1    XEQ "DRV"  >>>>  the successive approximations are displayed ( except the first one ) and eventually,
                                           f "(1) = 0.735758869
                              X<>Y   error estimate = 70 E-9               the error is actually  13 E-9  only  ( in absolute value )
 

Example2:     f(x) = exp(-x2).Ln(y)    x = 1 , y = 2   evaluate   f "xy = 2f / xy

  LBL "T"
  X^2
  CHS
  E^X
  X<>Y
  LN
  *
  RTN               "T"  ASTO 00     "MDV+"  ASTO 20

-With h0 = 0.3

   SF 02  CF 03
   0.3  ENTER^
    2    ENTER^
    1    XEQ "DRV"  >>>>  the successive approximations are displayed and eventually,
                                           f "xy = 2f / xy  =  -0.367879435
                                    X<>Y   error estimate = 35 E-9                   the error is in fact  7 E-9
 

Example3:     f(x) = exp(-x2).Ln(y2+z)    x = 1 , y = 2 , z = 3   Evaluate    Df = Laplacian ( f ) =   2f / x2 + 2f / y2 + 2f / z2

  LBL "T"
  X^2
  CHS
  E^X
  RDN
  X^2
  +
  LN
  R^
  *
  RTN               "T"  ASTO 00

  LBL "D3"
  XEQ "dF3"
  R^
  RTN               "D3"   ASTO 00

-With h0 = 0.3

   CF 02  SF 03
   0.3  ENTER^
    3    ENTER^
    2    ENTER^
    1    XEQ "DRV"  >>>>  the successive approximations are displayed and eventually,
                                   Df = Laplacian ( f ) =   2f / x2 + 2f / y2 + 2f / z2 =  1.409197493   ( exact value = 1.409197445 )
                                      X<>Y   error estimate = 105 E-9  which is too pessimistic!
 

Notes:

-Do not choose a too small h0-value.
-The results are not always so good, especially for expressions like a Laplacian, because the optimal h-value for a partial derivative with respect to x
  is not necessarily optimal for a partial derivative with respect to y or z!
 

5°) Extrapolation to the Limit
 

-Using high-order formulae requires many evaluations of the function and this may be time consuming for complicated functions, especially with "DRV".
-The following routine employs simpler formulas of order 2 and an extrapolation to the limit
-The method is quite similar to the Romberg integration

-Like "DRV", "DRV2" calculates a derivative d with an initial h0-value:               a00 = d(h0)  = d0
-Then, the calculation is performed again with  h1 = h0/µ                                     a10 = d(h1)
-But here, we use these 2 results to compute an estimation of higher order:         a11 = ( µ2 a10 - a00 ) / ( µ2 - 1 ) = d1

-Thus, we build the following array:          a00 = d(h0)  = d0
                                                               a10 = d(h0/µ)             a11 = ( µ2 a10 - a00 ) / ( µ2 - 1 ) = d1
                                                               a20 = d(h02)           a21 = ( µ2 a20 - a10 ) / ( µ2 - 1 )              a22 = ( µ4 a21 - a11 ) / ( µ4 - 1 )  = d2

                                                               .................................................................. and so on ...............................................................................

-In the following routine, µ = sqrt(2)  ( line 11-12 )

-The program stops when the difference | di+1 - di |  is no more decreasing.
-The result is the next to last di and   | di+1 - di |   is used as an error estimate.
 

Data Registers:           •  R00 = function name    •  R10 = derivative name     ( Registers R00 & R10 are to be initialized before executing "DRV2"  )

    function of 1 variable     R01 = x    R02 = h                                                    R09 = 2 , 3 or 4
    function of 2 variables   R01 = x    R02 = y    R03 = h                                    R11 = di   R12 = | di+1 - di |
    function of 3 variables   R01 = x    R02 = y    R03 = z    R04 = h                    R13 = µ    R14 = counter   R15 = 1 , µ2 , µ4 , µ6 , .....

                       Registers   R05 thru R08 are unused                                            R16 = ai0  R17 = ai1  R18 = ai2 ... etc ...
 

Flags:               Set flag F02 for a function of 2 variables              ( Clear these 2 flags for a function of 1 variable )
                         Set flag F03 for a function of 3 variables

Subroutines:   A program which computes f(x) [ or f(x,y)  or f(x,y,z) ]  assuming  x is in X-register  [ y in Y-register , z in Z-register ] upon entry
                         and a routine which computes the required derivative by a formula of the kind  d = ...... + O(h2)
                         and uses the data in registers R01  R02  R03  R04

-Such formulas and subroutines are listed below.
 
 

01  LBL "DRV2"
02  STO 01
03  RDN
04  STO 02
05  RDN
06  STO 03
07  X<>Y
08  STO 04
09   E99
10  STO 12
11  2
12  SQRT
13  STO 13
14  16
15  STO 14
16  2
17  FS? 02
18  3
19  FS? 03
20  4
21  STO 09
22  XEQ IND 10
23  STO 16
24  LBL 01
25  RCL 13
26  ST/ IND 09
27  RCL 14
28  INT
29   E3
30  /
31  16
32  +
33  STO 14
34  SIGN
35  STO 15
36  XEQ IND 10
37  LBL 02
38  RCL 15
39  RCL 13
40  ENTER^
41  *
42  *
43  STO 15
44  *
45  X<>Y
46  X<> IND 14
47  STO 11
48  -
49  RCL 15
50  1
51  -
52  /
53  ISG 14
54  GTO 02
55  STO IND 14
56  RCL 11
57  -
58  ABS
59  ENTER^
60  X<> 12
61  X>Y?
62  GTO 01
63  CLX
64  RCL 11
65  END

 
    ( 91 bytes )
 
 

      STACK       INPUTS1       INPUTS2       INPUTS3      OUTPUTS
           T             /             /             h0            /
           Z             /             h0             z            /
           Y             h0             y             y   error estimate
           X             x             x             x     derivative

   INPUTS1  for a function of 1 variable
   INPUTS2  for a function of 2 variables
   INPUTS3  for a function of 3 variables

Example:     f(x) = exp(-x2)    evaluate  f '(1)

  LBL "T"
  X^2
  CHS
  E^X
  RTN              "T"  ASTO 00

-For the first derivative, we can use the very simple formula  f '(x) = [ f(x+h) - f(x-h) ] / 2.h

  LBL "DR1"            STO 03               RCL 03           RTN
  RCL 01                  RCL 01               -
  RCL 02                  RCL 02               RCL 02
  -                             +                         ST+ X
  XEQ IND 00         XEQ IND 00       /                                         "DR1"   ASTO 10

-With h0 = 0.1

   CF 02  CF 03
   0.1  ENTER^
    1    XEQ "DRV2"  >>>>  f '(1) = -0.735758876
                                 X<>Y   error estimate = 49 E-9               the error is actually  7 E-9
 

-Do not choose too small h0-values.
-The method is sometimes very disappointing for complex expressions like a Laplacian: Calculate the partial derivatives separately and add these results.

-Here are a few formulas and subroutines needed by "DRV2"
 

Miscellaneous Formulas of order 2
 

1°) Functions of 1 variable

     a) First derivative:     f '(x) = [ f(x+h) - f(x-h) ] / 2.h + O(h2)

  LBL "DR1"             STO 03               RCL 03           RTN
  RCL 01                  RCL 01               -
  RCL 02                  RCL 02               RCL 02
  -                             +                         ST+ X
  XEQ IND 00         XEQ IND 00       /

     b) Second derivative:     f "(x) =  [ f(x+h) - 2.f(x) + f(x-h) ] / h2 + O(h2)

  LBL "DR2"           STO 03                 RCL 01                  +
  RCL 01                RCL 01                 RCL 02                  RCL 02
  RCL 02                XEQ IND 00         -                             X^2
  +                          ST+ X                   XEQ IND 00          /
  XEQ IND 00       ST- 03                   RCL 03                  RTN
 

2°) Functions of 2 variables
 

     a) First derivative:     fx '(x,y) =  f / x = [ f(x+h,y) - f(x-h,y) ] / 2.h + O(h2)

  LBL "DRX"         XEQ IND 00           +                            ST+ X
  RCL 02               STO 04                   XEQ IND 00          /
  RCL 01               RCL 02                   RCL 04                  RTN
  RCL 03               RCL 01                   -
  -                          RCL 03                   RCL 03

     b) First derivative:     fy '(x,y) =  f / y = [ f(x,y+h) - f(x,y-h) ] / 2.h + O(h2)

  LBL "DRY"         XEQ IND 00           RCL 01                  ST+ X
  RCL 02               STO 04                   XEQ IND 00          /
  RCL 03               RCL 02                   RCL 04                  RTN
  -                          RCL 03                   -
  RCL 01               +                             RCL 03

     c) Second derivative:     fxx "(x,y) =  2f / x2 = [ f(x+h,y) - 2.f(x,y) + f(x-h,y) ] / h2 + O(h2)

  LBL "DRXX"      XEQ IND 00           ST+ X                   -                             X^2
  RCL 02               STO 04                   ST- 04                   XEQ IND 00          /
  RCL 01               RCL 02                   RCL 02                  RCL 04                 RTN
  RCL 03               RCL 01                   RCL 01                  +
  +                         XEQ IND 00           RCL 03                  RCL 03

     d) Second derivative:     fyy "(x,y) =  2f / y2 = [ f(x,y+h) - 2.f(x,y) + f(x,y-h) ] / h2 + O(h2)

   LBL "DRYY"      XEQ IND 00           ST+ X                   RCL 01                  X^2
   RCL 02               STO 04                   ST- 04                  XEQ IND 00           /
   RCL 03               RCL 02                   RCL 02                 RCL 04                   RTN
   +                         RCL 01                   RCL 03                 +
   RCL 01               XEQ IND 00          -                            RCL 03

     e) Mixed derivative:     fxy "(x,y) =  2f / xy = [ f(x+h,y+h) - f(x+h,y-h) - f(x-h,y+h) + f(x-h,y-h) ] / 4.h2+ O(h2)

  LBL "DRXY"         +                             RCL 03                   RCL 02           XEQ IND 00        ST- Z                      RCL 03
  RCL 02                  XEQ IND 00          ST+ Z                     RCL 01            ST- 04                  -                             ST+ X
  RCL 01                  STO 04                   -                             RCL 03            RCL 02                XEQ IND 00          X^2
  RCL 03                  RCL 02                   XEQ IND 00         ST+ Z               RCL 01                RCL 04                   /
  ST- Z                     RCL 01                   ST+ 04                   +                      RCL 03                 -                             RTN
 

     f) Laplacian:   fxx "(x,y) +  fyy "(x,y) =  2f / x2 + 2f / y2 = [ f(x+h,y) + f(x,y+h) - 4.f(x,y) + f(x-h,y) + f(x,y-h) ] / h2 + O(h2)

  LBL "LAP2"       XEQ IND 00     RCL 01                  XEQ IND 00         RCL 01                RCL 02                   RCL 04             RTN
  RCL 02              STO 04              XEQ IND 00          4                           RCL 03                RCL 03                   +
  RCL 01              RCL 02              ST+ 04                   *                           -                           -                               RCL 03
  RCL 03              RCL 03              RCL 02                  ST- 04                   XEQ IND 00       RCL 01                    X^2
  +                        +                        RCL 01                  RCL 02                  ST+ 04                XEQ IND 00            /
 

3°) Functions of 3 variables
 

     a) First derivative:     fx '(x,y,z) =  f / x = [ f(x+h,y,z) - f(x-h,y,z) ] / 2.h + O(h2)

  LBL "DR3X"       -                             RCL 01                         -
  RCL 03               XEQ IND 00          RCL 04                        RCL 04
  RCL 02               STO 05                   +                                  ST+ X
  RCL 01               RCL 03                   XEQ IND 00                /
  RCL 04               RCL 02                   RCL 05                        RTN

     b) First derivative:     fy '(x,y,z) =  f / y = [ f(x,y+h,z) - f(x,y-h,z) ] / 2.h + O(h2)

  LBL "DR3Y"       RCL 01                   RCL 04                         -
  RCL 03               XEQ IND 00          +                                  RCL 04
  RCL 02               STO 05                   RCL 01                        ST+ X
  RCL 04               RCL 03                   XEQ IND 00                /
  -                          RCL 02                   RCL 05                        RTN

     c) First derivative:     fz '(x,y,z) =  f / y = [ f(x,y,z+h) - f(x,y,z-h) ] / 2.h + O(h2)

  LBL "DR3Z"       -                              RCL 01                         -
  RCL 03               XEQ IND 00          RCL 04                        RCL 04
  RCL 02               STO 05                   +                                  ST+ X
  RCL 01               RCL 03                   XEQ IND 00                /
  RCL 04               RCL 02                   RCL 05                        RTN

     d) Second derivative:     fxx "(x,y,z) =  2f / x2 = [ f(x+h,y,z) - 2.f(x,y,z) + f(x-h,y,z) ] / h2 + O(h2)

  LBL "DR3XX"       +                              RCL 01                    RCL 02                  RCL 05            RTN
  RCL 03                  XEQ IND 00           XEQ IND 00            RCL 01                  +
  RCL 02                  STO 05                    ST+ X                      RCL 04                  RCL 04
  RCL 01                  RCL 03                    ST- 05                      -                             X^2
  RCL 04                  RCL 02                    RCL 03                    XEQ IND 00          /

     e) Second derivative:     fyy "(x,y,z) =  2f / y2 = [ f(x,y+h,z) - 2.f(x,y,z) + f(x,y-h,z) ] / h2 + O(h2)

  LBL "DR3XX"       RCL 01                    RCL 01                    RCL 02                  RCL 05            RTN
  RCL 03                  XEQ IND 00           XEQ IND 00            RCL 04                  +
  RCL 02                  STO 05                    ST+ X                      -                             RCL 04
  RCL 04                  RCL 03                    ST- 05                     RCL 01                    X^2
  +                            RCL 02                    RCL 03                    XEQ IND 00          /

     f) Second derivative:     fzz "(x,y,z) =  2f / z2 = [ f(x,y,z+h) - 2.f(x,y,z) + f(x,y,z-h) ] / h2 + O(h2)

  LBL "DR3XX"       RCL 01                    RCL 01                    RCL 04                  RCL 05            RTN
  RCL 03                  XEQ IND 00           XEQ IND 00            -                             +
  RCL 04                  STO 05                    ST+ X                      RCL 02                  RCL 04
  +                            RCL 03                    ST- 05                      RCL 01                  X^2
  RCL 02                  RCL 02                    RCL 03                    XEQ IND 00          /

     g) Laplacian:   fxx "(x,y) +  fyy "(x,y) + fzz "(x,y,z) =  2f / x2 + 2f / y2 + 2f / z2
                                                                               = [ f(x+h,y,z) + f(x,y+h,z) + f(x,y,z+h) - 6.f(x,y,z) + f(x-h,y,z) + f(x,y-h,z) + f(x,y,z-h) ] / h2 + O(h2)
 

LBL "LAP3"   +                       RCL 04            RCL 03      XEQ IND 00    XEQ IND 00   RCL 02             ST+ 05     RCL 01            -                       +
RCL 03          XEQ IND 00     +                       RCL 04      ST+ 05             6                      RCL 01             RCL 03    XEQ IND 00    RCL 02           RCL 04
RCL 02          STO 05             RCL 01             +                RCL 03             *                      RCL 04             RCL 02    ST+ 05             RCL 01           X^2
RCL 01          RCL 03             XEQ IND 00     RCL 02     RCL 02             ST- 05             -                        RCL 04     RCL 03           XEQ IND 00    /
RCL 04          RCL 02             ST+ 05              RCL 01     RCL 01             RCL 03            XEQ IND 00    -                RCL 04            RCL 05           RTN
 

6°) Functions of n variables  ( n < 11 )
 

-The following routines generalize the programs above.
-They are limited to functions of 10 ( or less ) variables.
-If you want to study functions of more than 10 variables, simply shift registers R11, R12, ... in these listings.
 

     a) First Derivatives
 

Data Registers:   •  R00 = Function name
                              •  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn           ( These ( n+1 ) registers are to be initialized before executing "FD"  )

                                           R11 = h                 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
 
 

01  LBL "FD"
02  STO 13
03  X<>Y
04  STO 11
05  CLX
06  STO 14
07  RCL IND Y 
08  STO 15
09  XEQ 01
10  84
11  *
12  STO 12
13  XEQ 01
14  24
15  *
16  ST- 12
17  XEQ 01       
18  6
19  *
20  ST+ 12
21  XEQ 01
22  ST- 12
23  25
24  ST* 12
25  XEQ 01
26  ST+ X
27  ST+ 12
28  RCL 11
29  2520
30  *
31  ST/ 12
32  RCL 15
33  STO IND 13
34  RCL 12
35  RTN
36  LBL 01
37  RCL 11
38  ST+ 14
39  RCL 15
40  RCL 14
41  -
42  STO IND 13
43  XEQ IND 00
44  STO 16
45  RCL 14
46  RCL 15
47  +
48  STO IND 13
49  XEQ IND 00
50  RCL 16
51  -
52  END

 
   ( 91 bytes / SIZE 017 )
 
 

      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 )

  0.1   ENTER^
    2    XEQ "FD"   >>>>   f 'x2 = f / x2 = 0.367879440     ( The last decimal should be a 1 )

  0.1   ENTER^
    3    XEQ "FD"   >>>>   f 'x3 = f / x3 = 0.183939720     ( The last decimal should be a 1 )

-This program uses the formula of order 10 given in paragraph 1°) b)
 

     b) Second Derivatives
 

Data Registers:   •  R00 = Function name
                              •  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn      ( These ( n+1 ) registers are to be initialized before executing "SD"  )

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

01  LBL "SD"
02  CF 10
03  X=Y?
04  SF 10
05  STO 13
06  RDN
07  STO 14
08  X<>Y
09  STO 11
10  CLX
11  STO 15
12  RCL IND T  
13  STO 16
14  RCL IND Z
15  STO 17
16  XEQ 01
17  42
18  *
19  STO 12
20  XEQ 01
21  6
22  *
23  ST- 12
24  XEQ 01
25  ST+ 12
26  XEQ 01
27  8
28  /
29  ST- 12
30   E3
31  ST* 12
32  XEQ 01
33  8
34  *
35  ST+ 12
36  RCL 16
37  STO IND 13
38  XEQ IND 00
39  73766
40  *
41  ST- 12
42  25200
43  FC? 10
44  ST+ X
45  FC?C 10
46  CHS
47  RCL 11
48  X^2
49  *
50  ST/ 12
51  RCL 12        
52  RTN
53  LBL 01
54  RCL 11
55  ST+ 15
56  RCL 15
57  RCL 16
58  +
59  STO IND 13
60  XEQ IND 00
61  STO 18
62  RCL 16
63  RCL 15
64  -
65  STO IND 13
66  XEQ IND 00
67  ST+ 18
68  RCL 18
69  FS? 10
70  RTN
71  RCL 17
72  RCL 15
73  -
74  STO IND 14
75  XEQ IND 00
76  ST- 18
77  RCL 16
78  STO IND 13
79  XEQ IND 00
80  ST+ 18
81  RCL 15
82  RCL 17
83  +
84  STO IND 14
85  XEQ IND 00
86  ST+ 18
87  RCL 15
88  ST+ IND 13
89  XEQ IND 00
90  ST- 18
91  RCL 17
92  STO IND 14
93  RCL 18
94  END

 
   ( 169 bytes / SIZE 019 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             h             /
           Y             j             /
           X             i  f "xixj = 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 )

 0.1   ENTER^
   1    ENTER^
   2    XEQ "SD"  >>>>   f "xy = 2f / xy = -0.735758889  ( the last decimal should be a 2 )

Notes:

-If  i = j  this program uses the formula of paragraph 1°) b)
-If  i # j  -------------------------------------------  2°) b-2)
-Both formulas are of order 10 but the first one requires less evaluations of the function than the second one.

-The following routine may be used to compute the Laplacian
-Replace R20 by another register ( like synthetic register M ) or R10 ( if n < 10 )
  if you want to use "LAPL" as a subroutine of the following program "DV"
 
 

 01  LBL "LAPL"
 02  STO 19
 03  X<>Y
 04  STO 11
 05  CLX
 06  STO 20
 07  LBL 01
 08  RCL 11
 09  RCL 19
 10  ENTER^
 11  XEQ "SD"
 12  ST+ 20
 13  DSE 19
 14  GTO 01
 15  RCL 20
 16  END

 
 ( 35 bytes / SIZE 021 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             h             /
           X             n      Laplacian(f)

  where n is the number of variables

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

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

-With h = 0.1

   0.1   ENTER^
    4     XEQ "LAPL"  >>>>    Laplacian(f) =  2f / x2 + 2f / y2 + 2f / z2 + 2f / t2 = 0.673013929  ( the last 2 decimals should be 32 )
 

     c) Iterative Method
 

-This program employs the method described in §4
 

Data Registers:   •  R00 = Function name

                              •  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn      ( These ( n+1 ) registers and R20 are to be initialized before executing "DV"  )

                 R11 thru R16 ( or R18 ) are used by "FD" or "SD"                 •  R20 = "FD" or "SD" or another similar program
                                                                                                                  R21 = dk   R22 = dk+1
Flags: F09 - F10
                     F02:  Clear flag F02 to compute a first derivative
                                 Set flag F02 to compute a second derivative

Subroutines:  "FD" or "SD"
                         and a program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn

-If you plan to use "DV" with "FD" and "SD" only,
 add "FD"  FS? 02  "SD"  ASTO 20  after line 01 and register R20 is no more to be initialized.
 
 

01  LBL "DV"
02  XEQ IND 20
03  STO 21
04  SF 09
05  LBL 01
06  RCL 11
07  2
08  SQRT
09  /
10  FS? 02
11  RCL 14
12  RCL 13
13  XEQ IND 20
14  ENTER^
15  X<> 22
16  FS?C 09
17  GTO 01
18  VIEW X       
19  ST- Y
20  ENTER^
21  X<> 21
22  -
23  X=0?
24  GTO 02
25  /
26  X<0?
27  GTO 02       
28  1.1
29  X>Y?
30  GTO 01
31  LBL 02
32  2
33  SQRT
34  ST* 11
35  RCL 21
36  RCL 22       
37  -
38  ABS
39  RCL 21
40  END

 
   ( 67 bytes / SIZE 023 )
 
 

      STACK       INPUTS1       INPUTS2      OUTPUTS
           Z             /             h0            /
           Y             h0             j   error estimate
           X             i             i     derivative

  INPUTS1  for a first derivative - in this case, derivative =  f 'xi = f / xi
  INPUTS2  for a second derivative - in this case, derivative = f "xixj = 2f / xixj

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  h0 = 0.2

  CF 02  "FD"  ASTO 20
    0.2   ENTER^
     4     XEQ "DV"  >>>> the successive approximations are displayed ( except the first one ) and eventually,
                                         f 't = f / t = -0.254994598 ( the last decimal should be a 7 )
                                        X<>Y  errror estimate = E-9

  SF 02  "SD" ASTO 20
    0.2    ENTER^
     1      ENTER^
     3         R/S      >>>> the successive approximations are displayed and finally,
                                      f "xz = 2f / xz = -0.367879443 ( the last decimal should be a 1 )
                                      X<>Y  errror estimate = 80 E-9 which is too pessimistic!
 

     d) Extrapolation to the Limit
 

-This program employs the method described in §5
 

Data Registers:   •  R00 = Function name
                              •  R01 = x1 ,  •  R02 = x2 , .......... ,  •  Rnn = xn      ( These ( n+1 ) registers and R20 are to be initialized before executing "DV2"  )

                 R11 thru R14 ( or R16 ) are used by "FD2" or "SD2"             •  R20 = "FD2" or "SD2" or the name of another similar program

                                                                                                                  R21 = dk                      R24 = 26.0....
                                                                                                                  R22 =  | dk+1 - dk |       R25 = 1 , µ2 , µ4 , µ6 , .....
                                                                                                                  R23 = µ = sqrt(2)         R26 = ai0  R27 = ai1  R28 = ai2 ... etc ...
Flags:     F02:  Clear flag F02 to compute a first derivative
                             Set flag F02 to compute a second derivative

Subroutines:  "FD2" or "SD2" which are listed at the bottom of this page
                         and a program which computes  f(x1, ..... ,xn)  assuming  x1 is in R01 , ............ , xn is in Rnn
 
 

01  LBL "DV2"
02  XEQ IND 20
03  STO 26
04   E99
05  STO 22
06  2
07  SQRT
08  STO 23
09  26
10  STO 24
11  LBL 01
12  RCL 24
13  INT
14   E3
15  /
16  26
17  +
18  STO 24       
19  SIGN
20  STO 25
21  RCL 11
22  RCL 23
23  /
24  FS? 02
25  RCL 14
26  RCL 13
27  XEQ IND 20
28  LBL 02
29  RCL 25
30  RCL 23
31  ENTER^
32  *
33  *
34  STO 25
35  *
36  X<>Y
37  X<> IND 24
38  STO 21
39  -
40  RCL 25
41  1
42  -
43  /
44  ISG 24
45  GTO 02
46  STO IND 24
47  RCL 21
48  -
49  ABS
50  ENTER^
51  X<> 22
52  X>Y?
53  GTO 01       
54  CLX
55  RCL 21
56  END

 
   ( 93 bytes / SIZE 027+ ? )
 
 

      STACK       INPUTS1       INPUTS2      OUTPUTS
           Z             /             h0            /
           Y             h0             j   error estimate
           X             i             i     derivative

  INPUTS1  for a first derivative - in this case, derivative =  f 'xi = f / xi
  INPUTS2  for a second derivative - in this case, derivative = f "xixj = 2f / xixj

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  h0 = 0.1

  CF 02  "FD2"  ASTO 20   ( 1st derivative: load the subroutine "FD2" listed hereafter )

    0.1   ENTER^
     4     XEQ "DV2"  >>>>   f 't = f / t = -0.254994598 ( the last decimal should be a 7 )
                                            X<>Y  errror estimate =  3 E-9

  SF 02  "SD2" ASTO 20   ( 2nd derivative: load the subroutine "SD2" listed hereafter )

   0.1    ENTER^
     1      ENTER^
               R/S      >>>> f "xx = 2f / x2 = -0.509989079 ( the last 3 decimals should be 195 )
                                      X<>Y  errror estimate = 33 E-9  which is too optimistic!

    0.1    ENTER^
     1      ENTER^
     3         R/S      >>>>  f "xz = 2f / xz = -0.367879423 ( the last 2 decimals should be 41 )
                                      X<>Y  errror estimate = 61 E-9
 

     a) First derivatives

  LBL "FD2"     RCL IND Y     XEQ IND 00     +                            -                           RCL 11         RTN
  STO 13          STO 14            STO 12             STO IND 13          RCL 14                ST+ X
  X<>Y             X<>Y              RCL 11             XEQ IND 00         STO IND 13         /
  STO 11          ST- IND 13     RCL 14             RCL 12                  CLX                     STO 12

     b) Second derivatives

  LBL "SD2"     STO 11             ST+ IND 13         XEQ IND 00          RCL 11                  4                             RCL 11                   XEQ IND 00         /
  CF 10            RCL IND T       ST+ IND 14         ST- 12                    RCL 16                  ST/ 12                     ST+ IND 13            ST+ 12                 STO 12
  X=Y?            STO 15              XEQ IND 00        RCL 15                  +                            GTO 02                   XEQ IND 00          LBL 02                 END
  SF 10            RCL IND Z        STO 12                RCL 11                  STO IND 14          LBL 01                    ST+ 12                   RCL 15
  STO 13         STO 16              RCL 16                -                             XEQ IND 00          XEQ IND 00           RCL 15                  STO IND 13
  RDN             FS?C 10             RCL 11                STO IND 13          ST- 12                    ST+ X                     RCL 11                  RCL 12
  STO 14         GTO 01              -                           XEQ IND 00         RCL 16                   CHS                        -                             RCL 11
  X<>Y            RCL 11             STO IND 14         ST+ 12                  STO IND 14           STO 12                   STO IND 13          X^2
 

-These programs use formulas of order 2 and may be replaced by similar subroutines,
  provided the formulas are of the kind  . . . + O(h2)
 

References:

[1]  F.Scheid - "Numerical Analysis" - Mc Graw Hill - ISBN  07-055197-9
[2]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[3]  Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4