简介：基于三次插值法求解函数极值的MATLAB代码包，面向数值分析、优化设计方向的学生与工程师，核心功能是通过三个已知数据点构造三次插值多项式，逼近目标函数曲线，再对插值多项式的导数求根，进而定位候选极值点，并通过二阶导数符号判断极大值或极小值。压缩包共包含4个m脚本文件，整体大小仅1KB，代码结构清晰，主程序负责流程控制，目标函数、导数计算与搜索区间设置分别由独立脚本完成，便于替换目标函数后直接使用。已有355人浏览学习。虽然体积小巧，但完整覆盖了从函数定义、导数计算、三次插值构造到极值点搜索与类型判别的标准流程，适合数值方法与优化设计初学者对照学习，也可作为MATLAB中函数逼近与寻优的参考模板，迁移到实际工程参数优化场景中。
1. 三次插值法在优化设计里的位置：极值搜索的“最后一步”为什么值得单独抠
在优化设计里，求极值常常卡在一维搜索上：方向已经定了，步长怎么选，才能让目标函数下降最多。三次插值法就是用两个点的函数值和导数值构造三次多项式，通过求这个多项式的极值点来逼近原问题最优点的一种快速一维搜索方法。它解决的不是“能不能找到极值”，而是“要多少次迭代才能找到极值”。对结构优化、参数标定这类需要高频线搜索的场景，三次插值法能把迭代次数压到个位数，比黄金分割和二次插值明显更省。这篇文章写给手里有梯度、目标函数光滑，想把优化设计收敛速度再提一截的工程师。以下全部基于我实际调线搜索器的经验，包括那些不得不向玄学低头的时刻。
2. 从二次插值到三次插值：极值搜索效率差在哪
一个一维搜索方法快不快，看它每轮用了多少信息，以及这些信息怎么组合。黄金分割法每轮只用一个函数值，然后按固定比例收缩区间；二次插值法用三个函数值拟合抛物线；三次插值法更进一步，用两个点的函数值和导数值，凑出四个条件去构造一条三次曲线。信息量越多，逼近真实极值点的能力越强。但信息量多不代表没有麻烦：三次多项式会出现多余的驻点、分母为零、新点越界等问题，而这些恰好是工程里最容易翻车的环节。
2.1 二次插值（抛物线法）为何不适合“最后一公里”
二次插值法也叫抛物线法，是优化设计教科书里常见的一维搜索手段。给定区间内三个点，计算函数值，用这三个点唯一确定一条抛物线，再把抛物线顶点的横坐标当作新的步长估计。它最大的优势是实现简单，连导数都不用求。但劣势也正在这里：它只用函数值，相当于把目标函数局部强行看成抛物线。一旦目标函数在这个区间不是近似二次型，抛物线顶点的位置就会偏离真实极值点，尤其是三点间距严重不对称时，新点可能偏向区间的一端，导致下一轮几乎不更新。
实际问题里，连续函数在极值点附近通常可以用二次函数近似，但“附近”的范围有多宽？二次插值法只有在区间足够窄时才表现出超线性收敛。而在优化设计刚开始时，初始区间往往很宽，函数形态可能包含拐点或噪声，这时抛物线顶点连方向都容易给错。另一个更直接的短板是：二次插值完全丢掉了导数信息。在优化设计里，只要用的是带梯度的算法，如 SQP、内点法或者伴随敏感性分析，梯度本来就是现成的。放着 g0 和 g1 不用，只拿函数值去拟合抛物线，相当于你知道当前点是上坡还是下坡，却偏要猜坡顶在哪。所以我个人的习惯是：没有导数，或者导数噪声大到不可信，才退回二次插值或黄金分割；只要能算梯度，直接上三次插值。
2.2 三次插值法的四个边界条件：从 f, f' 构造三次多项式
三次插值法的数学基础很朴素：用两个端点的函数值和导数值，构造一个三次多项式 p(x)，让 p 在这两个位置上的取值和斜率，与原目标函数完全一致。假设左侧点 x0、右侧点 x1，令 h = x1 - x0，并记 t = x - x0，则三次多项式可以写成：
p(t) = A * t^3 + B * t^2 + C * t + D
四个边界条件是：
p(0) = f0 = f(x0) p(h) = f1 = f(x1) p'(0) = g0 = g'(x0) p'(h) = g1 = g'(x1)
分别代入后，可以解出系数。D 和 C 是直接得出来的：
D = f0 C = g0
A 和 B 需要消元。由 p(h) 得到 Ah^3 + Bh^2 + Ch + D = f1，再由 p'(h) 得到 3Ah^2 + 2B*h + C = g1。把 D、C 代进去，解这个二元一次方程组：
A = [h * (g0 + g1) - 2 * (f1 - f0)] / h^3 B = [g1 - g0 - 3 * A * h^2] / (2 * h)
到这里，三次多项式的系数就完全确定了。但这只是第一步，我们要找的是这个多项式的极小点，而不是随便一个驻点。
2.3 求导数方程的根：x_new 的选择与收敛条件
对 p(t) 求导，得到 p'(t) = 3At^2 + 2Bt + C，这是个二次方程，理论上在区间内可能有两个根。其中一个对应局部极大，另一个对应局部极小；我们要选的是那个让 p''(t) = 6At + 2*B > 0 的根。手工求两个根再判断区间很啰嗦，实际工程实现里更常用一个化简后的公式，直接给出新点估计：
z = 3 * (f0 - f1) / h + g0 + g1 w = sqrt(z^2 - g0 * g1) x_new = x1 - h * (g1 + w - z) / (g1 - g0 + 2 * w)
这个公式在不少优化库的线搜索器里都能看到，本质是把两次求根、二阶导符号判断、排除多余驻点这几件事一次做完。使用它的前提是区间端点满足 g0 < 0 < g1，也就是左端函数在下降、右端函数在上升，这样区间内才有极小点。如果用户传入的区间刚好相反，最简单的处理是交换 x0 和 x1，因为很多时候只是左右端点顺序搞反了。
算出 x_new 后，要在新点处计算 g_new = g(x_new)。如果 g_new < 0，说明极小点还在 x_new 右侧，把 x_new 更新为新的左端点；如果 g_new > 0，把 x_new 更新为新的右端点；如果 g_new 接近 0，就找到了极值点。这个更新逻辑和二分法类似，但区间收缩不是等比的，每次收缩量由三次多项式决定，通常几步就能把区间压得很小。
收敛条件建议设双重判据：一是梯度绝对值小于 tol，二是区间宽度小于 xtol。梯度判据在平坦函数上很挑，因为很多函数在极值点附近梯度变化很慢，绝对值一直停在 1e-8 以上；位置判据能识别“已经动不了了”的情况，直接把当前最优位置返回。后面第3章的代码里两条判据会同时放进去。
这里还要提一个数值尺度的问题。计算 z 时包含 (f0 - f1)/h，如果 f 的量级是 1e6 而导数的量级是 1e-3，z 的舍入误差会彻底盖过梯度信息。我在实际优化过重量、温度这类不同量纲的目标时，都先做变量归一化，或者至少保证 f 和 g 在同一个数量级附近，否则 w 和分母都会失真，迭代几步就开始乱跳。
3. 用 Python 跑通三次插值法：一维搜索的完整代码与参数说明
理论知识看一百遍，不如在桌面写一个能跑的例子。这一章我会用一个单峰函数演示三次插值法的最小可用实现。目标函数选 f(x)=exp(x)-3*x，它的极小值点解析解是 x = ln(3)，能直接对照结果，不存在“不知道正确答案”的问题。
3.1 目标函数与梯度：一个能验证到 ln3 的例子
这个函数本身很简单，但很适合用来验证线搜索器。它在 x=0 处导数为 1-3=-2，下降；在 x=2 处导数为 7.389-3=4.389，上升。让初始区间取 [0, 2]，天然满足三次插值法要求的左负右正。如果换成更复杂的优化设计，目标函数可能是仿真器返回的响应面，但接口形式完全一样：只要给我一个输入 x，返回函数值和梯度。
下面的代码只依赖 Python 标准库，不用 numpy 也能跑通：
import math def f_obj(x): """目标函数：exp(x) - 3*x，极值点 x = ln(3)""" return math.exp(x) - 3.0 * x def g_obj(x): """目标函数的一阶导数""" return math.exp(x) - 3.0 def cubic_search(f, g, x0, x1, tol=1e-8, xtol=1e-10, max_iter=50): """ 三次插值法一维搜索 f: 目标函数 g: 目标函数梯度 x0, x1: 区间端点，要求极小点位于区间内 tol: 梯度收敛阈值 xtol: 区间宽度收敛阈值 max_iter: 最大迭代次数 """ f0 = f(x0) f1 = f(x1) g0 = g(x0) g1 = g(x1) # 自动整理区间方向：保证左端导数为负，右端导数为正 if g0 > 0 and g1 < 0: x0, x1 = x1, x0 f0, f1 = f1, f0 g0, g1 = g1, g0 if not (g0 < 0 and g1 > 0): raise ValueError("区间两端需满足 f'(x0)<0<f'(x1)") # 保存历史最优值点，防止极端情况下的偶发错误点污染最终结果 best_x = (x0 + x1) / 2.0 best_f = min(f0, f1) for i in range(max_iter): h = x1 - x0 if abs(h) < xtol: break z = 3.0 * (f0 - f1) / h + g0 + g1 disc = z * z - g0 * g1 w = math.sqrt(disc) if disc > 0 else 0.0 denom = g1 - g0 + 2.0 * w if abs(denom) < 1e-15: # 分母退化：三次多项式在该区间没有可靠极值点，退化为二分 x_new = (x0 + x1) / 2.0 else: x_new = x1 - h * (g1 + w - z) / denom # 强制新点落在区间内部，离端点 10% 的位置至少保留一点空间 margin = 0.1 * h x_new = min(max(x_new, x0 + margin), x1 - margin) f_new = f(x_new) g_new = g(x_new) if f_new < best_f: best_x = x_new best_f = f_new # 收缩区间：看新点的导数符号 if g_new < 0: x0, f0, g0 = x_new, f_new, g_new else: x1, f1, g1 = x_new, f_new, g_new # 双重收敛判断 if abs(g_new) < tol or abs(x1 - x0) < xtol: break if i % 10 == 0: print(f"迭代 {i:2d} x={x_new:.10f} g={g_new:.3e} h={x1-x0:.3e}") return best_x x_opt = cubic_search(f_obj, g_obj, 0.0, 2.0) print(f"最优 x = {x_opt:.12f}") print(f"ln(3) = {math.log(3.0):.12f}") print(f"最小 f = {f_obj(x_opt):.12f}")3.2 迭代流程：区间收缩里每一步在做什么
这段代码的关键逻辑在三个地方。首先是自动整理区间方向。真实需求里用户不会在意端点顺序，只要传入一个包含极小点的区间就行，所以我在函数开头加了交换逻辑，让左端导数是负、右端导数是正。这个交换看起来简单，很多粗糙实现却漏掉了，结果就是直接报错或者迭代出负数步长。
其次是分母退化处理。denom 如果接近零，标准公式会给出一个无穷远点，这时候三次多项式在整个区间上可能没有一个内点的驻点，继续算下去没有任何意义。我选择退化为中点半分，相当于动能降到最低，但至少保存了区间信息。这个分支在简单函数上几乎不会触发，但一旦遇到仿真黑盒，几乎必然会发生。很多开源的线搜索器不做这个保护，就是留着一个暗坑。
最后是 best_x 的维护。由于我们强制把新点约束在区间内部，理论上每轮都应该比上一轮好，但数值误差和强制收缩可能让新点的函数值不如历史最优。在所有优化算法里，保留历史最优都是一个不费成本的“后悔药”，最后返回它，能避免偶发的一次坏点毁掉整个线搜索。
3.3 关键参数怎么给：容差、最大迭代与初始区间
参数设置不是玄学，但要配合目标函数的尺度来调。tol 默认给 1e-8，适合梯度量级在 1e-2 到 1e2 之间的普通问题。如果梯度本身就小，比如目标函数乘以 1e-6，1e-8 就太严格了，永远达不到；这时按梯度量级的 1/1000 去给更实际。xtol 则看你对步长精度的要求，物理参数通常不需要 1e-10 的精度，给 1e-6 就好，还能省两次迭代。
max_iter=50 在我的经验里绰绰有余，因为三次插值法对接近二次形态的局部区间收敛非常快。但如果目标函数是带噪声的黑盒模拟器，50 轮可能不够，因为每一步的新点可能都不稳定，区间收缩很慢。那时我通常把 max_iter 提到 100，并配合函数值缓存来压制噪声影响，这一点在第 5 章会详细讲。
初始区间是另一个容易出错的地方。三次插值法要求区间两端满足左负右正，这是保证区间内含极小点的前提。如果你只有两个点，并不知道极小点是不是在中间，可以先跑一小段黄金分割法把区间括出来，再交给三次插值法。这个组合是我在仿真优化里最常用的姿势：前面用黄金分割收粗，后面用三次插值精收。
4. 三次插值法的适用边界：黑盒目标函数与梯度噪声下还能用吗
三次插值法对函数形态和梯度质量的要求都不低。很多工程场景里，目标函数是仿真器输出的一个数，没有解析梯度，甚至函数本身都带随机噪声。这一章专门说清楚，什么条件下三次插值法值得用，什么条件下应该果断换方法。
4.1 没有解析梯度？有限差分梯度与步长设置
没有解析梯度时，最常见的做法是中心差分。方向导数可以用下面的公式估算：
g_cal = (f(x + eps) - f(x - eps)) / (2 * eps)
在一维线搜索里，x 就是步长 α，f(α) 是沿着搜索方向的目标函数值。eps 的选择直接决定梯度质量。过小，会撞上浮点舍入误差，尤其 f 量级很大时，两个大数相减会直接损失有效数字；过大，则差分结果包含太多高阶分量，和真实梯度明显偏离。我通常先按 eps = max(1e-6, sqrt(eps_machine) * abs(x)) 起步，再根据收敛行为放大到 1e-4 或 1e-3。对噪声显著的仿真黑盒，反而要把 eps 调大，取 1e-3 甚至 1e-2，因为完全忽略噪声里的小尺度抖动，反而更容易得到一个稳定的平均梯度。
注意，中心差分的代价是每算一个梯度要两次函数值。用三次插值法时每轮需要新点处的梯度，等于每轮额外多两次函数评估。如果函数每次评估要跑一次有限元或 CFD，这个成本就很可观。这时我会考虑用少量样本先拟合一个代理模型，在代理模型上算梯度和三次插值，而不是直接对黑盒差分。
4.2 三次插值与黄金分割、二次插值的对比：选型边界
一张对比表能说清楚它们的差异：
| 方法 | 所需信息 | 单轮函数评估成本 | 收敛特征 | 最适合的场景 |
|---|---|---|---|---|
| 黄金分割法 | 函数值 | 1 次 | 线性收敛，区间稳定缩小 | 导数不可用、区间可包住极值、函数不光滑 |
| 二次插值法 | 三点函数值 | 3 次（初始），每轮 1 次 | 超线性，但区间大时易振荡 | 导数不可用、函数光滑但计算便宜 |
| 三次插值法 | 两点函数值和导数值 | 每轮 f 和 g 各一次 | 超线性，迭代少 | 梯度可得、函数光滑、计算成本高 |
| 有限差分+三次插值 | 函数值 + 近似梯度 | 每轮 1 次 f + 2 次差分评估 | 取决于差分步长 | 黑盒但函数噪声可控 |
经验判断是这样：如果函数计算一次要十秒以上，而导数能伴随算出来，三次插值法几乎是必选。SQP、内点法这类外层优化算法本身已经提供了梯度和搜索方向，线搜索器再不用导数就太亏了。如果函数不光滑，比如目标是带罚函数或绝对值项的结果，黄金分割法更稳，因为三次插值碰到导数跳变会直接误判极小值位置。
4.3 在多维优化设计里做线搜索：方向与步长的接力
多维优化里的线搜索，本质上是在当前设计点 x_k 处给一个方向 d_k，把问题化为单变量函数 φ(α) = f(x_k + α * d_k) 的极小值问题。α 是步长。φ'(α) 可由链式法则得到，如果外层算法能提供梯度，φ'(α) = ∇f(x_k + α * d_k)^T d_k。这个式子直接把多维问题拉回一维，因此前面所有公式里的 f、g 都换成 φ 和 φ'。
在实现时，我习惯把线搜索器设计成只关心一维状态，不关心你的目标函数原来是几维。这样线搜索可以复用。方向 d_k 的质量也会影响三次插值的表现：如果 d_k 是负梯度，目标函数在 α 很小时可能接近线性；如果 d_k 是拟牛顿方向，φ(α) 在 α=1 附近很可能已经接近二次形态，三次插值会很快收敛到好的步长。这也是为什么很多拟牛顿算法里，一步三次插值线搜索就够用了。
需要警惕的是，每一步线搜索都要用到 ϕ 和 ϕ' 的关系。如果外层优化算法给你的梯度有噪声，差的线搜索器会把噪声放大到步长选择里，让整个优化路线变得歪歪扭扭。这种情况下我会在线搜索之前先把梯度做一步平滑，或者把三个连续步长的函数值做二次拟合后再交给三次插值。
5. 三次插值法避坑指南：5 个现场翻车案例与排查方法
做过几年优化器的人，都知道线搜索是充满暗坑的地方。这里整理 5 个最常见的翻车案例，每一条都是我或者同事调坏过的真实问题。现象、原因、解决三步拆开写。
5.1 极值点飞出区间：根的选择逻辑没写对
现象：新点 x_new 算出来比区间端点还大，或者跑到区间另一端，导致区间不但没收缩，反而被替换成一个更宽的范围。
原因：三次多项式的导数是二次方程，两个根中只有一个是极小点。如果实现时直接用通用的二次求根公式，然后随手取了一个根，很容易取到极大点。另一个原因是端点的函数值差异太大，标准公式里的 w 失去意义，计算出的 x_new 自然越界。
解决：在代码里对新点做硬性保护，让 x_new 必须落在 [x0 + 0.1h, x1 - 0.1h] 之间。这个技巧在我给的第 3 章代码里已经埋进去了。之前我有一版代码因为没做这个保护，在目标函数是“应力约束罚函数”时，三次插值给出的步长经常跨越整个变量边界，把设计点推到铁皮模型之外，后面所有计算全部崩掉。加了边界保护后，即使多项式有问题，区间也只会收缩，不会爆炸。
5.2 初始区间两端导数不满足“左负右正”
现象：程序报错说“区间两端需满足 f'(x0)<0<f'(x1)”，或者明明程序没报错，但迭代几轮后返回的 x 根本不是极小点。
原因：调用者随便给了两个点，假设极小点在中间，但实际上区间可能完全在极小点的左侧或右侧。比如目标函数单调，两个端点导数为负，区间内没有驻点，三次插值法天然没法工作。
解决：做两件事。第一，在代码里自动交换端点处理“顺序反了”的情况，这能救回一半的场景。第二，如果交换后仍不满足左负右正，就先用一个扫描程序从粗网格里找极值的大致范围，找到后再调用三次插值。另一个办法是用黄金分割法先跑 10 轮，把区间括住，再交棒给三次插值。我现在的流程里总是先跑一轮“括区间”函数，再决定用哪种线搜索，彻底告别“假设区间有效”的脆皮代码。
5.3 有限差分梯度噪声大：差分步长不是玄学
现象：迭代开始一两轮还挺正常，后面 x_new 开始来回跳，函数值并不持续下降，甚至最后停在明显不是极值点的地方。
原因：有限差分步长选得太小，或者目标函数本身有随机噪声。差分步长太小时，两个相近函数值的小数部分完全被舍入误差污染，算出的梯度符号都可能反。噪声则是另一个维度的问题：如果仿真器每次调用因为线程调度或随机种子得到不同结果，那么差分算出的梯度完全是来回收敛噪声，三次插值法就会像喝醉一样乱走。
解决：对步长，我用公式 eps = max(1e-3, (abs(f))^(1/3)) 做保底，再根据收敛日志下调。因为三次插值本意是找一个近似平滑的局部模型，容许差分梯度有一点偏移，不像牛顿法对 Hessian 那么敏感。对噪声，先把目标函数在同一个点重复评估 3 到 5 次取平均，或者对仿真器固定随机种子，保证函数值是确定性的。如果噪声太大无法消除，就不要用三次插值，换回黄金分割法更稳，黄金分割只比大小，噪声只要不改变单调性就不会出大错。
5.4 平坦区域迭代停滞：终止条件要补一个位置判据
现象：程序跑满 max_iter 才退出，但返回的 x 和初始值相差很小。画函数图发现，整个区间像一条平线，各处梯度绝对值都低于 1e-4，但继续缩小梯度阈值仍然不动。
原因：三次插值法依赖导数符号变化来收缩区间。如果函数在极值点附近非常平坦，梯度绝对值一直很小但符号可能因为数值误差来回变，导致区间两端不断交换位置，区间宽度也缩小到一个死区，无法进一步收敛。
解决：终止条件必须同时看梯度和区间宽度。我在第 3 章代码里用的abs(g_new) < tol or abs(x1-x0) < xtol就是专门对这种情况的双保险。另外，当检测到连续 5 轮区间宽度下降不到 10% 时，直接接受当前 best_x，比死等梯度容差强。平坦问题里的梯度容差纯粹是心理安慰，位置容差才是真正让程序停止的判据。
5.5 黑盒函数不“黑”但不确定：缓存与重复计算问题
现象：同样的 x，两次调用 f(x) 返回不同的函数值，三次插值的多项式系数一会儿大一会儿小，导致最终极值点完全随机。
原因：仿真优化里的“黑盒”根本不是黑盒，而是不确定盒。有限元软件在不同核数下做并行，结果有微小波动；CFD 迭代未完全收敛就停了；或者目标函数里有随机蒙特卡洛采样。梯度差分时，不同的 x±eps 对应不同的随机种子，算出来的梯度方向纯粹是伪随机。
解决：第一，对所有评估点做永久缓存，字典的键是浮点转成十六进制后的 bit 表示，避免同一个 x 重复计算。第二，如果条件允许，固定随机种子，让函数变成严格确定性。第三，实在无法固定，就在每个评估点调用多次求平均。我在优化一个注塑工艺参数时，为了让三次插值法正常工作，把每一步的目标函数固定到平均 5 次仿真结果，最终收敛曲线才稳定。这个代价换来的是线搜索不再随机漂移，整体反而节省了时间。
6. 验证与进阶：让三次插值法在你的优化流程里真正落地
线搜索器改完以后，不要直接塞进大优化工程，先做小规模验证。我给自己定了一个固定流程：用两条曲线验证，再配合函数值缓存和预搜索，确保它在真实业务里不出幺蛾子。
6.1 三条收敛曲线：梯度、步长、函数值一起看
跑完一次线搜索，我会把每次迭代的 x、g_new、f_new 打出来，然后画三条曲线：函数值下降路径、梯度绝对值随迭代变化、区间宽度变化。三者要一起看，只看函数值会掩盖停滞，只看梯度会忽略区间收缩异常。正确的形态是函数值单调下降或近乎单调下降，梯度绝对值快速趋近零，区间宽度指数级收缩。如果函数值下降但梯度绝对值上升，说明多项式选错根了；如果区间没收缩而梯度也没到容差，说明陷入了平坦区。
6.2 两个落地技巧：函数值缓存与区间预搜索
无论是不是黑盒，我都给线搜索器配一个全局缓存。因为外层优化算法会在多个地方调用同一个点，比如差分梯度时 f(x+eps) 和 f(x-eps) 很有可能是之前某个 α 计算过了。缓存命中一次，相当于省一次昂贵仿真。
第二个技巧是区间预搜索。三次插值法需要一个可靠的初始区间，最简单的做法是先按等比数列尝试步长 α=0.1, 0.3, 1.0, 3.0，直到发现目标函数从下降到上升，取相邻两个点作为区间。这个过程会多花几次函数评估，但能避免大量无效的迭代，总体仍然划算。
6.3 我踩过的最后一个坑：保护点在区间内
很多年前我第一次手搓三次插值法，觉得标准公式高不可攀，照着论文推导实现后，直接用生产数据跑，结果步长飞出边界，把整个优化方向带偏。后来在一行代码里加了边界限制 x_new = min(max(x_new, x0+0.1*（x1-x0)), x1-0.1*（x1-x0))，问题立刻消失。从那时候起，我所有优化器里的线搜索都保留“强制保护半径”，并且保留历史最优点作为返回结果。这个习惯救过我太多次，现在连给别人写插件都会不自觉地带上它。如果你只能从这篇文章带走一句话，那就是：三次插值法本身不玄，真正玄的是你给它的区间和导数，先把这两个稳住，再谈收敛速度。希望帮到你。
本文还有配套的精品资源，点击获取
