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

[Albert Chan]

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}]\)