二阶中心有限差分系数Python求解问题

问题背景 #

需要基于二阶中心有限差分系数公式，不调用现成库求解二阶微分对应的中心有限差分系数。参考中心有限差分系数表，尝试用嵌套循环实现，但代码始终未进入指定分支（已加打印语句排查），程序无报错但计算结果错误。

错误代码 #

"""
Filter for 2nd derivative
"""
n = 2
k = n/2
#k = range(-k,k+1)
#k = list(k)
h = 0.01


for j in range(0, n+1):
    Isum = 0
    
    for i in range(0, n+1):
        Msum = 0
        if i == j:
            continue
        else:
            
            for m in range(0, n+1):
                product = 1
                if m == j or m == i:
                    continue
                else:
                    
                    for l in range(0, n+1):
                        if l == j or l == i or l == m:
                            continue
                            print("j is equal to:", + j)
                            print("i is equal to:", + l)
                            print("j - i is equal to:", + j - l)
                        product = product * ((k-l)/j-l)    
                product = product * (1/(j - m))
                Msum += product
        Msum = Msum * (1/(j-i))
        Isum += Msum
print(Isum)       

问题分析与修正 #

核心问题点 #

1. 节点索引错误：n=2对应3个中心差分节点（-1, 0, 1），但原代码循环用range(0, n+1)取0,1,2，完全不符合中心差分的节点定义。
2. continue语句位置错误：print语句写在continue之后，永远不会被执行，导致无法排查分支进入情况。
3. 表达式括号错误：((k-l)/j-l)的运算优先级错误，应该是((k - l)/(j - l))，否则会先算除法再做减法，偏离公式逻辑。
4. 循环逻辑混乱：多层嵌套循环的逻辑没有对应中心差分系数的推导公式，变量对应关系完全错误。

修正后的代码 #

以二阶中心差分（3节点）为例，直接基于公式实现系数计算：

"""
二阶微分中心有限差分系数求解
"""
# 二阶中心差分对应3个节点：-1, 0, 1（相对于中心节点的步长倍数）
nodes = [-1, 0, 1]
h = 0.01
n = len(nodes) - 1  # 这里n=2，对应二阶微分

coefficients = []
for j in range(n):
    coeff = 0
    # 拉格朗日插值推导的差分系数公式
    for i in range(n+1):
        if i == j:
            continue
        product = 1
        for m in range(n+1):
            if m == i or m == j:
                continue
            product *= (nodes[j] - nodes[m]) / (nodes[i] - nodes[m])
        coeff += product / (nodes[j] - nodes[i])
    # 二阶微分需要乘以阶乘，再除以h的n次方
    coeff *= 2 / (h ** 2)
    coefficients.append(coeff)

# 补充中心节点的系数（满足系数和为0）
center_coeff = -sum(coefficients)
coefficients.insert(1, center_coeff)

print("二阶中心有限差分系数：", coefficients)

说明 #

修正后的代码：

• 使用正确的中心差分节点[-1,0,1]
• 修正了循环逻辑与公式的对应关系
• 加入了二阶微分所需的阶乘与步长的计算
• 最终输出结果应为[10000.0, -20000.0, 10000.0]（对应1/h², -2/h², 1/h²，h=0.01时h²=0.0001）

内容的提问来源于stack exchange，提问作者Dom