+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (35S) 1st order Derivative (at a point) (/thread-15882.html)
(35S) 1st order Derivative (at a point) - trojdor - 11-11-2020 11:11 AM
First order Numerical (explicit) Differentiation program for HP 35s
The companion program, 2nd order Derivative, can be found at https://www.hpmuseum.org/forum/thread-15875.html
This program is by Michael Lyons, and was inspired by the work by Alexis Zuleta.
OVERVIEW
This program is an alternative to the first order explicit derivative program found at https://www.hpmuseum.org/software/35derivp.htm
Zuleta uses this simple (forward distancing) approximation of the first derivative:
f'(x) ~ (f(x+h) - f(x)) / h (where h>0 : specifically, 1E-6 in his program)
D.Levy states that a centered differencing formula is more accurate, and may have practical advantages:
(It is generally the way some other calculators, such as TI-36x estimate f' at a point.)
f'(x) ~ (f(x+h)-f(x-h)) / 2h (D.Levy EQ 5.4 Introduction to Numerical Analysis)
I also chose to use an h value of 1E-6, which seems to work well.
This program is compact, uses a minimum of labels, and keeps the stack clean.
(The F label, shared by other programs, the D label and 1 variable 'D'.)
LISTING
F001 LBL F
F002 --> enter your function / eq here <-- (If F(x) has any trig func., be sure to set calc mode to Radians.)
F003 RTN
Be sure to have that RTN statement after the eq/function in F!
D001 LBL D
D002 INPUT D
D003 1E-6
D004 +
D005 XEQ D014
D006 RCL D
D007 1E-6
D008 -
D009 XEQ D014
D010 -
D011 2E-6 <-- note that is 2*h, not h (that's why it's different from step 3 and 7)
D012 /
D013 RTN
D014 STO X
D015 XEQ F001
D016 x<>y
D017 CLx
D018 R(down)
D019 RTN
EXAMPLE
If f(x)= x^4 + 5x^3 -21x^2
Calcuate f'(x) at x=4
Steps:
1. GTO F [ENTER]
2. -> [PRGM]
3. Go to line 2 of the F program and enter/edit the EQN to be x^4+5*x^3-21*x^2
4. Make sure there's only one EQN, and that you haven't accidentally deleted the RTN stmt.
5. Press [C/on] to exit program mode.
6. XEQ D [ENTER]
7. Displays D?
8. Enter the point to evaluate, in this case '4', then [R/S].
9. See RUNNING, then the result. In this case, 328.
Let's manually solve to confirm:
f'(x)=x^4 + 5x^3 -21x^2 | x=4
=4x^3 + 15x^2 -42x | x=4
=256 + 240 - 168
=328 <-- yes, it matches the calculator
Note:
This (and other) programs were written by me for the HP 35s to assist engineers in preparing/passing the NCEES FE/PE licensing exams.
The 35s is the only current HP/RPN allowed on the exam and competes directly with the TI-36x Pro and the Casio fx-991EX.
While each calculator has it's strengths, the 35s is the only programmable allowed. This, as well as RPN gives it a distinct potential advantage.
RE: (35S) 1st order Derivative (at a point) - Albert Chan - 11-11-2020 03:08 PM
Just to clarify where central difference formula come from ...
Let F(h) = f(x+h), we fit a quadratic thru 3 points: (-h, F(-h)), (0, F(0)), (h, F(h))
We assume fitted polynomial closely match F(ε), even its derivatives.
XCas> fitted(xs) := simplify(symb2poly(lagrange(xs, apply(F,xs))))
XCas> reverse(fitted([-h,0,h]))
\([F(0)\;,\;{F(h)-F(-h) \over 2h}\;,\;{F(h)+F(-h)-2 F(0) \over 2h^2}]\)
We reverse fitted polynomial coefficients, to match Taylor series of F(ε):
\(F(ε) = F(0) + ({ε\over1!})F'(0) + ({ε^2\over2!})F''(0)\;+\;... \)
Matching term by term, we have central difference formula for f'(x) and f''(x):
\(f'(x) = {f(x+h) - f(x-h) \over 2h}\)
\(f''(x) = {f(x+h) - 2f(x) + f(x-h) \over h^2}\)
---
We could automate building of derivatives formulas.
To make code efficient, we let F(k) = f(x+k*h), so fitted polynomial have no h symbol.
Because of the scaling, we need to unscale it for f derivatives.
XCas> make_diff(xs) := reverse(fitted(xs)) .* map(range(len(xs)), n->n!/h^n)
XCas> make_diff([-1,0,1])
\([F(0)\;,\;{F(1)-F(-1) \over 2h}\;,\;{F(1)+F(-1)-2 F(0) \over h^2}]\)
make_diff() returned estimate formula for [f(x), f'(x), f''(x), ...]
Note the numerator of derivatives: sum of coefficients = 0 (for flat line, all derivatives = 0)
Interestingly, cubic fit gives same formula for 2nd derivatives.
XCas> make_diff([-1,0,1,2])
\([F(0)\;,\;
{6 F(1) - 2 F(-1) - F(2) - 3 F(0) \over 6 h}\;,\;
{F(1) + F(-1) - 2 F(0) \over h^2}\;,\;
{-3 F(1) - F(-1) + F(2) + 3 F(0) \over h^3}]\)