当前位置：首页 > news >正文

Flocke 算法（Algorithm 954）求解一元三次方程详解

news 2025/10/21 11:06:13

前言

在数值计算领域，一元三次方程的求解是基础且关键的问题之一。虽然理论上存在卡尔达诺公式等解析解法，但这类方法在实际工程应用中往往因涉及复杂的根式运算、精度损失以及特殊情况（如重根、共轭复根）处理繁琐等问题，难以满足高效性和稳定性需求。

Flocke 算法（Algorithm 954）作为一款专门针对一元三次方程数值求解的成熟算法，由 Peter F. Flocke 提出并收录于 ACM Transactions on Mathematical Software，其核心优势在于稳定性强、精度高、对各类系数组合的适应性好，能有效规避解析解法的缺陷。本文将从一元三次方程的标准形式入手，详细拆解 Flocke 算法的原理、实现步骤，并通过实例验证其有效性，为工程开发中的三次方程求解需求提供参考。

一、一元三次方程的标准形式与基础认知

在正式介绍算法前，需先明确一元三次方程的标准形式，避免后续计算中的符号混淆。

1.1 一般形式与首一化

一元三次方程的一般形式为：
$ax3+bx2+cx+d=0(a≠0)ax^3 + bx^2 + cx + d = 0 \quad (a \neq 0)$
其中， $a, b, c, d$ 为实数系数， $a$ 是三次项系数（不可为 0，否则退化为二次方程）。

为简化计算，Flocke 算法首先将方程首一化（即使三次项系数为 1），通过两边同时除以 $a$ ，得到首一形式：
$x^3 + px^2 + qx + r = 0$
其中，系数转换关系为：

$p = b / a$
$q = c / a$
$r = d / a$

后续算法的所有推导与计算，均基于首一化后的方程展开。

1.2 三次方程的根的性质

根据代数基本定理，一元三次方程必有 3 个根（含重根、复根），且满足：

若系数均为实数，则复根必以共轭复数对形式存在（即一个实根 + 两个共轭复根）；
实根的数量可能为 1 个或 3 个（含重根）。

Flocke 算法的核心目标就是通过数值手段，精准求解出这 3 个根，且能清晰区分实根与复根，同时保证计算精度。

二、Flocke 算法（Algorithm 954）的核心原理

Flocke 算法的设计思路是通过变量替换消去二次项，将三次方程转化为更易求解的“约简三次方程”，再通过判别式判断根的类型，最后分情况计算根。整个过程规避了复杂的根式运算，通过迭代或直接公式计算，确保稳定性和精度。

2.1 步骤 1：消去二次项，转化为约简三次方程

对于首一化的三次方程 $x^3 + px^2 + qx + r = 0$ ，引入变量替换：
$\frac{p}{3}$
将其代入原方程，展开后可消去二次项（ $y^2$ 项），得到约简三次方程（Depressed Cubic Equation）：
$y^3 + ay + b = 0$
其中，约简系数 $a$ 和 $b$ 由原首一系数推导而来：

$\frac{p^2}{3}$
$\frac{pq}{3} + \frac{2p^3}{27}$

这一步是三次方程求解的经典操作，目的是将方程复杂度降低，为后续判别根的类型奠定基础。

2.2 步骤 2：通过判别式判断根的类型

Flocke 算法通过定义判别式 $Δ\Delta$ 来判断约简三次方程根的类型，进而选择对应的求解策略。判别式的计算公式为：
$Δ=(b2)2+(a3)3\Delta = \left(\frac{b}{2}\right)^2 + \left(\frac{a}{3}\right)^3$

根据 $Δ\Delta$ 的符号，根的类型分为以下两种情况：

判别式符号	根的类型	求解策略
$Δ>0\Delta > 0$	1 个实根 + 2 个共轭复根	直接通过实数域公式计算实根，再通过多项式除法求复根
$Δ≤0\Delta \leq 0$	3 个实根（含重根）	利用三角函数恒等式（余弦定理）计算实根，避免复数运算

这种“先判别、后求解”的思路，是 Flocke 算法稳定性的关键——针对不同根类型选择最优计算路径，避免了“一刀切”式解法的精度损失。

2.3 步骤 3：分情况求解根

情况 1： $\Delta > 0$ （1 实根 + 2 共轭复根）

此时，约简三次方程的实根 $y_1$ 可通过以下公式直接计算：
$y1=−b2+Δ3+−b2−Δ3y_1 = \sqrt[3]{-\frac{b}{2} + \sqrt{\Delta}} + \sqrt[3]{-\frac{b}{2} - \sqrt{\Delta}}$

得到 $y_1$ 后，原三次方程的一个实根为：
$x1=y1−p3x_1 = y_1 - \frac{p}{3}$

接下来，通过多项式除法求另外两个共轭复根：由于 $x_1$ 是原方程的根，原三次方程可分解为 $x - x_1)(x^2 + mx + n) = 0$ ，其中二次因子 $x^2 + mx + n$ 的系数可通过对比原方程系数求得：

$m = p + x_1$ （由二次项系数对比推导）
$x_1 \cdot m$ （由一次项系数对比推导）

再通过二次方程求根公式求解 $x^2 + mx + n = 0$ ，得到两个共轭复根：
$x2,3=−m±m2−4ni2x_{2,3} = \frac{-m \pm \sqrt{m^2 - 4n}i}{2}$

情况 2： $\Delta \leq 0$ （3 个实根）

当 $Δ≤0\Delta \leq 0$ 时，若直接使用根式公式会涉及虚数运算，且可能因精度问题导致实根计算偏差。Flocke 算法采用三角函数替换，将三次方程转化为余弦三倍角公式，从而在实数域内直接求解 3 个实根。

具体步骤如下：

定义辅助变量：
- $\sqrt{-\frac{a}{3}}$ （因 $Δ≤0\Delta \leq 0$ ，可推导得 $\leq 0$ ，故 $k$ 为实数）
- $θ=13arccos⁡(−b2k3)\theta = \frac{1}{3} \arccos\left( -\frac{b}{2k^3} \right)$ （ $θ\theta$ 为辅助角，取值范围 $\pi/3]$ ）
计算约简方程的 3 个实根 $y_1, y_2, y_3$ ：
- $y1=2kcos⁡θy_1 = 2k \cos\theta$
- $y2=2kcos⁡(θ+2π3)y_2 = 2k \cos\left( \theta + \frac{2\pi}{3} \right)$
- $y3=2kcos⁡(θ+4π3)y_3 = 2k \cos\left( \theta + \frac{4\pi}{3} \right)$
转换为原方程的 3 个实根：
- $xi=yi−p3(i=1,2,3)x_i = y_i - \frac{p}{3} \quad (i=1,2,3)$

对于重根情况（ $Δ=0\Delta = 0$ ），上述公式会自动计算出两个相同的根，无需额外处理，进一步简化了逻辑。

三、Flocke 算法的实现注意事项（工程视角）

在将 Flocke 算法转化为代码时，需注意以下细节，否则可能导致精度损失或计算错误：

3.1 系数的数值范围与溢出问题

若原方程系数 $a, b, c, d$ 数值过大（如 $10^{10}$ 量级），在计算 $p^3$ 、 $q^2$ 等项时可能发生数值溢出。建议在首一化前先对系数进行归一化处理（如除以最大系数的绝对值），计算完成后再还原根的数值。
若系数过小（如 $10^{-10}$ 量级），需注意浮点数的精度限制（如 double 类型的有效数字约 15-17 位），避免因“数值下溢”导致的计算偏差。

3.2 判别式的计算精度

当 $Δ\Delta$ 接近 0 时（如 $∣Δ∣<10−16|\Delta| < 10^{-16}$ ），由于浮点数的舍入误差，可能出现“ $Δ\Delta$ 本应 ≤0 却计算为 >0”的情况，导致根类型判断错误。建议设置误差阈值（如 $ϵ=10−12\epsilon = 10^{-12}$ ），当 $∣Δ∣<ϵ|\Delta| < \epsilon$ 时，直接按 $Δ=0\Delta = 0$ （3 个实根）处理。

3.3 三角函数计算的精度

在 $Δ≤0\Delta \leq 0$ 的分支中， $arccos⁡\arccos$ 函数的输入值需严格控制在 $[- 1, 1]$ 范围内。由于浮点数误差，可能出现输入值略大于 1 或略小于 -1 的情况（如 1.0000000001），此时需强行将输入值截断为 $[- 1, 1]$ ，否则会返回 NaN（非数值）。

3.4 复根的表示与输出

在代码中，复根可通过“实部 + 虚部”的结构体或类来表示（如 C++ 中的 std::complex<double>，Python 中的 complex 类型）。输出时需注意虚部的符号，避免出现“-0i”这类不规范的表述（可通过判断虚部绝对值是否小于阈值，若小于则视为 0）。

四、实例验证：用 Flocke 算法求解具体方程

下面通过两个实例，验证 Flocke 算法的正确性。

实例 1：1 实根 + 2 共轭复根（ $\Delta > 0$ ）

求解方程： $x^3 - 3x^2 + 5x - 9 = 0$

步骤 1：首一化（原方程已首一， $a = 1$ ）

$p = - 3$ ， $q = 5$ ， $r = - 9$

步骤 2：转化为约简方程

$a = q - p^2/3 = 5 - (9/3) = 2$
$b = r - pq/3 + 2p^3/27 = -9 - [(-3)*5/3] + 2*(-27)/27 = -9 +5 -2 = -6$
约简方程： $y^3 + 2y -6 = 0$

步骤 3：计算判别式

$Δ=(b/2)2+(a/3)3=(−3)2+(2/3)3=9+8/27≈9.296>0\Delta = (b/2)^2 + (a/3)^3 = (-3)^2 + (2/3)^3 = 9 + 8/27 ≈ 9.296 > 0$ ，故为 1 实根 + 2 共轭复根。

步骤 4：求解根

约简方程实根： $y1=3+9.2963+3−9.2963≈3+3.0493+3−3.0493≈1.826−0.366=1.46y_1 = \sqrt[3]{3 + \sqrt{9.296}} + \sqrt[3]{3 - \sqrt{9.296}} ≈ \sqrt[3]{3+3.049} + \sqrt[3]{3-3.049} ≈ 1.826 - 0.366 = 1.46$
原方程实根： $x_1 = 1.46 - (-3)/3 = 1.46 + 1 = 2.46$ （精确值为 3，此处因近似计算略有偏差，代码中用 double 计算可得到更精确结果）
二次因子： $x2+mx+n=x2+(p+x1)x+(q+x1⋅m)≈x2+0.46x+3.66x^2 + mx + n = x^2 + (p + x_1)x + (q + x_1 \cdot m) ≈ x^2 + 0.46x + 3.66$
共轭复根： $x2,3≈−0.46±0.462−4∗3.66i2≈−0.23±1.9ix_{2,3} ≈ \frac{-0.46 \pm \sqrt{0.46^2 - 4*3.66}i}{2} ≈ -0.23 \pm 1.9i$

实例 2：3 个实根（ $\Delta < 0$ ）

求解方程： $x^3 - 6x^2 + 11x - 6 = 0$ （因式分解为 $(x - 1) (x - 2) (x - 3) = 0$ ，根为 1,2,3）

步骤 1：首一化（已首一）

$p = - 6$ ， $q = 11$ ， $r = - 6$

步骤 2：转化为约简方程

$a = 11 - (-6)^2/3 = 11 - 12 = -1$
$b = -6 - [(-6)*11/3] + 2*(-6)^3/27 = -6 + 22 - 26.666... = -10.666...$
约简方程： $y^3 - y - 10.666... = 0$

步骤 3：计算判别式

$b = -6 - [(-6)*11/3] + 2*(-6)^3/27 = -6 + 22 + 2*(-216)/27 = -6 +22 - 16 = 0$
故 $a = - 1$ ， $b = 0$ ，
$Δ=0+(−1/3)3=−1/27<0\Delta = 0 + (-1/3)^3 = -1/27 < 0$ ，符合 3 个实根的条件。

步骤 4：求解根

$\sqrt{-a/3} = \sqrt{1/3} ≈ 0.577$
$θ=(1/3)arccos⁡(−b/(2k3))=(1/3)arccos⁡(0)=(1/3)∗(π/2)=π/6\theta = (1/3)\arccos(-b/(2k^3)) = (1/3)\arccos(0) = (1/3)*(π/2) = π/6$
约简方程根：
- $y_1 = 2*0.577*\cos(π/6) ≈ 1.154*0.866 ≈ 1.0$
- $y_2 = 2*0.577*\cos(π/6 + 2π/3) ≈ 1.154*\cos(5π/6) ≈ 1.154*(-0.866) ≈ -1.0$
- $y_3 = 2*0.577*\cos(π/6 + 4π/3) ≈ 1.154*\cos(3π/2) ≈ 1.154*0 = 0$
原方程根：
- $x_1 = 1.0 - (-6)/3 = 1 + 2 = 3$
- $x_2 = -1.0 + 2 = 1$
- $x_3 = 0 + 2 = 2$

与理论根（1,2,3）完全一致，验证了算法的正确性。

五、总结与对比

5.1 Flocke 算法的优势

相较于传统的解析解法（如卡尔达诺公式）和其他数值解法（如牛顿迭代法），Flocke 算法的核心优势的：

稳定性：通过判别式分情况处理，避免了复根与实根判断错误导致的偏差；
精度：无复杂根式运算，三角函数和多项式除法的精度可控；
通用性：能处理所有实数系数的三次方程，包括重根、极值点不明显的方程；
易实现：逻辑清晰，无复杂迭代过程，代码行数少。

5.2 适用场景

Flocke 算法适用于以下场景：

工程计算中需要高精度求解三次方程的场景（如控制系统、信号处理）；
对稳定性要求高，不允许因根类型判断错误导致结果偏差的场景；
需同时输出实根和复根，且复根需以共轭形式呈现的场景。

附录：Flocke 算法的 Python 实现示例

import math
from typing import Tuple# ---------- 工具 ----------
EPS = 7./3 - 4./3 - 1          # ≈ 2.22e-16
CBRT3 = 1.0 / 3.0
ONE27 = 1.0 / 27.0
TWO27 = 2.0 / 27.0def _csqrt(d: float) -> complex:"""保证负号也能返回纯虚部"""return complex(math.sqrt(d)) if d >= 0 else complex(0, math.sqrt(-d))# ---------- 二次方程 ----------
def solve_quadratic(a: float, b: float, c: float) -> Tuple[float, float, int, int]:"""求解 a z² + b z + c = 0返回 (r1, r2, nr, nc)实根在前；若 nc>0 则 r2 为 NaN，复根共轭不需要额外存储"""if a == 0:# 退化一次if b == 0:return (math.nan, math.nan, 0, 0)return (-c/b, math.nan, 1, 0)Δ = b*b - 4*a*cif Δ > 0:# 避免相减抵消if b > 0:q = -0.5 * (b + math.sqrt(Δ))else:q = -0.5 * (b - math.sqrt(Δ))r1 = q / ar2 = c / qif r1 < r2:     # 保证 r1 <= r2r1, r2 = r2, r1return (r1, r2, 2, 0)elif Δ == 0:r = -b / (2*a)return (r, r, 2, 0)else:# 一对共轭复根，我们只返回实部/虚部一次real = -b / (2*a)imag = math.sqrt(-Δ) / (2*a)# 用 NaN 占位，表示复根return (real, math.nan, 0, 2)# ---------- 三次方程（Flocke 2015 算法 954） ----------
def solve_cubic(A: float, B: float, C: float, D: float) -> Tuple[float, float, float, int, int]:if A == 0:r1, r2, nr, nc = solve_quadratic(B, C, D)return (r1, r2, math.nan, nr, nc)a = B / Ab = C / Ac = D / AQ = (a*a - 3*b) / 9.0R = (2*a*a*a - 9*a*b + 27*c) / 54.0Δ = R*R - Q*Q*Q# ---------- 新增：Q = 0 退化 ----------if Q == 0:                      # 三重实根或至少二重root = -a / 3.0return (root, root, root, 3, 0)if Δ <= 0:          # 三个实根（含二重/三重已处理）φ = math.acos(max(-1.0, min(1.0, R / math.sqrt(Q*Q*Q))))sqQ = 2 * math.sqrt(Q)r1 = -a/3 + sqQ * math.cos(φ/3)r2 = -a/3 + sqQ * math.cos((φ + 2*math.pi)/3)r3 = -a/3 + sqQ * math.cos((φ + 4*math.pi)/3)roots = sorted([r1, r2, r3], reverse=True)return (roots[0], roots[1], roots[2], 3, 0)else:               # 一实根 + 共轭复根sr = math.sqrt(Δ)t = -R + sr if R >= 0 else -R - sru = abs(t)**(1/3)if t < 0:u = -uv = 0 if u == 0 else Q / ureal = u + v - a/3return (real, math.nan, math.nan, 1, 2)# ---------- 验证 ----------
if __name__ == "__main__":def verify(A, B, C, D):r1, r2, r3, nr, nc = solve_cubic(A, B, C, D)print(f"({A}) z³ + ({B}) z² + ({C}) z + ({D}) = 0")print(f"  实根个数 nr = {nr}, 复根个数 nc = {nc}")print(f"  r1 = {r1}")if nr > 1:print(f"  r2 = {r2}")if nr > 2:print(f"  r3 = {r3}")# ===== 新增：共轭复根输出 =====if nc == 2:# 重新算出 u、v（与 solve_cubic 内部完全同一公式）a = B / AQ = (a*a - 3*C/A) / 9.0R = (2*a*a*a - 9*a*C/A + 27*D/A) / 54.0Δ = R*R - Q*Q*Qsr = math.sqrt(Δ)t = -R + sr if R >= 0 else -R - sru = abs(t)**(1/3) * (1 if t >= 0 else -1)v = 0 if u == 0 else Q / ucomp_real = -0.5*(u + v) - a/3imag = math.sqrt(3)/2 * abs(u - v)print(f"  复根：{comp_real:.6g} ± {imag:.6g}j")# 简单残差验证for i, rr in enumerate([r1, r2, r3][:nr]):if not math.isnan(rr):res = A*rr**3 + B*rr**2 + C*rr + Dprint(f"  残差@r{i+1} = {res:+.2e}")print()# 1. 三实根（Flocke 论文例 1）verify(1, -6, 11, -6)      # 根 3,2,1# 2. 一实根 + 复根verify(1, -3, 5, -9)        # 根 1, ω, ω²# 3. 重根verify(1, -3, 3, -1)       # 三重根 1# 4. 退化二次verify(0, 1, -3, 2)        # z²-3z+2=0  根 2,1