I'm trying to calculate the weights for a finite-difference approximation of the first derivative f′(x)f'(x)f′(x) using the Taylor series expansion. I'm solving for weights a,b,c,d,ea, b, c, d, ea,b,c,d,e such that:
af(x+2Δx)+bf(x+Δx)+cf(x)+df(x−Δx)+ef(x−2Δx)a f(x+2\Delta x) + b f(x+\Delta x) + c f(x) + d f(x-\Delta x) + e f(x-2\Delta x)af(x+2Δx)+bf(x+Δx)+cf(x)+df(x−Δx)+ef(x−2Δx)
approximates f′(x)f'(x)f′(x). Here's what I did:
I used the Taylor series expansion for f(x±kΔx)f(x \pm k\Delta x)f(x±kΔx), where k=1,2k = 1, 2k=1,2.
I built a system of linear equations to enforce the following conditions:
Coefficients of f(x): a+b+c+d+e=0
Coefficients of f′(x): 2a+b−d−2e=1
Coefficients of f′′(x): 4a+b+d+4e=0
Coefficients of f(3)(x): 8a+b−d−8e=0
Coefficients of f(4)(x): 16a+b+d+16e=0
I implemented the matrix equation A⋅z=bA \cdot z = bA⋅z=b in Python:
import numpy as np A = np.array([ [1, 1, 1, 1, 1], # Coefficients of f(x) [2, 1, 0, -1, -2], # Coefficients of f'(x) [4, 1, 0, 1, 4], # Coefficients of f''(x) [8, 1, 0, -1, 8], # Coefficients of f'''(x) [16, 1, 0, 1, 16] # Coefficients of f''''(x) ]) b = np.array([0, 1, 0, 0, 0]) # Targeting the first derivative z = np.linalg.solve(A, b) print(z)
The Issue:
The output I'm getting is:
[0.25,0,-0,0,-0.25]
However, the expected weights for the first derivative should be something like:
-
[-1/12,2/3,0,-2/3,1/12]
What I Tried:
Double-checked the coefficients in matrix A to ensure they match the Taylor series expansion.
Verified that the right-hand side vector b=[0,1,0,0,0] is correct for approximating f′(x).
Despite this, the weights are incorrect. Am I missing something in the matrix setup or the Python implementation?
Expected Behavior:
I want the solution to match the theoretical weights for a central finite difference approximation of the first derivative f′(x) using five points.
I'm trying to calculate the weights for a finite-difference approximation of the first derivative f′(x)f'(x)f′(x) using the Taylor series expansion. I'm solving for weights a,b,c,d,ea, b, c, d, ea,b,c,d,e such that:
af(x+2Δx)+bf(x+Δx)+cf(x)+df(x−Δx)+ef(x−2Δx)a f(x+2\Delta x) + b f(x+\Delta x) + c f(x) + d f(x-\Delta x) + e f(x-2\Delta x)af(x+2Δx)+bf(x+Δx)+cf(x)+df(x−Δx)+ef(x−2Δx)
approximates f′(x)f'(x)f′(x). Here's what I did:
I used the Taylor series expansion for f(x±kΔx)f(x \pm k\Delta x)f(x±kΔx), where k=1,2k = 1, 2k=1,2.
I built a system of linear equations to enforce the following conditions:
Coefficients of f(x): a+b+c+d+e=0
Coefficients of f′(x): 2a+b−d−2e=1
Coefficients of f′′(x): 4a+b+d+4e=0
Coefficients of f(3)(x): 8a+b−d−8e=0
Coefficients of f(4)(x): 16a+b+d+16e=0
I implemented the matrix equation A⋅z=bA \cdot z = bA⋅z=b in Python:
import numpy as np A = np.array([ [1, 1, 1, 1, 1], # Coefficients of f(x) [2, 1, 0, -1, -2], # Coefficients of f'(x) [4, 1, 0, 1, 4], # Coefficients of f''(x) [8, 1, 0, -1, 8], # Coefficients of f'''(x) [16, 1, 0, 1, 16] # Coefficients of f''''(x) ]) b = np.array([0, 1, 0, 0, 0]) # Targeting the first derivative z = np.linalg.solve(A, b) print(z)
The Issue:
The output I'm getting is:
[0.25,0,-0,0,-0.25]
However, the expected weights for the first derivative should be something like:
-
[-1/12,2/3,0,-2/3,1/12]
What I Tried:
Double-checked the coefficients in matrix A to ensure they match the Taylor series expansion.
Verified that the right-hand side vector b=[0,1,0,0,0] is correct for approximating f′(x).
Despite this, the weights are incorrect. Am I missing something in the matrix setup or the Python implementation?
Expected Behavior:
I want the solution to match the theoretical weights for a central finite difference approximation of the first derivative f′(x) using five points.
2 Answers
Reset to default 2You have a typo in the 4th line of the matrix A, the last element should be -8 instead of 8:
A = np.array([
[1, 1, 1, 1, 1], # Coefficients of f(x)
[2, 1, 0, -1, -2], # Coefficients of f'(x)
[4, 1, 0, 1, 4], # Coefficients of f''(x)
[8, 1, 0, -1, -8], # Coefficients of f'''(x)
[16, 1, 0, 1, 16] # Coefficients of f''''(x)
])
Write the Taylor series out carefully:
Hence,
where
Inverting:
The weights on the points (x-2Δx,x-Δx,x,x+Δx,x+2Δx) (note that these are the opposite way round to yours) are the relevant rows in the weight matrix
For example (as in the code below) you get (Δx)f' from row number 1 (counting from 0) of the matrix W:
[ 8.33333333e-02 -6.66666667e-01 2.96059473e-16 6.66666667e-01
-8.33333333e-02]
which is, essentially, [1/12,-2/3,0,2/3,-1/12] (reminder again: my points are the opposite way round to yours, so reversing the list).
import numpy as np
A = np.array( [ # coefficients in f = A (deltax fderiv)
[ 1, -2, 2, -4/3, 2/3 ],
[ 1, -1, 1/2, -1/6, 1/24 ],
[ 1, 0, 0, 0, 0 ],
[ 1, 1, 1/2, 1/6, 1/24 ],
[ 1, 2, 2, 4/3, 2/3 ]
] )
W = np.linalg.inv( A )
derivative = 1 # this is (deltax)f'
print( W[derivative,:] )
Az=b
given thez
that was found in this script. However, ifc=[-1/12,2/3,0,-2/3,1/12]
, thenAc
is not equal tob
. Have you written your weights correctly in your script? – Slavensky Commented Jan 18 at 13:29