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 / ¶x¶y

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 / ¶x¶y

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 / ¶xi¶xj

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

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

-With  h = 0.1

0.1  ENTER^
1    ENTER^
XEQ "SD"  >>>>   f "xx = 2f / x2 = 0.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:

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