(35S) 2nd order Derivative (at a point)

11092020, 09:51 PM
(This post was last modified: 11132020 10:02 AM by trojdor.)
Post: #1




(35S) 2nd order Derivative (at a point)
2nd order Numerical (explicit) Differentiation program for HP 35s
This program is by Michael Lyons, and was inspired by the work by Alexis Zuleta. OVERVIEW This program is an addition to the first order explicit derivative program found at https://www.hpmuseum.org/software/35derivp.htm By assigning this second order program to the G key, just below the D key used for the 1st order pgm, (and sharing the same F / function,) it's easy to remember and use. 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) (It's surprisingly efficient. Note that while D.Levy states that a centered differencing formula is more accurate, in the real world testing I performed with a CD program I wrote, I could see no real benefit to replacing the Zuleta pgm.) My d" program is based on the following 2nd order approximation of the second derivative: f"(x) ~ (f(x+h)2f(x)+f(xh)) / h^2 (D.Levy EQ 5.6 Introduction to Numerical Analysis) Note that the h value needs to be larger for the 2nd derivative approximation (I use 1E3), otherwise the calculator's precision limitation will cause the program to fail. 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! G001 LBL G G002 INPUT G G003 STO X G004 XEQ F001 G005 2 G006 * G007 +/ G008 x<>y G009 1E3 G010 + G011 XEQ G019 G012 RCL G G013 1E3 G014  G015 XEQ G019 G016 1E6 < note that is h^2, not h (that's why it's different from step 9 and 13) G017 / G018 RTN G019 STO X G020 XEQ F001 G021 x<>y G022 CLx G023 R(down) G024 + G025 RTN EXAMPLE (Updated for better example per AC suggestion...) 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 G [ENTER] 7. Displays G? 8. Enter the point to evaluate, in this case '4', then [R/S]. 9. See RUNNING, then the result. In this case, 270. Let's manually solve to confirm: f'(x)=x^4 + 5x^3 21x^2 =4x^3 + 15x^2 42x f"(x)= =12x^2 + 30x 42  x=4 =192 + 120  42 =270 < 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 HP is the only programmable allowed. This, as well as RPN gives it a distinct potential advantage. ENTER > = 

11102020, 02:20 PM
Post: #2




RE: (35S) 2nd order Derivative (at a point)
(11092020 09:51 PM)trojdor Wrote: EXAMPLE This is a bad example to test effect of h on 2nd derivative estimate. XCas> f(x) := horner([a4,a3,a2,a1],x) XCas> simplify( (f(x+h)  2*f(x) + f(xh)) / (h*h) ) → 2*a3 + 6*a4*x For cubics, central difference 2nd derivatives is independent of h. XCas> f(x):= x^3 + 5*x^2  21*x XCas> (f(x+h)  2*f(x) + f(xh)) / (h*h)  x=4, h=1000 → 34 BTW, for decimal calculator, we might get slightly better accuracy by splitting f(x): (this avoids possible rounding errors of 2*f(x)) c = f(x) f''(x) ≈ ((f(x+h)  c)  (c  f(xh))) / (h*h) Quote:Zuleta uses this simple (forward distancing) approximation of the first derivative: For fair comparison, let's try both ways for f'(4), with same h = 1E3: XCas> (f(x+h)  f(x)) / h  x=4, h=1E3 → 67.017001 XCas> (f(x+h)  f(xh)) / (2*h)  x=4, h=1E3 → 67.000001 We can confirm all derivatives by expanding f(x+4): XCas> expand(f(x+4)) → x^3 + 17*x^2 + 67*x + 60 f(4) = 60 f'(4) = 67*1! = 67 f''(4) = 17*2! = 34 f'''(4) = 1*3! = 6 

11102020, 11:23 PM
(This post was last modified: 11112020 04:39 AM by trojdor.)
Post: #3




RE: (35S) 2nd order Derivative (at a point)
(11102020 02:20 PM)Albert Chan Wrote:(11092020 09:51 PM)trojdor Wrote: EXAMPLE Albert, thank you for your response, and for pointing that out. What would you suggest as a good example? (11102020 02:20 PM)Albert Chan Wrote: For fair comparison, let's try both ways for f'(4), with same h = 1E3: Well, to be really fair, let's try on the 35s with h=1E6 as used in his pgm: (f(x+h)  f(x)) / h  x=4, h=1E6 → 67.0000000000 (f(x+h)  f(xh)) / (2*h)  x=4, h=1E6 → 67.0000000000 Which is why I stated that I could see no real/practical benefit to replacing the smaller Zuleta pgm with a larger program that ends up with the same practical accuracy (to 10 decimal places) on the HP 35s. This program was written for the HP 35s to assist in the NCEES / Fundamentals of Engineering Licensing Exam. (The only HP/RPN calculator allowed on the exam is the HP 35s.) Thus my focus tends to be more engineering/results oriented, as opposed to the theoretical. However, I truly do appreciate your response/advice/feedback....thanks so much! (And I'm always up for learning something new.) ...now I'm off to experiment with your idea of splitting f(x) to minimize rounding errors... ENTER > = 

11112020, 12:35 AM
Post: #4




RE: (35S) 2nd order Derivative (at a point)
(11102020 11:23 PM)trojdor Wrote: What would you suggest as a good example? We wanted a test function with errors dependent on size of h. For polynomials, use quartic or higher degree. XCas> f(x) := x*(x^3 + 5*x^2  21*x) XCas> expand(f(x+4)) → x^4+21*x^3+135*x^2+328*x+240 f(4) = 240 f'(4) = 328*1! = 328 f''(4) = 135*2! = 270 // Estimate f'(4), with different h's XCas> (f(x+h)  f(x)) / h  x=4, h=1e3 → 328.135021001 XCas> (f(x+h)  f(x)) / h  x=4, h=1e6 → 328.000134999 XCas> (f(x+h)  f(xh)) / (2h)  x=4, h=1e3 → 328.000021 XCas> (f(x+h)  f(xh)) / (2h)  x=4, h=1e6 → 328.00000001 // Estimate f''(4), with different h's XCas> t1 := (f(x+h)  2*f(x) + f(xh))/(h*h)  x=4, h=1e1 → 270.02 XCas> t2 := (f(x+h)  2*f(x) + f(xh))/(h*h)  x=4, h=1e2 → 270.0002 XCas> t3 := (f(x+h)  2*f(x) + f(xh))/(h*h)  x=4, h=1e3 → 270.000002132 h=1e3 is slightly too small. Without rounding errors, t3 = 270.000002 (exactly) Richardson Extrapolation from t1, t2 gives excellent estimated for f''(4): XCas> t2 + (t2t1) / (1001) → 270.0 

11112020, 04:11 AM
(This post was last modified: 11112020 04:32 AM by trojdor.)
Post: #5




RE: (35S) 2nd order Derivative (at a point)
Thank you, kind sir.
(Interestingly, the rounding errors of the calculator seem to work in favor of the algorithms chosen. When I use my 35s programs to solve the 4th deg problem for both f' and f", I get 328.000000000 and 270.000000000 respectively.) You may also have convinced me that my other, center differencing program for the 1st derivative est. may have some practical value after all. Rather than toss it out, I'll type it up and post it later. It's pretty compact. Thanks again for all the time/effort, mike ENTER > = 

« Next Oldest  Next Newest »

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