The Museum of HP Calculators
The Museum of HP Calculators
This program is Copyright © 2007 by Jean-Marc Baillard and is used here by permission.
This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.
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+"
or 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+"
or 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
-0.509989195
RDN f 'y = ¶f
/ ¶y =
0.183939720 the exact results
are 0.183939721
RDN f "xx =
¶2f /
¶x2 =
0.509989206
0.509989195
RDN f "yy =
¶2f /
¶y2 =
-0.091969921
-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 / ¶x¶y = (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 /
¶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 = ¶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 /
¶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 =
¶2f /
¶x¶y =
(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 /
¶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 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 /
¶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 =
¶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+"
or 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+"
or 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+"
or 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 |
where 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 "zz =
¶2f /
¶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"
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
or another number slighly greater than 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 / ¶x¶y
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 /
¶x¶y
= -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(h0/µ2)
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 / ¶x¶y = [ 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 /
¶xi¶xj
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 /
¶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
01 LBL "LAPL"
02 STO 19
03 X<>Y
04 STO 11
05 CLX
06 STO
20
Choose another register ( like synthetic register M ) or R10 ( if n <
10 )
07 LBL
01
if you want to use "LAPL" as a subroutine of the following program "DV"
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
01 LBL "DV"
02 XEQ IND
20
If you plan to use "DV" with "FD" and "SD" only,
03 STO
21
add "FD" FS? 02 "SD" ASTO 20 after line 01 and register
R20 is no more to be initialized.
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 /
¶xi¶xj
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 /
¶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 = 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 /
¶xi¶xj
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 /
¶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(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
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall