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

[trojdor]

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, 1E-6 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(x-h)) / 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 1E-3),
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 1E-3
G010 +
G011 XEQ G019
G012 RCL G
G013 1E-3
G014 -
G015 XEQ G019
G016 1E-6 <-- 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^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 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 TI-36x Pro and the Casio fx-991EX.
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.

[Albert Chan]

(11-09-2020 09:51 PM)trojdor Wrote:  EXAMPLE
If f(x)= x^3 + 5*x^2 - 21*x
Calcuate f"(x) at x=4

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(x-h)) / (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(x-h)) / (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(x-h))) / (h*h)

Quote: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)

(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.)

For fair comparison, let's try both ways for f'(4), with same h = 1E-3:

XCas> (f(x+h) - f(x)) / h | x=4, h=1E-3                → 67.017001
XCas> (f(x+h) - f(x-h)) / (2*h) | x=4, h=1E-3       → 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

[trojdor]

(11-10-2020 02:20 PM)Albert Chan Wrote:  

(11-09-2020 09:51 PM)trojdor Wrote:  EXAMPLE
If f(x)= x^3 + 5*x^2 - 21*x
Calcuate f"(x) at x=4

This is a bad example to test effect of h on 2nd derivative estimate.

Albert, thank you for your response, and for pointing that out.
What would you suggest as a good example?

(11-10-2020 02:20 PM)Albert Chan Wrote:  For fair comparison, let's try both ways for f'(4), with same h = 1E-3:

XCas> (f(x+h) - f(x)) / h | x=4, h=1E-3                → 67.017001
XCas> (f(x+h) - f(x-h)) / (2*h) | x=4, h=1E-3       → 67.000001

Well, to be really fair, let's try on the 35s with h=1E-6 as used in his pgm:
(f(x+h) - f(x)) / h | x=4, h=1E-6                         → 67.0000000000
(f(x+h) - f(x-h)) / (2*h) | x=4, h=1E-6                → 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...

[Albert Chan]

(11-10-2020 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=1e-3             → 328.135021001
XCas> (f(x+h) - f(x)) / h | x=4, h=1e-6             → 328.000134999

XCas> (f(x+h) - f(x-h)) / (2h) | x=4, h=1e-3      → 328.000021
XCas> (f(x+h) - f(x-h)) / (2h) | x=4, h=1e-6      → 328.00000001

// Estimate f''(4), with different h's

XCas> t1 := (f(x+h) - 2*f(x) + f(x-h))/(h*h) | x=4, h=1e-1       → 270.02
XCas> t2 := (f(x+h) - 2*f(x) + f(x-h))/(h*h) | x=4, h=1e-2       → 270.0002
XCas> t3 := (f(x+h) - 2*f(x) + f(x-h))/(h*h) | x=4, h=1e-3       → 270.000002132

h=1e-3 is slightly too small. Without rounding errors, t3 = 270.000002 (exactly)

Richardson Extrapolation from t1, t2 gives excellent estimated for f''(4):

XCas> t2 + (t2-t1) / (100-1)       → 270.0

[trojdor]

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