(35S) 1st order Derivative (at a point)

11112020, 11:11 AM
(This post was last modified: 11132020 09:52 AM by trojdor.)
Post: #1




(35S) 1st order Derivative (at a point)
First order Numerical (explicit) Differentiation program for HP 35s
The companion program, 2nd order Derivative, can be found at https://www.hpmuseum.org/forum/thread15875.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, 1E6 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 TI36x estimate f' at a point.) f'(x) ~ (f(x+h)f(xh)) / 2h (D.Levy EQ 5.4 Introduction to Numerical Analysis) I also chose to use an h value of 1E6, 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 1E6 D004 + D005 XEQ D014 D006 RCL D D007 1E6 D008  D009 XEQ D014 D010  D011 2E6 < 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^321*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 TI36x Pro and the Casio fx991EX. 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. ENTER > = 

11112020, 03:08 PM
Post: #2




RE: (35S) 1st order Derivative (at a point)
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(xh) \over 2h}\) \(f''(x) = {f(x+h)  2f(x) + f(xh) \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}]\) 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)