前阵子帮一位做运动控制的同行看代码他用一段自研的“一元三次方程求解器”去算轨迹曲线与某个平面的交点结果在某个区间算出来的值总是偏那么一点点肉眼看不出来但下游的插补算法会因此抖一下。他把问题定位到求根模块上翻来覆去看那段从某处抄来的“牛莱公式”越看越糊涂——公式里既没“牛”也没“莱”只有一堆立方根和虚数单位。这就是“一元三次方程的解-牛莱公式”这个说法在工程圈里流传的典型状态名字叫得很响内容却经常被拼错、被简化、被误用。这篇文章想把这件事一次说透一元三次方程的闭式解到底长什么样牛顿迭代法在求根里扮演什么角色微积分那套分析手段为什么绕不开以及在真实代码里怎么把这几个东西组合成一个能上生产的求解器。不管你是刚学到三次方程的学生还是需要在项目里手写一个数值求根模块的工程师下面这些内容都能直接拿去用。1. 拆标题这几个名字凑在一起到底在说什么1.1 三次方程的闭式解并不叫“牛莱公式”先把最容易出错的地方点破。一元三次方程的求根公式标准叫法是卡尔达诺公式Cardanos formula也有人按推导路径叫它“根式解”。它给的是 ax³bx²cxd0 在复数域上的三个根形式上由平方根和立方根嵌套构成。而“牛顿-莱布尼茨公式”是微积分基本定理写出来是 ∫ₐᵇ f(x)dx F(b) − F(a)讲的是定积分与不定积分的关系它跟解方程求根压根不是一个问题域。那为什么网上会冒出“牛莱公式解三次方程”这种说法我观察下来有三种来源。第一种是口口相传时的缩写把“牛顿迭代法”简称为“牛顿法”再被人补上半个“莱布尼茨”就成了四不像的组合词。第二种是把微积分里分析函数形态的那套工具求导、判断单调性、找极值、用介值定理定位零点统称为“牛顿-莱布尼茨那一套”确实没有这些工具你连根有几个都说不清。第三种纯粹是关键词拼凑搜索词库里“三次方程”“牛顿迭代”“公式法”经常同框出现久而久之被误当成一个专有名词。所以我在实际工作里给出的定义是这样的一元三次方程的解 闭式解卡尔达诺公式负责精确形式与复数根 数值解牛顿迭代负责工程精度与通用性 微积分分析导数与极值负责根的定位与个数判定。这三块拼起来才是那个被叫歪了的“牛莱公式”真正想表达的东西。理解这一点很关键因为你在代码里如果只用闭式解会遇到浮点抵消只用迭代法会遇到初值发散两样都不做导数分析你连迭代区间都划不出来。1.2 牛顿迭代法在求根流程里的实际位置牛顿迭代法的核心就一行xₙ₊₁ xₙ − f(xₙ)/f′(xₙ)。几何意义是把当前点处的切线延长到与 x 轴相交用交点作为下一次猜测。它的好处在收敛速度——只要初值落在根附近且导数不为零误差是平方级缩减的也就是说有效数字大概每迭代一次翻一倍。但它在三次方程上有两个必须知道的脾气。第一它只保证局部收敛。三次函数 f(x)x³pxq 在 p0 时是“先增后减再增”的波浪形从不同初值出发会掉进不同的根从极值点附近出发甚至会来回震荡不收敛。第二遇到重根会退化成线性收敛。单根处 f′(r)≠0收敛快二重根处 f′(r)0牛顿法的收敛因子退化为 1/2 左右迭代次数直接翻好几倍而且精度上不去。这就解释了为什么工程代码里很少让牛顿法“裸奔”。常见的做法是先用微积分手段把实轴切成若干段单调区间确认每段最多一个根再在每个区间内用牛顿法或二分法收口。这个“先定位后精化”的流程才是把牛顿法用稳的关键也是后面第 3 章要给完整代码的原因。1.3 导数分析为什么是绕不开的第一步很多人拿到三次方程就直奔公式其实在动手之前用导数花三十秒做一次形态分析能省掉后面半小时的调试。对首一化之后的 f(x) x³ px q二次项已经消掉这一步第 2 章会讲怎么消导数 f′(x) 3x² p 的符号决定了整个函数的长相当p ≥ 0时f′ 恒非负函数在整个实轴上严格单调递增那么有且只有一个实根。这种情况最省心随便取一个足够远的初值做牛顿迭代都能收敛。当p 0时f′ 有两个零点 x ±√(−p/3)函数在这两点分别取得极大值和极小值于是根的个数就取决于这两个极值是否跨越 x 轴。跨了就三个实根正好相切就两个含一个二重根都在同侧就只剩一个实根。这套分析看起来是教科书内容但它在代码里的价值极高它直接告诉你“该求几个根”“每个根大致在哪个区间”从而把牛顿法的初值选择从“猜”变成“算”。我见过太多实现是把初值硬编码成 0 或 1然后在某些系数下直接不收敛报出来的错还特别难查。2. 数学底座把三次方程拆到能下手为止2.1 首一化与消二次项韦达代换的推导标准形式 ax³ bx² cx d 0第一步除以 aa≠0得到 x³ bx² cx d 0这里 b、c、d 已经是被除过的系数了。此时二次项还在公式会变得非常臃肿所以要把它消掉。做法是平移代换 x t − b/3。为什么是 b/3因为对 (t − h)³ b(t − h)² 展开t² 项的系数是 −3h b让它等于零就得到 h b/3。这个思路和二次方程配方是一回事只不过二次配方是为了凑完全平方三次平移是为了消掉次高项本质都是在“对齐”多项式的对称中心。代入展开后得到 t³ pt q 0其中p c − b²/3q 2b³/27 − bc/3 d这两条式子必须背准符号错一个后面全废。我自己的记忆方式是p 是“c 减去 b 平方的三分之一”q 是“2b³/27 减 bc/3 加 d”。如果你不放心代入 t 0即 x −b/3验证一下f(−b/3) 应该等于 q这个恒等式可以用来单测你的系数换算函数非常好用。得到 t³ pt q 0 之后三个根 t₁、t₂、t₃ 与原方程的三个根只差一个平移量 −b/3最后统一加回去就行。这一步的意义在于之后的判别式、求根公式、三角解法全部只依赖两个参数 p 和 q代码实现会清爽非常多。2.2 判别式一眼看出根的类型对 t³ pt q 0定义判别量 Δ (q/2)² (p/3)³。注意这个 Δ 的符号约定在不同教材里不一样有的用 −(4p³27q²)两者相差一个正的常数因子和一个负号判断结论是对应的别混用就行。我这里统一用 Δ (q/2)² (p/3)³Δ 的符号根的情况处理策略Δ 0一个实根 一对共轭复根直接用立方根公式数值稳定Δ 0至少两个根重合三重根或一单一双单独分支避免除以零Δ 0三个互不相同的实根转三角解法避开复数开立方这个判别式不是天上掉下来的它可以由 1.3 节的极值分析直接推出来这个推导我在下面 2.4 节会完整走一遍看完你就再也不会记错符号了。实际写代码时有个细节必须提醒Δ 0 的判断不能用 0。因为 p、q 都是浮点数理论上为零的量在机器里可能是 1e-18 也可能是 −3e-17。正确做法是构造一个与量纲匹配的容差具体怎么定在 4.1 节讲。2.3 卡尔达诺公式的完整推演有了 t³ pt q 0核心技巧是设 t u v。代入(uv)³ p(uv) q u³ v³ 3uv(uv) p(uv) q 0整理成 u³ v³ q (3uv p)(uv) 0。这时候我们手里有两个自由变量 u 和 v可以额外加一个约束条件让括号里那项消失即令3uv p 0也就是 uv −p/3。剩下的部分就要求 u³ v³ −q。于是 u³ 和 v³ 成了一个新的二次方程的两个根设 U u³、V v³则 U V −qUV (uv)³ −p³/27。由韦达定理U、V 是 z² qz − p³/27 0 的两根。解这个二次方程z [−q ± √(q² 4p³/27)] / 2 −q/2 ± √((q/2)² (p/3)³) −q/2 ± √Δ到这里一切顺畅取 U −q/2 √Δ、V −q/2 − √Δ然后 u ∛U、v ∛Vt u v 就是一个实根。剩下两个根呢利用立方根在复数域上有三个分支这个事实把 u 乘上三次单位根 ω (−1 i√3)/2同时 v 乘上 ω² (−1 − i√3)/2保证 uv 乘积不变也就是约束条件 3uvp0 始终成立就得到全部三组解t₁ u vt₂ ωu ω²vt₃ ω²u ω⁴vω⁴ ω这个推导的价值在于它把“三次”问题降成了“二次”问题——这就是为什么三次方程能用根式解而五次及以上不能伽罗瓦理论的地盘这里不展开。对写代码的人来说记住U、V 是同一个二次方程的两根这一点就够了因为所有数值稳定性问题都出在这对根的相减上。2.4 不可约情形三角解法的来龙去脉上一节最后一步有个陷阱当 Δ 0 时√Δ 是纯虚数U 和 V 变成一对共轭复数然后你要对复数开立方再相加。数学上没问题结果一定是实数但在浮点运算里两个模长很大、虚部相反的复数相加实部会发生灾难性抵消有效数字成片丢失。这就是著名的“不可约情形”历史上逼着数学家正视复数现在逼着程序员改算法。出路是换一个表达形式。既然 Δ0 意味着 p0那么 −p/3 0可以设 t 2√(−p/3)·cosθ。令 m 2√(−p/3)代入 t³ ptt³ pt m³cos³θ pm·cosθ由 m² −4p/3 得 p −3m²/4于是上式 m³cos³θ − (3m³/4)cosθ (m³/4)(4cos³θ − 3cosθ) (m³/4)cos3θ用的是三倍角恒等式 cos3θ 4cos³θ − 3cosθ。原方程 t³ptq0 于是变成 (m³/4)cos3θ q 0即cos3θ −4q/m³m³ 8(−p/3)^{3/2}代进去化简得到 cos3θ −q / (2(−p/3)^{3/2})等价写成 cos3θ (3q)/(2p)·√(−3/p)两种写法我都见过代码里我用前者的变形 3q/(p·m)。剩下就是反解θ 有三组值 θ (arccos(...) 2kπ)/3k0,1,2代回 t m·cosθ 就得到三个实根。全程只有实数运算没有虚数没有抵消精度稳稳的。顺手把 2.2 节的判别式推一遍三个实根的条件是极大值大于零、极小值小于零。由 f′(x)0 得 x±m/2不对重算——f(x)x³pxqf′(x)3x²p0 的解是 x ±√(−p/3)注意 m 2√(−p/3)所以临界点是 ±m/2。代入f(m/2) (m/2)³ p(m/2) q用 p −3m²/4 得 (m³/8) − (3m³/8) q q − m³/4 f(−m/2) −(m³/8) (3m³/8) q q m³/4三个实根要求 q − m³/4 0 q m³/4即 |q| m³/4 2(−p/3)^{3/2}。两边平方q² 4(−p/3)³ −4p³/27即 q²/4 p³/27 0也就是 Δ (q/2)² (p/3)³ 0。推导闭环判别式的符号再也不会记反了。3. 代码实操三条实现路径与取舍3.1 路径一按判别式分流的闭式解这是最直接的一条路把第 2 章的三套分支翻译成代码。要点有三个判别式的容差要与量纲匹配、立方根要自己写以处理负数、输出统一成 (实部, 虚部) 的形式方便调用方使用。import math def _cbrt(x): # math.pow 对负数开立方会返回 nan这里手动处理符号 return math.copysign(abs(x) ** (1.0 / 3.0), x) def solve_cubic(a, b, c, d, rel_tol1e-13): 返回三个根 [(re, im), ...]a 为三次项系数 if a 0.0: raise ValueError(a 不能为 0请退化为二次方程处理) b, c, d b / a, c / a, d / a # 首一化 p c - b * b / 3.0 # 消二次项 q 2.0 * b ** 3 / 27.0 - b * c / 3.0 d shift -b / 3.0 # 最后平移回去的量 hq q / 2.0 tp p / 3.0 delta hq * hq tp ** 3 # 容差按量纲缩放delta 的量纲是 p^3 或 q^2 scale max(abs(hq * hq), abs(tp ** 3), 1e-300) eps rel_tol * scale if delta eps: # 一实两复 s math.sqrt(delta) u _cbrt(-hq s) v _cbrt(-hq - s) t1 u v re -(u v) / 2.0 im math.sqrt(3.0) / 2.0 * (u - v) return [(t1 shift, 0.0), (re shift, im), (re shift, -im)] if delta -eps: # 三个不同实根三角解法 m 2.0 * math.sqrt(-tp) arg 3.0 * q / (p * m) # 等价于 (3q)/(2p)*sqrt(-3/p) arg max(-1.0, min(1.0, arg)) # 夹紧防止 acos 域外报错 theta math.acos(arg) / 3.0 out [] for k in range(3): t m * math.cos(theta - 2.0 * math.pi * k / 3.0) out.append((t shift, 0.0)) return sorted(out) # 重根分支delta 落在容差内 u _cbrt(-hq) return [(2.0 * u shift, 0.0), (-u shift, 0.0), (-u shift, 0.0)]这段代码有几个刻意设计的地方值得说明。arg那一行做了max(-1.0, min(1.0, arg))的夹紧因为当判别式非常接近零时浮点误差可能让 arg 略微超出 [−1,1]math.acos会直接抛ValueError这在线上是致命的。_cbrt单独抽出来是因为 Python 的x ** (1/3)对负底数返回复数或 nan取决于版本和写法自己在外面处理符号最稳。三角解法那三个根算完之后我用sorted排了序这是为了满足调用方“根按从小到大排列”的隐含约定这个约定在几何求交、区间判定场景里非常重要。3.2 路径二牛顿迭代配合单调区间分段闭式解虽然精确但系数病态时依然会掉精度而且不容易处理“只想要实根”“只要某个区间内的根”这类需求。这时候用迭代法更合适。完整流程是先求导找临界点把实轴切段逐段用介值定理判断有没有根再用牛顿法精化。def newton(f, df, x0, tol1e-14, maxit200): x x0 for i in range(maxit): fx, dfx f(x), df(x) if dfx 0.0: x 1e-8 * max(1.0, abs(x)) # 导数踩零小扰动后重试 continue step fx / dfx x - step # 相对判据兼顾大根和小根 if abs(step) tol * max(1.0, abs(x)): return x, i 1, True return x, maxit, False def real_roots(a, b, c, d, tol1e-14): if a 0.0: raise ValueError(a 不能为 0) b, c, d b / a, c / a, d / a f lambda x: ((x b) * x c) * x d # 秦九韶少两次乘法 df lambda x: (3.0 * x 2.0 * b) * x c # 根的界|x| 1 max(|b|,|c|,|d|)取对称区间 bound 1.0 max(abs(b), abs(c), abs(d)) inner b * b - 3.0 * c # 临界点判别 if inner 0.0: # 单调唯一实根 x, _, ok newton(f, df, bound) return [x] if ok else [] sq math.sqrt(inner) crit sorted([(-b - sq) / 3.0, (-b sq) / 3.0]) # 两个临界点 cuts [-bound] crit [bound] roots [] for i in range(3): # 三个单调区间 lo, hi cuts[i], cuts[i 1] flo, fhi f(lo), f(hi) if abs(flo) tol: roots.append(lo); continue if flo * fhi 0: # 同号区间内无根 continue # 有根用右端点做牛顿初值失败则回退二分 x, _, ok newton(f, df, hi) if not ok or not (lo - 1e-9 x hi 1e-9): lo_, hi_ lo, hi for _ in range(200): mid 0.5 * (lo_ hi_) if f(lo_) * f(mid) 0: hi_ mid else: lo_ mid x 0.5 * (lo_ hi_) roots.append(x) return sorted(roots)f用秦九韶Horner形式写比直接x**3 b*x**2 c*x d少两次乘方对需要迭代几百次的场景是实打实的提速。临界点公式由 3x²2bxc0 解出 x (−b ± √(b²−3c))/3注意这里是inner b² − 3c而不是判别式 b²−4ac因为二次项系数是 3 不是 1这一步很多人写错。牛顿法之外保留了二分回退是因为牛顿法在某些区间端点附近可能跳出当前区间收敛到别的根而二分法虽然慢但永远不会跑偏作为兜底手段非常合适。3.3 路径三先求一个实根再降阶三次方程必然至少有一个实根这是个很好用的性质。先用上面两种方法任意一种求出一个实根 r然后做多项式除法把三次降成二次剩下两个根用求根公式一步到位。综合除法的系数递推是设 f(x) ax³bx²cxd除以 (x−r) 后得到 ax² (arb)x (ar²brc)余项是 ar³br²crd理论上为 0实际会是个小残差。这一步的好处是把“求三个根”变成“求一个根 解一个二次方程”二次求根公式的数值稳定性大家都很熟了。降阶法最大的风险是误差累积如果 r 的精度只有 1e-8降阶后的二次方程系数误差会被放大剩下两个根可能只有 1e-5 的精度。我的经验是只在对精度要求不高、或者只需要那一个实根的场景用降阶法要全部根且要精度还是走 3.1 或 3.2。3.4 用韦达定理做回代校验不管是哪条路径根算出来之后都应该免费做一次校验成本几乎为零收益极大。首一化后的三个根必须满足r₁ r₂ r₃ −br₁r₂ r₂r₃ r₃r₁ cr₁r₂r₃ −ddef vieta_check(roots, b, c, d, tol1e-8): re [r[0] for r in roots] s1 sum(re) s2 re[0]*re[1] re[1]*re[2] re[2]*re[0] s3 re[0]*re[1]*re[2] err max(abs(s1 b), abs(s2 - c), abs(s3 d)) return err tol * max(1.0, abs(b), abs(c), abs(d)), err这里用对称多项式的相对误差做判据比逐个代入 f(x) 更好用——因为 f(x) 在重根附近对扰动不敏感本来就是零而韦达关系对误差是放大的更容易暴露问题。我在项目里就是靠这个校验抓出过一次系数换算函数的符号错误当时单测全过但韦达残差一直在 1e-3 量级下不去一看就是 q 的公式里2b³/27写成了b³/27。4. 参数与精度那些让结果翻车的细节4.1 立方根分支与复数抵消的处理卡尔达诺公式最大的数值隐患在 Δ 接近零的时候。当 Δ 是个很小的正数√Δ ≈ 0U 和 V 都近似等于 −q/2如果 q 本身很大那么 −q/2 √Δ 和 −q/2 − √Δ 两个数的相对差极小各自开立方之后的相对误差会被放大到原来的三分之一——听起来还好但如果后面还要做 uv 的相加两个几乎相等的数相加倒不会抵消同号相加是安全的真正危险的是 Δ0 走到复数域那条路实部相减会产生灾难性抵消。所以我的实现里坚持一条Δ0 一律走三角解法绝不允许进入复数分支。这条规则让代码少了一整类线上事故。反过来Δ0 时用立方根公式是安全的因为此时 U 和 V 一正一负或者都接近零但不会相互抵消开立方后符号保持相加是真正的加法而不是抵消。还有一个容易被忽略的点_cbrt的精度。abs(x) ** (1/3)在 x 很大或很小时误差会放大因为 1/3 在二进制里不能精确表示。要求高的时候可以用math.exp(math.log(abs(x))/3)或者直接调用numpy.cbrt后者内部做了优化精度更好。我在一个小量级的物理仿真里对比过**版本在 x≈1e-12 时的相对误差能到 1e-10 量级换numpy.cbrt之后降到 1e-16。4.2 迭代停止条件怎么定牛顿法的停止条件看似简单其实决定了结果质量。常见的三种写法停止条件写法适用与风险绝对步长abs(step) 1e-12根很大时过早停止根很小时永远停不下来相对步长abs(step) tol*abs(x)x 接近 0 时分母趋零可能死循环混合判据abs(step) tol*max(1, abs(x))推荐兼顾大小量级我一般用第三种tol 取 1e-14双精度机器精度约 2.2e-16留两个数量级的余量。另外一定要加最大迭代次数200 次足够了——牛顿法二次收敛正常情况下 5 到 8 次就到极限精度跑满 200 次说明初值选得很差或者在震荡这时候与其继续算不如回退二分。还有一个隐蔽的坑判断收敛不能用 f(x) 是否接近零。三次函数在极值点附近的斜率很小f(x) 可以很小但 x 离根还很远此时停止会得到一个大偏差的解。用步长判据最可靠因为步长直接反映 x 的移动量。4.3 重根与近重根的敏感性重根是三次方程求根里最麻烦的情况。当 Δ 严格等于零时公式给出二重根但浮点环境下 Δ 永远不会严格为零于是你会得到一个 Δ 1e-18 的“近重根”公式算出来是两个相差 1e-9 的根而真值应该是两个完全相同的根。这种伪分裂在后续计算里会被放大比如做根的敏感性分析或者判断根的个数时会误导结论。处理办法有两条。一是用容差合并近根算完之后检查任意两根的距离如果小于某个阈值比如 1e-8 乘以根的尺度就认为它们是重根取平均合并。二是用 f 与 f′ 的最大公因式判断重根三次方程有重根的充要条件是判别式为零工程上等价于 f 和 f′ 有公共根数值上可以检查 f 在两个临界点处的值是否接近零。我在代码里用的是第一条简单可靠。还有一种更阴险的情况叫“近重根”两个根距离 1e-4 但确实不同。这时候闭式解和迭代法都可能只给出一个根另一个被吞掉因为函数在两根之间的区域几乎是平的。解决办法是提高精度用高精度浮点或者用韦达关系反推已知一个根和根的和另一个根可以由 −b − r₁ 得到比直接在平面上找更稳。5. 常见问题速查与排查思路5.1 典型异常与对应处理现象可能原因处理方式acos 抛 domain error判别式接近零时 arg 略超 [−1,1]对 arg 做夹紧负系数开立方得到 nan用了x**(1/3)而非带符号的 cbrt手写_cbrt或换 numpy三个根里出现 NaNp0 且走了三角分支除以零p≥0 必须先判单调分支迭代 200 次不收敛初值落在极值点附近或导数接近零换区间端点做初值或直接二分韦达残差偏大系数换算公式错尤其 q用 f(−b/3)q 单测系数函数重根算成两个近根浮点判别式不可能严格为零用容差合并近根大系数下结果偏差首一化时除法的舍入误差累积先按最大系数做尺度归一化表格里最后一条值得单独说。当系数跨度很大比如 a1e-8、d1e8首一化之后系数能到 1e16平方之后直接溢出。解决办法是在首一化之前先把所有系数除以 max(|a|,|b|,|c|,|d|)只改变尺度不改变根但能避免中间量溢出。这个技巧我在处理实验数据拟合时救过好几次场。5.2 我自己踩过的几个坑第一个坑是早期的实现里我把判别式写成了q*q/4 p*p*p/27看起来没问题但当 p 是 −1e5 量级时p³ −1e15而 q²/4 可能只有 1e-3两者相加时小的那个直接被吞掉。这不是 bug 是浮点特性只能靠量纲分析和容差设计规避。后来我改成先算tp p/3用tp**3计算虽然还是同一个量但避免了 p³ 的中间溢出。第二个坑是排序。有一次几何求交的代码偶发算错交点查了两天才发现是三角解法给出的三个根顺序不固定取决于 arccos 的分支调用方按“第一个实根”去取结果拿到了最大的那个后面的判断全反了。从此我的所有求根函数都强制排序并且把这条写进了接口文档。第三个坑是复数根的处理。调用方如果只关心实根闭式解返回的复根会被他直接丢掉但如果遇到 Δ 略大于零而真实情况是三个实根系数本身有测量误差丢掉的两个“复根”其实对应真实存在的两个实根只是实部相同、虚部很小。这种情况在实验数据处理里非常常见我的处理经验是当复根的虚部相对于实部小于 1e-6 时把它当作实根输出并标记为“疑似重根”让人工或下游算法决定怎么用。6. 应用场景这玩意儿在工程里到底用来干嘛6.1 曲线拟合与几何求交三次方程在实际项目里出现频率最高的场景是求交。直线与三次贝塞尔曲线的交点、两条参数曲线的交点、圆与抛物线的交点参数化之后都会落到三次方程上。以贝塞尔曲线为例P(t) (1−t)³P₀ 3(1−t)²tP₁ 3(1−t)t²P₂ t³P₃要和一条直线求交把直线方程代入并整理成 t 的三次方程求 t∈[0,1] 内的根。这种场景对求解器的要求很特别只需要区间内的根区间外的根一律丢弃。所以我在这一类项目里不会用完整的闭式解而是用 3.2 的分段牛顿法把搜索范围限制在 [0,1]再用单调性判断最多有几个根效率比全解高一个数量级。另一个经验是 t 的定义域天然有界这让根的界估计变得非常容易初值可以直接取端点收敛一般不超过 6 次迭代。在数据拟合里三次多项式回归是最常用的非线性拟合模型之一拟合出系数之后经常需要求极值点或者零点这时候又是解方程——极值点解二次方程零点是三次方程。模具设计、材料应力曲线、经济学的供需均衡模型都会在这里碰到同一个求解器。6.2 控制系统与状态方程的稳定性判定在控制领域系统的特征方程常常是三次的。判断一个三阶系统是否稳定本质上就是判断特征方程的三个根是否都落在复平面的左半部分也就是所有根的实部都小于零。这种场景不关心根的具体数值只关心符号分布所以用劳斯判据或者直接用判别式分析更快。但如果要做参数整定需要知道根随参数变化的轨迹根轨迹那就必须真求根了。工程实践里的做法是遍历参数网格每个格点解一次三次方程把所有根画在复平面上。这种场景对单次求解速度要求不高但对稳定性要求极高——一旦某个格点因为数值问题返回了错误的根整条根轨迹就会出现跳变图上会出现一根莫名其妙的断线。我的经验是这种批量计算要用统一的分支逻辑绝对不要在循环里根据结果动态切换算法否则相邻格点用了不同路径结果的系统性偏差会画出假的轨迹分叉。三次特征方程还会出现在电路分析里三个储能元件的网络其暂态响应的时间常数就是特征方程三个根的倒数。这个场景里根必须是负数否则系统自激发散所以求解之后直接检查实部符号就行比完整的复数运算简单很多。6.3 物理与工程中的临界条件计算三次方程在物理里最经典的出场是运动学问题。抛体运动中求“何时达到某个高度”如果加上线性空气阻力方程会变成指数形式但如果加上与速度平方成正比的阻力位置方程对时间求逆就会落到三次上。结构力学里算梁的最大挠度位置挠度方程是三次的求导得二次但要算某个指定挠度对应的位置又要解三次。热力学里某些状态方程比如三次型的状态方程求体积也是同一个套路。这类场景的共同特点是物理量有明确单位量纲一致性是查错的第一手段。如果算出来的根数量级明显不对先别怀疑算法先检查系数换算的单位——我见过把毫米当米代入、把度当弧度代入的情况算出来的根看着“很合理”实际差了三个数量级。所以在这种场景下我习惯在求解之前先把所有系数归一化到 O(1) 量级解出来之后再乘回去这样物理直觉和数值稳定性都能兼顾。另外提一句这个求解器在图像处理里也有用武之地。三次样条插值的分段函数求反函数、伽马校正的某些反解、颜色空间转换里的多项式变换都会出现三次方程。这类应用对精度的要求通常在 1e-6 量级用迭代法配合严格的最大迭代次数就够了不需要动用闭式解。7. 一个可以直接抄走的组合方案把前面所有内容压成一个可落地的策略我平时就是这么用的系数先做尺度归一化除以最大绝对值系数首一化消二次项得到 p、q判断 p ≥ 0单调用牛顿法在柯西界处起迭代只输出一个实根判断 p 0算判别式Δ eps 用闭式解Δ −eps 用三角解法落中间用重根分支算完统一做韦达校验残差超阈值就切到分段牛顿法重算近根按容差合并标记疑似重根输出前统一排序。这套组合在实际项目里跑了大半年覆盖了从 1e-6 到 1e8 量级的系数量级没出过精度问题。要说唯一的代价就是代码比“一行公式”长了不少但换来的是任何一组系数都能给出可解释的结果——这在工程上比省几行代码重要得多。最后再分享一个小技巧如果你只是想快速验证一组系数对不对不要急着写整个求解器先把系数代进f(-b/3)看是否等于 q再把三个根代回韦达关系看残差。这两个检查五秒钟就能写完能提前拦掉绝大多数符号和换算错误。我现在的习惯是每写一个新的多项式处理函数都先把这两个断言加上再往下写逻辑省下来的调试时间远比写断言的时间多。
