Approximating function derivatives

[Albert Chan]

Let F(x) = f(x0 + k*h), so that F arguments are nice integers ±1, ±2
Let D = derivatives of F(x)'s, at x=0

F(1) ≈ F(0) + F'(0)*1 + F''(0)*1²/2! + F'''(0)*1³/3! --> coef = 1, 1, 1/2, 1/6
F(-1)                                                     symmetry --> coef = 1, -1, 1/2, -1/6
F(2) ≈ F(0) + F'(0)*2 + F''(0)*2²/2! + F'''(0)*2³/3! --> coef = 1, 2, 2, 4/3
F(-2)                                                     symmetry --> coef = 1, -2, 2, -4/3

Apply some symmetry before matrix inversion:

Code:

| 2 0 1  0  |           | F(1) + F(-1) |    
| 0 2 0 1/3 | × D = B = | F(1) - F(-1) |
| 2 0 4  0  |           | F(2) + F(-2) |    
| 0 4 0 8/3 |           | F(2) - F(-2) |

XCas> R := inv([[2,0,1,0],[0,2,0,1/3],[2,0,4,0],[0,4,0,8/3]])

\(\left(\begin{array}{cccc}
\frac{2}{3} & 0 & \frac{-1}{6} & 0 \\
0 & \frac{2}{3} & 0 & \frac{-1}{12} \\
\frac{-1}{3} & 0 & \frac{1}{3} & 0 \\
0 & -1 & 0 & \frac{1}{2}
\end{array}\right) \)

XCas> f, x0, h := ln, 3., 0.001/2
XCas> F(x) := f(x0 + x*h)
XCas> X := [1, -1, 2, -2]
XCas> B := map(F, X)

[1.09877894145, 1.09844560811, 1.09894556646, 1.09827889977]

XCas> B := [B[0]+B[1], B[0]-B[1], B[2]+B[3], B[2]-B[3]]

[2.19722454956, 0.00033333333642, 2.19722446623, 0.000666666691358]

XCas> D := R * B; /* = transpose(R * transpose(B) */

[1.09861228867, 0.000166666666667, -2.7777779632e-08, 9.25926002537e-12]

Scale Derivatives of F(x)|x=0 to f(x)|x=x0

XCas> D / h^range(len(D))

[1.09861228867, 0.333333333333, -0.111111118528, 0.074074080203]