Approximating function derivatives

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