Approximating function derivatives

[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)