Python3非等间隔采样数据的二阶导数估算、定积分计算与函数验证

核心计算逻辑（适配非等间隔采样） #

• 导数计算：采用非均匀节点二阶精度三点拉格朗日差商法，不预设固定步长。对每个采样点取相邻三个有效节点构造局部二次插值多项式，求导得到该点的g'(x)估算值，全节点的计算精度满足二阶要求。首点取前三个节点做前向三点插值，末点取最后三个节点做后向三点插值，无硬编码初始导数值的逻辑。
• 定积分计算：采用非均匀分段二次积分法，不要求采样点等间隔。基于相邻三点构造的二次插值多项式，逐段累加区间内的积分贡献，达到二阶积分精度，比传统梯形法的一阶精度更稳定。

可直接运行的Python3实现代码 #

import numpy as np
import matplotlib.pyplot as plt

def load_dataset(file_path="data.txt"):
    """读取双列data.txt数据集，自动按x值升序排序，适配非等间隔采样"""
    data = np.loadtxt(file_path, skiprows=1)
    x_raw = data[:, 0]
    g_raw = data[:, 1]
    sort_index = np.argsort(x_raw)
    return x_raw[sort_index], g_raw[sort_index]

def calc_second_order_derivative(x, g):
    """计算非均匀采样下的g'(x)二阶精度估算值"""
    node_num = len(x)
    derivative = np.zeros(node_num)
    # 首节点：前向三点插值求导
    x0, x1, x2 = x[0], x[1], x[2]
    g0, g1, g2 = g[0], g[1], g[2]
    derivative[0] = (
        g0 * (2*x0 - x1 - x2) / ((x0 - x1)*(x0 - x2))
        + g1 * (x0 - x2) / ((x1 - x0)*(x1 - x2))
        + g2 * (x0 - x1) / ((x2 - x0)*(x2 - x1))
    )
    # 内部节点：相邻三点插值求导
    for i in range(1, node_num-1):
        x_im1, x_i, x_ip1 = x[i-1], x[i], x[i+1]
        g_im1, g_i, g_ip1 = g[i-1], g[i], g[i+1]
        derivative[i] = (
            g_im1 * (x_i - x_ip1) / ((x_im1 - x_i)*(x_im1 - x_ip1))
            + g_i * (2*x_i - x_im1 - x_ip1) / ((x_i - x_im1)*(x_i - x_ip1))
            + g_ip1 * (x_i - x_im1) / ((x_ip1 - x_i)*(x_ip1 - x_im1))
        )
    # 末节点：后向三点插值求导
    x_nm2, x_nm1, x_n = x[-3], x[-2], x[-1]
    g_nm2, g_nm1, g_n = g[-3], g[-2], g[-1]
    derivative[-1] = (
        g_nm2 * (x_n - x_nm1) / ((x_nm2 - x_n)*(x_nm2 - x_nm1))
        + g_nm1 * (x_n - x_nm2) / ((x_nm1 - x_n)*(x_nm1 - x_nm2))
        + g_n * (2*x_n - x_nm2 - x_nm1) / ((x_n - x_nm2)*(x_n - x_nm1))
    )
    return derivative

def calc_second_order_integral(x, g):
    """计算非均匀采样下[a,b]区间g(x)的二阶精度定积分值"""
    node_num = len(x)
    integral = 0.0
    # 逐段基于三点二次插值计算积分
    for i in range(0, node_num - 2):
        x0, x1, x2 = x[i], x[i+1], x[i+2]
        g0, g1, g2 = g[i], g[i+1], g[i+2]
        h0 = x1 - x0
        h1 = x2 - x1
        segment_int = (h0 + h1)/6 * (
            g0*(2*h0 - h1)/h0
            + g1*(h0 + h1)**2/(h0*h1)
            + g2*(2*h1 - h0)/h1
        )
        # 每两个相邻区间合并累加，避免重复计算
        if i % 2 == 0:
            integral += segment_int
    # 奇数节点时，最后一个单区间用梯形公式补全
    if node_num % 2 == 0:
        integral += (x[-1] - x[-2]) * (g[-1] + g[-2]) / 2
    return integral

if __name__ == "__main__":
    # 加载数据
    x, g = load_dataset()
    interval_left, interval_right = x.min(), x.max()
    # 数值计算
    g_prime = calc_second_order_derivative(x, g)
    integral_value = calc_second_order_integral(x, g)
    # 绘制双y轴关系曲线
    fig, ax_left = plt.subplots(figsize=(10, 6))
    ax_left.plot(x, g, color="#1f77b4", linewidth=2, label="g(x)")
    ax_left.set_xlabel("x")
    ax_left.set_ylabel("g(x)", color="#1f77b4")
    ax_left.tick_params(axis="y", labelcolor="#1f77b4")
    ax_right = ax_left.twinx()
    ax_right.plot(x, g_prime, color="#ff7f0e", linestyle="--", linewidth=2, label="g'(x)")
    ax_right.set_ylabel("g'(x)", color="#ff7f0e")
    ax_right.tick_params(axis="y", labelcolor="#ff7f0e")
    plt.title(f"g(x)与g'(x)关系曲线（积分区间[{interval_left:.2f}, {interval_right:.2f}]）")
    plt.show()
    # 打印积分结果
    print(f"区间[{interval_left:.4f}, {interval_right:.4f}]上g(x)的定积分计算结果：{integral_value:.6f}")
    # 曲线形态判断原函数类型
    derivative_std = np.std(g_prime)
    derivative_mean = np.mean(np.abs(g_prime))
    if derivative_std / derivative_mean < 0.05:
        print("原函数类型判断：g(x)为线性函数，g'(x)近似为常数")
    elif np.std(np.diff(g_prime)) / np.mean(np.abs(np.diff(g_prime))) < 0.05:
        print("原函数类型判断：g(x)为二次多项式函数，g'(x)为线性函数")
    elif np.corrcoef(g, g_prime)[0, 1] > 0.98:
        print("原函数类型判断：g(x)为指数函数，g'(x)与g(x)线性正相关")
    else:
        print("原函数类型判断：g(x)为三角函数类，g'(x)为同频偏移的周期函数")

已知问题的修正说明 #

• 移除硬编码固定步长h=0.05的逻辑，所有步长计算完全基于实际读取的x值间隔，支持任意随机非等间隔采样场景
• 移除硬编码首个导数估算值的逻辑，首末点自动通过相邻三点插值计算导数，无手动赋值内容
• 修正三角函数拟合的逻辑缺陷：不强制对全采样点做全局三角函数拟合，而是基于数值计算得到的曲线特征匹配函数类型，避免全局拟合的过拟合、龙格现象问题
• 修正牛顿差商插值的逻辑缺陷：不使用全节点高阶差商，仅使用局部三点的低阶插值，数值计算稳定，精度可达二阶要求

数值结果验证方法 #

根据判断出的原函数类型，代入对应精确导数、精确积分解析式做定性对比：

• 线性函数g(x)=mx+c：精确导数为常数m，精确积分为0.5*m*(b²-a²)+c*(b-a)，对比数值计算结果的残差，残差小于5%即验证有效
• 二次多项式g(x)=px²+qx+r：精确导数为2px+q，精确积分为p*(b³-a³)/3 + q*(b²-a²)/2 + r*(b-a)，对比残差验证
• 正弦函数g(x)=A*sin(ωx+φ)：精确导数为Aω*cos(ωx+φ)，精确积分为A/ω*(-cos(ωb+φ)+cos(ωa+φ))，对比残差验证

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