Derivative on HP-41

This routine will approximate the first derivative of a function at a point in one of two ways. A quick four point polynomial estimate may be made using a step size that is provided by the user, or a more precise adaptive procedure may be use which automatically searches for the optimal step size. The code originates PPC ROM Module [1] released by the international group People Programming Computers on December 1981.

Instructions

In the adaptive routine setting a flag allows the user to view convergence of the optimal step size. The adaptive procedure delivers an error estimate along with the derivative, but 6.5 decimal digits is all that can be trusted in any case.

Both versions of FD sample on only one side of the evaluation point so that FD may be used for one-sided derivatives. Discontinuities or singularities may be avoided by selecting an increment with the appropriate sign. The routine may also be used to compute partial derivatives.

1. First select a SIZE which must be a minimum of 018 depending on the number of actual additional registers required by the function.
2. The function must be programmed as a subroutine in program memory with a global label name of six of less characters. The input to this subroutine, namely X, will come from a data register pointed to by R11. The output, namely ƒ(x), should be left in the X-register. ƒ(x) should not modify registers R11 through R17, flag F09 or flag F10.
3. Store the global function name from step 2 in register R10.
4. Store the register to hold X in the ƒ(x) subroutine in R11. This register may normally be any register other than any of R10 though R17.
5. If the calculation of ƒ'(x) is desired, store the initial X-value “a” in the register named in step 4.
6. Store an initial step size in register R12.
7. Set or reset F09 to select for receptively the quick approximation, or the derivative procedure. Note that the initial step size should be large enough to the algorithm to iterate at least three times before terminating.
8. If the adaptive procedure was selected in step 7 then the convergence may be viewed as an option (set F10 to do so).
9. XEQ “FD”. An estimate of the first derivative will be left in the X-register. In addition, if the adaptive procedure was used an error estimate will appear in the Y-register.

Examples

1. Approximate ƒ'(x) where ƒ(x) = 3.x³ – 4.x² + 5.x + 6
   :
   Establish SIZE018., and place the subroutine listed below in program memory.

   01	LBL "EX1"
   02	RCL IND 11
   03	ENTER^
   04	ENTER^
   05	ENTER^
   06	3
   07	*
   08	4
   09	-
   10	*
   11	5
   12	+
   13	*
   14	6
   16	RTN

   Store the function name “EX1” in register R10.
   Set flag F09 to select the quick approximation method for this example.
   Store 18 in R11 (pointer to register which holds the X value)
   Store .1 in R12 (step size).
   Store 2 in R18 (X value)XEQ “FD” the gives the answer 25, this is exact because the polynomial formula is exact for polynomials of degree three or less.
2. Find the gradient of ƒ(x,y) = (x + ln(y))² at the point P(2,1)
   :
   The gradient is simply a vector whose components are the partial derivatives of ƒ(x,y), so this example will require the calculation of two partial derivatives.
   :
   Create ƒ(x,y) in program memory with x in register R20 and yin register R21.

   01	LBL "FXY2"
   02	RCL 20
   03	RCL 21
   04	LN
   05	+
   06	X^2
   07	RTN

   The following auxiliary program will be used to setup and execute FD.

   08	LBL "GRAD"
   09	"FXY2"
   10	ASTO 10
   11	20
   12	STO 11
   13	.1
   14	STO 12
   15	2
   16	STO 20
   17	1
   18	STO 21
   19	CF 09
   20	SF 10
   21	XEQ "FD"
   22	STO 22
   23	ISG 11
   24	ABS
   25	XEQ "FD"
   26	STO 23
   27	RTN

   After executing GRAD the gradient will be in memory with dƒ/dx in R22 (4.000000000) and dƒ/dy in R23 (3.999999379).

For more examples refer to [1].

Method

The quick 4-point polynomial approximation uses the following formula:

ƒ'(x) = [-11.ƒ(x) + 18.ƒ(x+h) – 9.ƒ(x+2.h) + 2.ƒ(x+3.h)] / 6.h

where h is the current step size.

The above formula may be found in the Handbook of Mathematical Formulas, table 25.2. The formula samples the user’s function on only one size of the point where the derivative is being evaluated; on the right if h>0 and on the left if h<0. Thus a discontinuity or singularity may be avoided by selecting an increment with the appropriate sign. The formula will be exact for polynomials of degree 3 or less, and the average of estimates will be exact for polynomials of degree 4 or less. For most functions, however, the error of the polynomial is of the order of h to the 4^th power.

The adaptive routine works by evaluating the polynomial formula for successively smaller h values. The sequence of estimates D_i should be monotonic and the sequence | D_i+1 – D_i | should be monotonic and decreasing. If either monotonicity condition is violated, numerical truncation is causing error in D_i+1 and D_i is delivered to the user as the final output. | D_i+1 – D_i | is placed in the stack as an estimate of the error. Note however that 6 1/2 digits is the most you can depend on regardless of the error estimate. The final h value is left in R12 so that any subsequent derivative evaluations at nearby points may be made by the quick polynomial formula using the nearby optimum increment value h.

Listing

Available through the repository

• Requires
  • X-Functions module on the HP-41cv
• Available as
  • source code
  • raw code for the V41 emulator
  • bar code for the HP Wand (HP82153A)
• Note
  • For line by line comments refer to [1,2,3].

01	LBL "FD"
02	FS? 09		; quick polynomial?
03	GTO 03

04	17		; store display mode in R17 (PPC SD-function)
05	SIGN		;  for details refer to [3]
06	RDN
07	RCL d
08	STO M
09	>"\00\00"
10	X<> M
11	"*"
12	X<> M
13	STO N
14	ASTO IND L
15	RDN

16	SCI 1		; change display mode
17	2 E-3		; R14 is the loop counter
18	STO 14

19	LBL 00		; both part of init and part of main loop
20	RCL 12		;   see [1]
21	.7
22	*
23	RDN
24	STO 12
25	XEQ 08
26	ENTER^
27	X<> 16
28	-
29	ENTER^
30	FS? 10
31	VIEW X
32	X<> 15
33	ISG 14
34	GTO 00

35	LASTX		; main loop cont'd
36	RDN		;   see [1]
37	X=0?
38	GTO 02
39	/		; test monotonicity
40	 E		;   see [1]
41	X<>Y
42	X<0?
43	GTO 01
44	X<Y?
45	GTO 00

46	LBL 01		; setup stack contents
47	X<> L		;  see [1]
48	LBL 02
49	R^
50	.7
51	ST/ 12		; restore the previous h value
52	CLX 		;   without lifting the stack

53	17		; recall display mode from R17 "PPC RD-function"
54	SIGN		;   for details refer to [2]
55	ARCL IND L
56	RDN
57	RCL d
58	STO N
59	>"**"
60	X<> O
61	STO N
62	>"*****"
63	X<> N
64	STO d
65	RDN
66	CLA
67	RTN		; return

68	LBL 03		; 4-point polynomial estimate
69	.
70	STO 13
71	XEQ IND 10
72	11
73	XEQ 09
74	-18
75	XEQ 09
76	9
77	XEQ 09
78	ST+ X
79	RCL 13
80	-
81	RCL 12
82	3
83	*
84	ST- IND 11
85	ST+ X
86	/
87	ENTER^
88	RTN

89	LBL 09
90	*
91	ST+ 13
92	RCL 12
93	ST+ IND 11
94	GTO IND 10
95	END

References