不使用sympy求解三次方程，可否通过numpy或纯Python实现？

错误原因排查 #

你用numpy得到的结果不对，核心是系数构造错误：

• 原方程中num = 30 + 4.44 = 34.44，方程形式为 0.0509x³ + 0.0192x² + 3.68x - num = 0，按x降序排列的系数最后一位是常数项-34.44
• 你写的--4.44等价于+4.44，且完全没有用到计算得到的num值，相当于求解的是完全不同的三次方程，结果自然不一致。

正确numpy实现代码 #

import numpy as np
xp = 30
num = xp + 4.44
# 系数按x³、x²、x、常数项的顺序排列
coeff = [0.0509, 0.0192, 3.68, -num]
roots = np.roots(coeff)
print(roots)

运行输出结果和sympy的计算结果一致：

[ 6.07118098+0.j        -3.2241956 +10.05248912j -3.2241956 -10.05248912j]

无第三方库实现（基于卡尔达诺公式） #

如果不想依赖任何第三方库，可以使用Python标准库cmath手动实现三次方程求根逻辑：

import cmath

def solve_cubic(a, b, c, d):
    # 归一化最高次项系数为1
    b /= a
    c /= a
    d /= a
    # 消去二次项：做变量替换x = y - b/3
    p = c - b**2 / 3
    q = d - b*c/3 + 2*b**3 / 27
    # 计算三次方程判别式
    delta = (q/2)**2 + (p/3)**3
    # 求解三次根
    u = ((-q/2) + cmath.sqrt(delta)) ** (1/3)
    v = -p/(3*u) if u != 0 else ((-q/2) - cmath.sqrt(delta)) ** (1/3)
    # 生成三个根（复数单位根转换）
    omega = complex(-0.5, cmath.sqrt(3)/2)
    roots = [
        u + v - b/3,
        u*omega + v*omega**2 - b/3,
        u*omega**2 + v*omega - b/3
    ]
    return roots

# 代入参数计算
xp = 30
num = xp + 4.44
roots = solve_cubic(0.0509, 0.0192, 3.68, -num)
for r in roots:
    print(r)

运行后输出的根和sympy、numpy的计算结果一致。

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