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(h⁴)
d²f/dx² = (1/12h²).[ - f(x-2h) + 16.f(x-h) - 30.f(x) +16.f(x+h) - f(x+2h) ] + O(h⁴)

-These formulas ( without + O(h⁴) ) 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(-x²) 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.( f₁ - f_-1 ) - 600.( f₂ - f_-2 ) + 150.( f₃ - f_-3 ) - 25.( f₄ - f_-4 ) + 2.( f₅ - f_-5 ) ] + O(h¹⁰)
d²f/dx² = (1/25200.h²).[ -73766 f₀ + 42000.( f₁ + f_-1 ) - 6000.( f₂ + f_-2 ) + 1000.( f₃ + f_-3 ) - 125.( f₄ + f_-4 ) + 8.( f₅ + f_-5 ) ] + O(h¹⁰)

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

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

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

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 = ¶²f / ¶x²
R02 = y R04= f(x,y) R06 = f '_y = ¶f / ¶y R08 = f "_yy = ¶²f / ¶y² 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 = ¶²f / ¶y²
Z	h	f "_xx = ¶²f / ¶x²
Y	y	f '_y = ¶f / ¶y
X	x	f '_x = ¶f / ¶x

Example: f(x) = exp(-x²).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 = ¶²f / ¶x² = 0.509989206      the exact result is     0.509989195
                               RDN f "_yy = ¶²f / ¶y² = -0.091969921      the exact result is    -0.091969860

-Whence Laplacian ( f ) = ¶²f / ¶x² + ¶²f / ¶y² = 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 = ¶²f / ¶x¶y = (1/24.h²).[ 30 f₀₀ + 16.( f₁₁ + f_-1-1 - f₁₀ - f_-10 - f₀₁ - f_0-1 ) + ( - f₂₂ - f_-2-2 + f₂₀ + f_-20 + f₀₂ + f_0-2 ) ] + O(h⁴)

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

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

Example: f(x) = exp(-x²).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 = ¶²f / ¶x¶y = -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 = ¶²f / ¶x¶y = (1/50400.h²).[ 73766 f₀₀ + 42000.( f₁₁ + f_-1-1 - f₁₀ - f_-10 - f₀₁ - f_0-1 ) + 6000.( - f₂₂ - f_-2-2 + f₂₀ + f_-20 + f₀₂ + f_0-2 )
                                                         + 1000.( f₃₃ + f_-3-3 - f₃₀ - f_-30 - f₀₃ - f_0-3 ) + 125.( - f₄₄ - f_-4-4 + f₄₀ + f_-40 + f₀₄ + f_0-4 )
                                                               + 8.( f₅₅ + f_-5-5 - f₅₀ - f_-50 - f₀₅ - f_0-5 ) ] + O(h¹⁰)

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

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

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

Example: f(x) = exp(-x²).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 = ¶²f / ¶x¶y = -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 = ¶²f / ¶x² ; f "_yy = ¶²f / ¶y² ; f "_zz = ¶²f / ¶z²

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 = ¶²f / ¶x²
                                         R02 = y        R05= f(x,y,z)      R07 = f '_y = ¶f / ¶y        R10 = f "_yy = ¶²f / ¶y²          R12 & R13: temp
                                         R03 = z                                   R08 = f '_z = ¶f / ¶z         R11 = f "_zz = ¶²f / ¶z²

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

Example: f(x) = exp(-x²).Ln(y²+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 = ¶²f / ¶x² = 1.413720635
                                RDN   f '_y = ¶f / ¶y =    0.210216822                      R10 = f "_yy = ¶²f / ¶y² = -0.015015516
                                RDN   f '_z = ¶f / ¶z   =    0.052554204                      R11 = f "_zz = ¶²f / ¶z² = -0.007507659
    RDN    ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z² =     1.409197460

-In this example, execution time = 50 seconds

4°) An Iterative Method

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

-In other words, the iteration stops when the relation 0 <= ( d_i+1 - d_i ) / ( d_i - d_i-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 d_i and | d_i+1 - d_i | 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 = d_i    R22 = d_i+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	/	/	h₀	/
Z	/	h₀	z	/
Y	h₀	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(-x²) 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 h₀ = 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(-x²).Ln(y) x = 1 , y = 2 evaluate f "_xy = ¶²f / ¶x¶y

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

-With h₀ = 0.3

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

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

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 h₀ = 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 ) =   ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z² = 1.409197493   ( exact value = 1.409197445 )
                                      X<>Y   error estimate = 105 E-9 which is too pessimistic!

Notes:

-Do not choose a too small h₀-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 h₀-value:               a₀₀ = d(h₀) = d₀
-Then, the calculation is performed again with h₁ = h₀/µ                                     a₁₀ = d(h₁)
-But here, we use these 2 results to compute an estimation of higher order:         a₁₁ = ( µ² a₁₀ - a₀₀ ) / ( µ² - 1 ) = d₁

-Thus, we build the following array:          a₀₀ = d(h₀) = d₀
                                                               a₁₀ = d(h₀/µ)             a₁₁ = ( µ² a₁₀ - a₀₀ ) / ( µ² - 1 ) = d₁
                                                               a₂₀ = d(h₀/µ²)           a₂₁ = ( µ² a₂₀ - a₁₀ ) / ( µ² - 1 )              a₂₂ = ( µ⁴ a₂₁ - a₁₁ ) / ( µ⁴ - 1 ) = d₂

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

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

-The program stops when the difference | d_i+1 - d_i | is no more decreasing.
-The result is the next to last d_i and | d_i+1 - d_i | 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 = d_i   R12 = | d_i+1 - d_i |
    function of 3 variables   R01 = x    R02 = y    R03 = z    R04 = h                    R13 = µ    R14 = counter   R15 = 1 , µ² , µ⁴ , µ⁶ , .....

Registers R05 thru R08 are unused R16 = a_i0 R17 = a_i1 R18 = a_i2 ... 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(h²)
                         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	/	/	h₀	/
Z	/	h₀	z	/
Y	h₀	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(-x²) 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 h₀ = 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 h₀-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(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       /

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

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: f_x '(x,y) = ¶f / ¶x = [ f(x+h,y) - f(x-h,y) ] / 2.h + O(h²)

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: f_y '(x,y) = ¶f / ¶y = [ f(x,y+h) - f(x,y-h) ] / 2.h + O(h²)

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: f_xx "(x,y) = ¶²f / ¶x² = [ f(x+h,y) - 2.f(x,y) + f(x-h,y) ] / h² + O(h²)

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: f_yy "(x,y) = ¶²f / ¶y² = [ f(x,y+h) - 2.f(x,y) + f(x,y-h) ] / h² + O(h²)

   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: f_xy "(x,y) = ¶²f / ¶x¶y = [ f(x+h,y+h) - f(x+h,y-h) - f(x-h,y+h) + f(x-h,y-h) ] / 4.h²+ O(h²)

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: f_xx "(x,y) + f_yy "(x,y) = ¶²f / ¶x² + ¶²f / ¶y² = [ f(x+h,y) + f(x,y+h) - 4.f(x,y) + f(x-h,y) + f(x,y-h) ] / h² + O(h²)

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: f_x '(x,y,z) = ¶f / ¶x = [ f(x+h,y,z) - f(x-h,y,z) ] / 2.h + O(h²)

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: f_y '(x,y,z) = ¶f / ¶y = [ f(x,y+h,z) - f(x,y-h,z) ] / 2.h + O(h²)

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: f_z '(x,y,z) = ¶f / ¶y = [ f(x,y,z+h) - f(x,y,z-h) ] / 2.h + O(h²)

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: f_xx "(x,y,z) = ¶²f / ¶x² = [ f(x+h,y,z) - 2.f(x,y,z) + f(x-h,y,z) ] / h² + O(h²)

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: f_yy "(x,y,z) = ¶²f / ¶y² = [ f(x,y+h,z) - 2.f(x,y,z) + f(x,y-h,z) ] / h² + O(h²)

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: f_zz "(x,y,z) = ¶²f / ¶z² = [ f(x,y,z+h) - 2.f(x,y,z) + f(x,y,z-h) ] / h² + O(h²)

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: f_xx "(x,y) + f_yy "(x,y) + f_zz "(x,y,z) = ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z²
= [ 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) ] / h² + O(h²)

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 = x₁ , • R02 = x₂ , .......... , • Rnn = x_n ( These ( n+1 ) registers are to be initialized before executing "FD" )

                                           R11 = h                 R13 = i                        R15 = x_i
                                           R12 = ¶f / ¶x_i       R14 = 0 , h , 2h , ....     R16: temp
Flags: /
Subroutine:   A program which computes f(x₁, ..... ,x_n) assuming x₁ is in R01 , ............ , x_n is in Rnn

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 / ¶x_i

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

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

-With h = 0.1

0.1 ENTER^
1 XEQ "FD" >>>> f '_x1 = ¶f / ¶x₁ = -0.509989198 ( The last decimal should be a 5 )

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

0.1 ENTER^
3 XEQ "FD" >>>> f '_x3 = ¶f / ¶x₃ = 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 = x₁ , • R02 = x₂ , .......... , • Rnn = x_n ( These ( n+1 ) registers are to be initialized before executing "SD" )

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

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

-With h = 0.1

0.1 ENTER^
1 ENTER^
XEQ "SD" >>>> f "_xx = ¶²f / ¶x² = 0.509989206 ( exact derivative = 0.509989195 )

0.1   ENTER^
   1    ENTER^
   2    XEQ "SD" >>>>   f "_xy = ¶²f / ¶x¶y = -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(-x².t).Ln(y²+z) x = 1 , y = 1 , z = 1 , t = 1

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

-With h = 0.1

0.1 ENTER^
4 XEQ "LAPL" >>>> Laplacian(f) = ¶²f / ¶x² + ¶²f / ¶y² + ¶²f / ¶z² + ¶²f / ¶t² = 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 = x₁ , • R02 = x₂ , .......... , • Rnn = x_n ( 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 = d_k   R22 = d_k+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(x₁, ..... ,x_n) assuming x₁ is in R01 , ............ , x_n 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	/	h₀	/
Y	h₀	j	error estimate
X	i	i	derivative

INPUTS1 for a first derivative - in this case, derivative = f '_xi = ¶f / ¶x_i
INPUTS2 for a second derivative - in this case, derivative = f "_xixj = ¶²f / ¶x_i¶x_j

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

-With h₀ = 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 = ¶²f / ¶x¶z = -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 = x₁ , • R02 = x₂ , .......... , • Rnn = x_n ( 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 = d_k                      R24 = 26.0....
                                                                                                                  R22 = | d_k+1 - d_k |       R25 = 1 , µ² , µ⁴ , µ⁶ , .....
                                                                                                                  R23 = µ = sqrt(2)         R26 = a_i0 R27 = a_i1 R28 = a_i2 ... 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(x₁, ..... ,x_n) assuming x₁ is in R01 , ............ , x_n 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	/	h₀	/
Y	h₀	j	error estimate
X	i	i	derivative

INPUTS1 for a first derivative - in this case, derivative = f '_xi = ¶f / ¶x_i
INPUTS2 for a second derivative - in this case, derivative = f "_xixj = ¶²f / ¶x_i¶x_j

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

-With h₀ = 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 = ¶²f / ¶x² = -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 = ¶²f / ¶x¶z = -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(h²)

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

hp41programs

Numerical Differentiation for the HP-41