Approximating function derivatives

[Namir]

Albert,

I translated your last pseudo-code into Excel VBA. The code works well and gives accurate results With x0 = 6, the code takes 8 iterations (like my Newto'ns version).

The real savings in function callis adapting your fx averaging approach to Halley's method. This way the fx averaging approach will use two function calls per iteration as opposed to 3 function calls when I calculate f(x), f(x+h) and f(x-h).

Namir

[Pekis]

(01-27-2024 08:04 AM)Pekis Wrote:  I've found a more iterative way, and it seems interesting:

f'(x0) ~= a * f(x0-2h) + b * f(x0-h) + c * f(x0+h) + d * f(x0+2h):
=> a = 1/(12h), b = -8/(12h), c = 8/(12h), d = -1/(12h)

But after all, if the first derivative f'(x0) ~= a * f(x0-2h) + b * f(x0-h) + c * f(x0+h) + d * f(x0+2h) (already very good precision), then the 2nd derivative f"(x0) ~= a * f'(x0-2h) + b * f'(x0-h) + c * f'(x0+h) + d * f'(x0+2h) will have good precision, and so on, with the same a, b, c, d coefficients from above:

So, if you want it, f"(x0) ~= a * f'(x0-2h) + b * f'(x0-h) + c * f'(x0+h) + d * f'(x0+2h):
the f'(x0-2h) term has to be calculated as a * f(x0-4h) + b * f(x0-3h) + c * f(x0-h) + d * f(x0) (since f'(x0) ~= a * f(x0-2h) + b * f(x0-h) + c * f(x0+h) + d * f(x0+2h))
the f'(x0-h) term has to be calculated as a * f(x0-3h) + b * f(x0-2h) + c * f(x0) + d * f(x0+h)
the f'(x0+h) term has to be calculated as a * f(x0-h) + b * f(x0) + c * f(x0+2h) + d * f(x0+3h)
the f'(x0+2h) term has to be calculated as a * f(x0) + b * f(x0+h) + c * f(x0+3h) + d * f(x0+4h)
We can reuse previous calculations but also have to add new ones:
f(x0-4h), f(x0-3h), f(x0+3h), f(x0+4h) [and f(x0)].

(01-27-2024 03:34 PM)Albert Chan Wrote:  Hi, Pekis

Yes, your way work!

I think you just re-discovered divided finite-difference (see post #7)
...

Well, that's more a Franken Central Finite Difference calculation :
f"(x0) is actually evaluated as 1/(144h^2) * (f(x0-4h) - 16 * f(x0-3h) + 64 * f(x0-2h) + 16 * f(x0-h) - 130 * f(x0) + 16 * f(x0+h) + 64 * f(x0+2h) - 16 * f(x0+3h) + f(x0+4h))

while the real Central Finite difference for f"(x0) up to 4h (nine points, including x0) (see Derivative 2, Accuracy 8) is
1/(5040h^2) * (-9 * f(x0-4h) + 128 * f(x0-3h) -1008 * f(x0-2h) + 8064 * f(x0-h) - 14350 * f(x0) + 8064 * f(x0+h) - 1008 * f(x0+2h) + 128 * f(x0+3h) - 9 * f(x0+4h))

In order to find the real Central Finite Difference coefficients up to 4h (nine points, including x0), you have to solve this system, based on Taylor expansion similar to what I presented in the first post (except that I added a coefficient for f(x0))

For f"(x0): (See Derivative 2 , Accuracy 8)
a + b + c + d + e + f + g + h + i = 0
-4a - 3b -2c - d + 0 + f + 2g + 3h + 4i = 0
16a + 9b + 4c + d + 0 + f + 4g + 9h + 16i = 2!(h^2)
-64a - 27b - 8c - d + 0 + f + 8g + 27h + 64i = 0
256a + 81b + 16c + d + 0 + f + 16g + 81h + 256i = 0
-1024a - 243b -32c - d + 0 + f + 32g + +243h + 1024i = 0
4096a + 729b + 64c + d + 0 + f + 64g + 729h + 4096i = 0
-16384a - 2187b - 128c - d + 0 + f + 128g + 2187h + 16384i = 0
65536a + 6561b + 256c +d + 0 + f + 256g + 6561h + 65536i = 0

giving 1/(5040h^2) * (-9 * f(x0-4h) + 128 * f(x0-3h) -1008 * f(x0-2h) + 8064 * f(x0-h) - 14350 * f(x0) + 8064 * f(x0+h) - 1008 * f(x0+2h) + 128 * f(x0+3h) - 9 * f(x0+4h))

For f'(x0): (See Derivative 1 , Accuracy 8)
a + b + c + d + e + f + g + h + i = 0
-4a - 3b -2c - d + 0 + f + 2g + 3h + 4i = 1/h
16a + 9b + 4c + d + 0 + f + 4g + 9h + 16i = 0
-64a - 27b - 8c - d + 0 + f + 8g + 27h + 64i = 0
256a + 81b + 16c + d + 0 + f + 16g + 81h + 256i = 0
-1024a - 243b -32c - d + 0 + f + 32g + +243h + 1024i = 0
4096a + 729b + 64c + d + 0 + f + 64g + 729h + 4096i = 0
-16384a - 2187b - 128c - d + 0 + f + 128g + 2187h + 16384i = 0
65536a + 6561b + 256c +d + 0 + f + 256g + 6561h + 65536i = 0

giving 1/(840h) * (3 * f(x0-4h) - 32 * f(x0-3h) + 168 * f(x0-2h) - 672 * f(x0-h) + 0 * f(x0) + 672 * f(x0+h) - 168 * f(x0+2h) + 32 * f(x0+3h) - 3 * f(x0+4h))

[Albert Chan]

(01-30-2024 10:32 AM)Pekis Wrote:  Well, that's more a Franken Central Finite Difference calculation :

I like your setup better than the "proper" central difference fit.
Spreadsheet is trivial to implement, copy *same* derivative formula to cells!
Neighoring nearness of slopes is reassuring, gives us a heads up for typos.
It is also easy to understand, slope of slope of slope ...

Doing the calculation by fitting all points may not get much better.
Bigger N-th difference is just noise ... garbage in, garbage out.

(01-29-2024 02:59 AM)Paul Dale Wrote:  Everyone needs to remember that numerical derivative algorithms are generally very unstable.

As a reminder, function is ln, x0 = 3, h = 0.001/2, and the d(k)'s k-th differences.

Code:

 P     X          Y           d1          d2          d3            d4          d5          d6
                                                             
 -3  2.9985 1.0981121636                                     
                          1.6674E-04                         
 -2   2.999 1.0982788998             -2.7796E-08             
                          1.6671E-04              9.2661E-12 
 -1  2.9995 1.0984456081             -2.7787E-08             -4.6629E-15
                          1.6668E-04              9.2615E-12              2.2204E-16
  0       3 1.0986122887  1.6667E-04 -2.7778E-08  9.2593E-12 -4.4409E-15  1.1102E-16 -2.2204E-16
                          1.6665E-04              9.2570E-12              0.0000E+00
  1  3.0005 1.0987789414             -2.7769E-08             -4.4409E-15
                          1.6663E-04              9.2526E-12
  2   3.001 1.0989455665             -2.7759E-08
                          1.6660E-04
  3  3.0015 1.0991121637

I added mean of odd central difference to row center (P=0).
This allowed direct computation of derivatives estimate.
Source: Fundamentals of Numerical Analysis, 1975, Kellison, p. 152

(hD) = δ - 1/24*δ^3 + 3/640*δ^5 - ...

What is nice is we can treat operator as expression, and get other derivatives!

(hD)^2 = (δ - 1/24*δ^3 + 3/640*δ^5 - ...)^2 = δ^2 - 1/12*δ^4 + 1/90*δ^6 - ...
(hD)^3 = (δ - 1/24*δ^3 + 3/640*δ^5 - ...)^3 = δ^3 - 1/8 * δ^5 + ...
(hD)^4 = (δ - 1/24*δ^3 + 3/640*δ^5 - ...)^4 = δ^4 - 1/6 * δ^6 + ...

(hD)^k f(0) is simply dot product of center row and corresponding coefficients.

Code:

  D                       1.0000E+00             -4.1667E-02              4.6875E-03
  1   DOT = 0.3333333356  3.3333E-01  0.0000E+00 -7.7161E-10  0.0000E+00  1.0408E-15  0.0000E+00  
                                                             
                                      1.0000E+00             -8.3333E-02              1.1111E-02
  2         -0.111111111  0.0000E+00 -1.1111E-01  0.0000E+00  1.4803E-09  0.0000E+00 -9.8686E-12
                                                             
                                                  1.0000E+00             -1.2500E-01
  3         0.0740739692  0.0000E+00  0.0000E+00  7.4074E-02  0.0000E+00 -1.1102E-07  0.0000E+00
                                                             
                                                              1.0000E+00             -1.6667E-01
  4         -0.070462154  0.0000E+00  0.0000E+00  0.0000E+00 -7.1054E-02  0.0000E+00  5.9212E-04

From size of numbers, (hD)^k ≈ δ^k, the others contributed very little.

BTW, this is done in Lotus 123! I still retained muscle memory

Correction:

Fundaments of Numerical Analysis, p154 Wrote:However, there is a problem in using these central difference formulas to find the odd derivatives
at the tabulated points of collocations, since the required differences do not appear in the difference table ...

Author provided alternative formula, with mean operator µ (exactly what I did to fill missing holes)

(hD) = µ*(δ - 1/6*δ^3 + 1/30*δ^5 - ...)

Other odd derivatives can be built the same way

(hD)^3
= (hD) * (hD)^2
= µ*(δ - 1/6*δ^3 + 1/30*δ^5 - ...) * (δ^2 - 1/12*δ^4 + 1/90*δ^6 - ...)
= µ*(δ^3 - 1/4*δ^5 + ...)

With these adjusted coefficients updated to spreadsheet, I get

D¹ f(0): 0.3333333356 --> 0.3333333333, much better
D³ f(0): 0.0740739692 --> 0.0740738582, slightly worse

3rd derivative getting worse because µδ^5 is garbage.
With only µδ^3 term, D³ f(0) = 0.0740740802

This showed more points may make it worse. Less is more.
(or, spread the points apart, with bigger h)

[Albert Chan]

We can get all derivative formulas in XCas. Here is same code for HP Prime Cas
Note that code does not care about symmetry. Use whatever points you had available.

Cas> fitted(xs) := simplify(symb2poly(lagrange(xs, apply(f,xs))))
Cas> make_diff(xs) := reverse(fitted(xs)) .* map(range(len(xs)), n->n!/h^n)

Cas> make_diff([-1,0,1])
Cas> simplify(Ans)

\(\displaystyle [ f(0), \frac{f(1) - f(-1)}{2h}, \frac{f(1)-2f(0)+f(-1)}{h^2}] \)

With many points, resulting derivatives may be messy.

For central differences, we may use symmetry to make formula more compact.
We know sum of f's gives even derivatives, differences gives odd.

Let g^n = (f(n) ± f(-n)), except for g^0 = 1 = f(0)
Sign depends on parity of derivatives. (odd = -, even = +)

Cas> make_diff([-2,-1,1,2]) | f = (x-> when(x<0,0,g^x))
Cas> simplify(Ans)

\( \displaystyle [
\frac{-g^2+4g}{6},
\frac{-g^2+8g}{12h},
\frac{g^2-g}{3 h^2},
\frac{g^2-2g}{2 h^3}
] \)

Undo g's, we get back f's

(01-26-2024 08:33 AM)Thomas Klemm Wrote:  \(
\begin{bmatrix}
f(x_0) \\
{f}''(x_0) h^2 \\
\end{bmatrix}
\approx
\begin{bmatrix}
\frac{2}{3} & - \frac{1}{6} \\
- \frac{1}{3} & \frac{1}{3} \\
\end{bmatrix}
\cdot
\begin{bmatrix}
f_{+1} + f_{-1} \\
f_{+2} + f_{-2} \\
\end{bmatrix}
\)

\(
\begin{bmatrix}
{f}'(x_0) h \\
{f}'''(x_0) h^3 \\
\end{bmatrix}
\approx
\begin{bmatrix}
\frac{2}{3} & - \frac{1}{12} \\
- 1 & \frac{1}{2} \\
\end{bmatrix}
\cdot
\begin{bmatrix}
f_{+1} - f_{-1} \\
f_{+2} - f_{-2} \\
\end{bmatrix}
\)