Approximating function derivatives
|
01-29-2024, 08:43 PM
(This post was last modified: 01-29-2024 08:46 PM by Namir.)
Post: #21
|
|||
|
|||
RE: Approximating function derivatives
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 |
|||
01-30-2024, 10:32 AM
Post: #22
|
|||
|
|||
RE: Approximating function derivatives
(01-27-2024 08:04 AM)Pekis Wrote: I've found a more iterative way, and it seems interesting: (01-27-2024 03:34 PM)Albert Chan Wrote: Hi, Pekis 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)) |
|||
01-30-2024, 10:08 PM
(This post was last modified: 02-03-2024 09:04 PM by Albert Chan.)
Post: #23
|
|||
|
|||
RE: Approximating function derivatives
(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 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 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 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) |
|||
01-30-2024, 11:06 PM
Post: #24
|
|||
|
|||
RE: Approximating function derivatives
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: \( |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)