如何在Python中对非均匀网格使用有限差分法求导？

针对非均匀的能量(e)-压强(p)网格，要计算导数$f = dp/de$，不能全用简单的一阶向前差分，得按点的位置选择合适的差分格式来保证精度，尤其是前几个点的精细计算需求，具体方案如下：

分点处理策略 #

1. 第一个点（i=0）：二阶向前差分 #

当网格点数≥3时，用二阶公式比一阶精度高，公式为：

$f_0 = \frac{-(3p_0 - 4p_1 + p_2)}{2(e_1 - e_0)}$
如果只有2个点，退化为一阶向前差分。

2. 中间点（0 < i < n-1）：二阶中心差分（适配非均匀网格） #

利用相邻三个点的间距加权计算，避免均匀网格的限制，公式为：

$f_i = \frac{p_{i+1}(e_i - e_{i-1}) - p_i(e_{i+1} - e_{i-1}) + p_{i-1}(e_{i+1} - e_i)}{(e_{i+1} - e_i)(e_i - e_{i-1})(e_{i+1} - e_{i-1})/2}$
本质是用三个点拟合二次曲线后求导，比一阶差分精度更高。

3. 最后一个点（i=n-1）：二阶向后差分 #

和首点对应，利用倒数三个点计算，公式为：

$f_{n-1} = \frac{3p_{n-1} - 4p_{n-2} + p_{n-3}}{2(e_{n-1} - e_{n-2})}$

Python代码实现 #

import numpy as np

def calculate_dp_de(e_vals, p_vals):
    n_points = len(e_vals)
    dp_de = np.zeros(n_points)
    
    # 处理第一个点
    if n_points >= 3:
        dp_de[0] = -(3 * p_vals[0] - 4 * p_vals[1] + p_vals[2]) / (2 * (e_vals[1] - e_vals[0]))
    elif n_points == 2:
        dp_de[0] = (p_vals[1] - p_vals[0]) / (e_vals[1] - e_vals[0])
    
    # 处理中间点
    for i in range(1, n_points - 1):
        h_prev = e_vals[i] - e_vals[i-1]
        h_next = e_vals[i+1] - e_vals[i]
        numerator = p_vals[i+1] * h_prev - p_vals[i] * (h_prev + h_next) + p_vals[i-1] * h_next
        denominator = (h_prev * h_next * (h_prev + h_next)) / 2
        dp_de[i] = numerator / denominator
    
    # 处理最后一个点
    if n_points >= 3:
        dp_de[-1] = (3 * p_vals[-1] - 4 * p_vals[-2] + p_vals[-3]) / (2 * (e_vals[-1] - e_vals[-2]))
    elif n_points == 2:
        dp_de[-1] = (p_vals[-1] - p_vals[-2]) / (e_vals[-1] - e_vals[-2])
    
    return dp_de

# 示例测试
e = np.array([0.5, 1.2, 2.0, 3.5, 4.8])
p = np.array([1.1, 2.3, 3.0, 4.7, 6.2])
result = calculate_dp_de(e, p)
print("dp/de计算结果：", result)

可选方案：三次样条插值求导 #

如果对精度要求更高，或者网格点分布极不均匀，可以用三次样条先拟合p(e)曲线，再求导，代码更简洁且稳定：

from scipy.interpolate import CubicSpline

# 拟合样条曲线
spline = CubicSpline(e, p)
# 计算各点的一阶导数
dp_de_spline = spline(e, 1)
print("样条插值求导结果：", dp_de_spline)

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