如何求解一维三次插值函数的最大值？能否通过求导找极值？

当然可以通过求导的方式来找这个三次插值函数的最大值！你当前用插值数组找最大值的方法虽然简单，但求导找极值的方式会更精准——尤其是当极值点不在你采样的x_new点上时。

原理说明 #

scipy.interpolate.interp1d的kind='cubic'生成的是分段三次Hermite插值多项式，它在每个子区间上都是三次多项式，因此它的导数是分段二次多项式，完全可以通过找导数为零的点来定位极值点。

具体实现步骤 #

我们可以按以下步骤来做：

1. 定义插值函数的导数函数；
2. 在插值的整个区间内，找到导数为零的所有点；
3. 计算这些点以及区间端点的函数值，找到其中的最大值。

完整代码示例 #

import numpy as np
from scipy.interpolate import interp1d
from scipy.optimize import root_scalar
import matplotlib.pyplot as plt

# 用户原始数据
x_axis = np.array([2.14414414, 2.15270826, 2.16127238, 2.1698365, 2.17840062, 
                   2.18696474, 2.19552886, 2.20409298, 2.2126571, 2.22122122])
y_axis = np.array([0.67958442, 0.89628424, 0.78904004, 3.93404167, 6.46422317, 
                   6.40459954, 3.80216674, 0.69641825, 0.89675386, 0.64274198])

# 创建三次插值函数
f = interp1d(x_axis, y_axis, kind='cubic')

# 定义导数函数：用中心差分求数值导数，精度足够
def df(x):
    # 取x附近三个点计算梯度，得到中心差分的导数
    x_samples = np.linspace(x-1e-6, x+1e-6, 3)
    return np.gradient(f(x_samples))[1]

# 遍历每个子区间，寻找导数零点
extrema_points = []
for i in range(len(x_axis)-1):
    x_left = x_axis[i]
    x_right = x_axis[i+1]
    try:
        # 用brentq方法在区间内找根，要求区间两端导数值符号相反
        sol = root_scalar(df, bracket=[x_left, x_right], method='brentq')
        if sol.converged:
            extrema_points.append(sol.root)
    except ValueError:
        # 区间内没有导数零点，直接跳过
        continue

# 把区间端点加入候选点，避免漏掉端点极值
candidate_points = np.concatenate([[x_axis[0]], extrema_points, [x_axis[-1]]])
candidate_values = f(candidate_points)

# 找到最大值对应的点和值
max_idx = np.argmax(candidate_values)
max_x = candidate_points[max_idx]
max_y = candidate_values[max_idx]

print(f"最大值点：x={max_x:.6f}, y={max_y:.6f}")

# 绘图展示结果
x_new = np.linspace(x_axis[0], x_axis[-1], 200)
fig, ax = plt.subplots()
ax.plot(x_new, f(x_new), label='三次插值曲线')
ax.scatter(max_x, max_y, color='red', s=100, marker='*', label='最大值点')
ax.scatter(x_axis, y_axis, color='blue', label='原始数据点')
ax.legend()
plt.show()

更简洁的替代方案 #

如果你不想手动处理求导和找根的细节，也可以直接用scipy.optimize.minimize_scalar来最大化插值函数——本质是通过最小化-f(x)来实现，内部会自动处理求导和极值查找：

from scipy.optimize import minimize_scalar

# 最小化-f(x)等价于最大化f(x)
res = minimize_scalar(lambda x: -f(x), bounds=(x_axis[0], x_axis[-1]), method='bounded')
max_x = res.x
max_y = f(max_x)
print(f"最大值点：x={max_x:.6f}, y={max_y:.6f}")

这个方法代码更简洁，适合快速实现需求。

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