Why does my Python code produce incorrect finite difference weights from a Taylor series expansion?

Below is the key insight: one of the rows in your system is set up with an incorrect sign. In your derivation, the Taylor expansion for f(x – 2Δx) contributes a term to f‴(x) with a negative sign, but in your matrix you used +8 for that coefficient.

Derivation Quick Review

Expanding about x:

• f(x + 2h) = f(x) + 2h·f′(x) + (4h²/2)·f″(x) + (8h³/6)·f‴(x) + (16h⁴/24)·f⁽⁴⁾(x) + …
• f(x + h) = f(x) + h·f′(x) + (h²/2)·f″(x) + (h³/6)·f‴(x) + (h⁴/24)·f⁽⁴⁾(x) + …
• f(x – h) = f(x) – h·f′(x) + (h²/2)·f″(x) – (h³/6)·f‴(x) + (h⁴/24)·f⁽⁴⁾(x) + …
• f(x – 2h) = f(x) – 2h·f′(x) + (4h²/2)·f″(x) – (8h³/6)·f‴(x) + (16h⁴/24)·f⁽⁴⁾(x) + …

After clearing the denominators by multiplying the appropriate rows, you want your system to satisfy:

1. f(x) coefficient: a + b + c + d + e = 0
2. f′(x) coefficient: 2a + b – d – 2e = 1
3. f″(x) coefficient (after multiplying by 2): 4a + b + d + 4e = 0
4. f‴(x) coefficient (after multiplying by 6): 8a + b – d – 8e = 0
5. f⁽⁴⁾(x) coefficient (after multiplying by 24): 16a + b + d + 16e = 0

What’s Wrong in Your Code

In your Python matrix, the fourth row is given as

[8, 1, 0, -1, 8] 

But based on the derivation it should be

[8, 1, 0, -1, -8]

That sign on the coefficient multiplying e must be negative, otherwise the system is not representing the Taylor expansion correctly and you end up with the wrong weights.

Corrected Code

Below is the corrected Python snippet:

PYTHONCopied
import numpy as np

A = np.array([
    [1,  1, 1,  1,  1],     # f(x) coefficients
    [2,  1, 0, -1, -2],     # f'(x) coefficients
    [4,  1, 0,  1,  4],     # f''(x) coefficients (multiplied by 2)
    [8,  1, 0, -1, -8],     # f'''(x) coefficients (multiplied by 6)
    [16, 1, 0,  1, 16]      # f⁽⁴⁾(x) coefficients (multiplied by 24)
])

b = np.array([0, 1, 0, 0, 0])  # Targeting the first derivative

z = np.linalg.solve(A, b)
print(z)

When you run this corrected code, you should obtain the weights:

  a = -1/12, b = 2/3, c = 0, d = -2/3, e = 1/12

These match the theoretical five-point central difference scheme for the first derivative:

  f'(x) ≈ (-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h))⁄(12h)

Summary

Your original matrix had the wrong sign in the row corresponding to the third derivative (f‴(x)). Correcting the sign to –8 in that row will yield the expected finite-difference weights.