How do i fix this code to find the central finite difference coefficients for the second order?

Question

I am relatively new to coding and am trying to code a subroutine that takes in a grid spacing h and a polynomial degree n and returns the n + 1 numbers that correspond to the symmetric finite difference stencil/filter for the operator d2/dx2, i.e write some code that returns the finite difference coefficients found on wikipedia. I am close but not getting the correct values for 4th, 6th and 8th degree accuracy. Any help would be greatly appreciated!

I've coded this Lagrange polynomial, but it doesnt seem to work (Edit: The code doesn't return any errors, but does not return the correct coefficients for the second order central finite difference found here https://en.wikipedia.org/wiki/Finite_difference_coefficient. Thanks for the help!) Edit: The formula i am trying to implement is which i understand is a Lagrange polynomial that works as a 'filter', however i do not know why it produces the finite central difference coefficients. :

def central_difference_coefficients(order, n):
    coefficients = []

    for j in range(n+1):
        coefficient = 0.0
        for i in range(n+1):
            if i != j:
                term1 = 1.0 / (j - i)
                product = 1.0
                for m in range(n+1):
                    if m != j and m != i:
                        denominator = j - m
                        if denominator != 0:
                            product *= 1.0 / denominator
                            for l in range(n+1):
                                if l != j and l != i and l != m:
                                    denominator = j - l
                                    if denominator != 0:
                                        product *= (n / 2 - l) /denominator
                term2 = term1 * product
                coefficient += term2
        if j != 0:
            coefficient *= 1.0 / (order**2)
        coefficients.append(coefficient)
    return coefficients
# Example usage:
order = 2
n = 4
coefficients_2nd_order = central_difference_coefficients(order, n)
print(f"Central Difference Coefficients (2nd order) for n={n}:\n{coefficients_2nd_order}")

For an example of the answer it gives for n=2 is [1.0, -0.5, 0.0625], which should be [1,-2,1] Let me know where I've gone wrong!

Answer 1

I spotted at least one bug in your code: the resetting line product = 1.0 should be inside the loop indexed by m, because there is a different product for every m. Instead you have put this line just outside the loop over m, so that subsequent terms of the m sum will be affected by previous products.

In fact, you could remove line product = 1.0 completely, and replace product *= 1.0 / denominator by product = 1.0 / denominator, which would also reset the product correctly.

To avoid these kinds of bugs, python conveniently offers a sum and a prod function which you can use to calculate a sum or a product with a for-loop without manually writing product = 1.0 and product *= ....

Here are two ways to rewrite your function using sum and prod:

from math import prod

def C(h, n):
    return [
        sum(
            sum(
                prod(
                    (n / 2 - l) / (j - l)
                    for l in range(n+1)
                    if l not in (i,j,m)
                ) / (j - m)
                for m in range(n+1)
                if m not in (i, j)
            ) / (j - i)
            for i in range(n+1)
            if i != j
        ) / (h**2)
        for j in range(n+1)
    ]

def K(h, n):
    def factor_l(n, j, l):
        return (n/2 - l) / (j - l)
    def term_m(n, j, i, m):
        return (1 / (j - m)) * prod(factor_l(n, j, l)
                                    for l in range(n+1)
                                    if l not in (j,i,m))
    def term_i(n, j, i):
        return (1 / (j - i)) * sum(term_m(n, j, i, m)
                                   for m in range(n+1)
                                   if m not in (j, i))
    def coeff_j(n, j):
        return sum(term_i(n, j, i) for i in range(n+1) if i != j)
    return [coeff_j(n, j) / (h**2) for j in range(n+1)]

sum and prod are particularly convenient when implementing formula; you can compare directly the code with the formula and don't have to worry about maintaining variables and resetting them correctly.

These kinds of "inline for-loops" differ from usual for-loops and are called "generator-expressions" in python, and you can use them with sum, prod, max, min, or even with [ ] to create a list like I did with the loop indexed by j (an expression to create a list like this is called a "list comprehension").

You can test that C(1, 2) == K(1, 2) == [1.0, -2.0, 1.0] as expected.