python - Why does my finite difference weights calculation for Taylor series give incorrect results? - Stack Overflow

admin管理员组

文章数量:1320627

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:

1. 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.

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

3. 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)
   
   

4. 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:

5. [-1/12,2/3,0,-2/3,1/12]
   

6. What I Tried:

   1. Double-checked the coefficients in matrix A to ensure they match the Taylor series expansion.

   2. 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?

7. 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:

1. 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.

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

3. 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)
   
   

4. 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:

5. [-1/12,2/3,0,-2/3,1/12]
   

6. What I Tried:

   1. Double-checked the coefficients in matrix A to ensure they match the Taylor series expansion.

   2. 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?

7. 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.

• python
• numpy
• linear-algebra
• derivative
• taylor-series

Share Improve this question Follow asked Jan 18 at 13:20 Bhavninder Singh VirdiBhavninder Singh Virdi 2333 bronze badges 2

• 1 It is a bit difficult to read the math in this question. However, it can easily be verified that Az=b given the z that was found in this script. However, if c=[-1/12,2/3,0,-2/3,1/12], then Ac is not equal to b. Have you written your weights correctly in your script? – Slavensky Commented Jan 18 at 13:29
• Please correct your formulas at the beginning of your question (they appear 3 times each). Also the \Delta, \cdot, ... needs updating. – rehaqds Commented Jan 18 at 16:21

Add a comment  | 

2 Answers 2

Sorted by: Reset to default 2

You 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,:] )

本文标签：