c - 找到（三次）多项式实根的简单方法是什么？

Question

这对我来说似乎是一个显而易见的问题，但我在 SO 的任何地方都找不到它。我有一个三次多项式，我需要找到函数的实根。这样做的方法是什么？

我找到了几个三次函数根的封闭形式公式，但它们都使用复数或大量测角函数，我不喜欢它们（也不知道该选择哪一个）。

我需要一些简单的东西；越快越好；而且我知道我最终将需要求解高阶多项式，因此拥有一个数值求解器可能也会有所帮助。我知道我可以使用一些库来为我做艰苦的工作，但可以说我想把它作为一个练习来做。

我正在用 C 编码，所以import magic_poly_solver请不要。

额外的问题：我如何只找到给定区间内的根？

score 10 · Accepted Answer

对于三次多项式，有封闭形式的解，但它们并不是特别适合数值微积分。

对于三次情况，我会执行以下操作：任何三次多项式至少有一个实根，您可以使用牛顿法轻松找到它。然后，您使用通货紧缩来求解剩余的二次多项式，请参阅我的答案以了解如何正确执行后一步。

需要注意的一点是：如果判别式接近于零，就会有一个数值上的多重实根，牛顿的方法将惨遭失败。此外，由于在根附近，多项式类似于 (x - x0)^2，因此您将丢失一半有效数字（因为一旦 x - x0 < sqrt(epsilon，P(x) 将 < epsilon ））。因此，您可能希望排除这种情况并在这种特殊情况下使用封闭形式的解决方案，或者求解导数多项式。

如果要在给定区间内求根，请查看Sturm 定理。

用于通用多项式求解的更通用（复杂）算法是Jenkins-Traub 算法。这显然是矫枉过正，但它适用于立方。通常，您使用第三方实现。

既然你做 C，那么使用GSL肯定是你最好的选择。

另一种通用方法是使用例如找到伴随矩阵的特征值。平衡 QR 分解，或简化为 Householder 形式。这是 GSL 采用的方法。

score 3 · Accepted Answer

为了用简单的 C 代码求解三次方程，我发现QBC著名数字专家教授 William Kahan 的求解器足够稳健、相当快速且相当准确：

威廉·卡汉（William Kahan），“求解一个真正的三次方程。” PAM-352，加州大学伯克利分校纯数学和应用数学中心。1986 年 11 月 10 日。（在线，在线）

这使用基于导数的迭代方法来找到实根，基于此简化为二次方程，最后使用数值稳健的二次方程求解器来找到剩余的两个根。通常，迭代求解器需要大约五到十次迭代才能收敛到结果。通过明智地使用融合乘加(FMA) 运算，这两种求解器都可以提高精度和性能，这些运算在 ISO C99 中通过fma()标准数学函数提供。

对二次求解器的准确性至关重要的是判别式的计算。为此，我根据最近的研究使用以下代码：

/* Compute B*B - A*C, accurately 

   Claude-Pierre Jeannerod, Nicolas Louvet, and Jean-Michel Muller, 
   "Further Analysis of Kahan's Algorithm for the Accurate Computation 
   of 2x2 Determinants". Mathematics of Computation, Vol. 82, No. 284, 
   Oct. 2013, pp. 2245-2264

   https://www.ams.org/journals/mcom/2013-82-284/S0025-5718-2013-02679-8/S0025-5718-2013-02679-8.pdf
*/
double DISC (double A, double B, double C)
{
    double w = C * A;
    double e = fma (-C, A, w);
    double f = fma (B, B, -w);
    double r = f + e;
    return r;
}

使用双精度算术，Kahan 的求解器不能总是产生精确到双精度的结果。Kahan 的论文中提供的一个测试用例说明了为什么会这样：

658x³ - 190125x² + 18311811x - 587898164

使用任意精度数学库，我们发现这个三次方程的根如下：

96.229639346592182_18...
96.357064825184152_07... ± i * 0.069749752043689625_43...

QBC使用双精度算术计算根为

96.2296393 50445893
96.35706482 3257289 ± i * 0.0697497 48521837268

原因是围绕实根的函数评估在计算的函数值中存在高达 60% 的误差，从而阻止了迭代求解器更接近根。通过更改函数和导数评估以使用双双计算进行中间计算（计算成本很高），我们可以解决这个问题。

/* Data type for double-double computation */
typedef struct {
    double l;  // low / tail
    double h;  // high / head
} dbldbl;

dbldbl make_dbldbl (double head, double tail);
double get_dbldbl_head (dbldbl a);
double get_dbldbl_tail (dbldbl a);
dbldbl add_dbldbl (dbldbl a, dbldbl b);
dbldbl mul_dbldbl (dbldbl a, dbldbl b);

void EVAL (double X, double A, double B, double C, double D, 
           double * restrict Q, double * restrict Qprime, 
           double * restrict B1, double * restrict C2)
{
#if USE_DBLDBL_EVAL
    dbldbl AA, BB, CC, DD, XX, AX, TT, UU;
    AA = make_dbldbl (A, 0);
    BB = make_dbldbl (B, 0);
    CC = make_dbldbl (C, 0);
    DD = make_dbldbl (D, 0);
    XX = make_dbldbl (X, 0);
    AX = mul_dbldbl (AA, XX);
    TT = add_dbldbl (AX, BB);
    *B1 = get_dbldbl_head (TT) + get_dbldbl_tail(TT);
    UU = add_dbldbl (mul_dbldbl (TT, XX), CC);
    *C2 = get_dbldbl_head (UU) + get_dbldbl_tail(UU);
    TT = add_dbldbl (mul_dbldbl (add_dbldbl (AX, TT), XX), UU);
    *Qprime = get_dbldbl_head (TT) + get_dbldbl_tail(TT);
    UU = add_dbldbl (mul_dbldbl (UU, XX), DD);
    *Q = get_dbldbl_head (UU) + get_dbldbl_tail(UU);
#else // USE_DBLDBL_EVAL
    *B1 = fma (A, X, B); 
    *C2 = fma (*B1, X, C);
    *Qprime = fma (fma (A, X, *B1), X, *C2);
    *Q = fma (*C2, X, D);
#endif // USE_DBLDBL_EVAL
}

/* Construct new dbldbl number. |tail| must be <= 0.5 ulp of |head| */
dbldbl make_dbldbl (double head, double tail)
{
    dbldbl z;
    z.l = tail;
    z.h = head;
    return z;
}

/* Return the head of a double-double number */
double get_dbldbl_head (dbldbl a)
{
    return a.h;
}

/* Return the tail of a double-double number */
double get_dbldbl_tail (dbldbl a)
{
    return a.l;
}

/* Add two dbldbl numbers */
dbldbl add_dbldbl (dbldbl a, dbldbl b)
{
    dbldbl z;
    double e, q, r, s, t, u;

    /* Andrew Thall, "Extended-Precision Floating-Point Numbers for GPU 
       Computation." 2006. http://andrewthall.org/papers/df64_qf128.pdf
    */
    q = a.h + b.h;
    r = q - a.h;
    t = (a.h + (r - q)) + (b.h - r);
    s = a.l + b.l;
    r = s - a.l;
    u = (a.l + (r - s)) + (b.l - r);
    t = t + s;
    s = q + t;
    t = (q - s) + t;
    t = t + u;
    z.h = e = s + t;
    z.l = (s - e) + t;

    /* For result of zero or infinity, ensure that tail equals head */
    if (isinf (s)) {
        z.h = s;
        z.l = s;
    }
    if (z.h == 0) {
        z.l = z.h;
    }

    return z;
}

/* Multiply two dbldbl numbers */
dbldbl mul_dbldbl (dbldbl a, dbldbl b)
{
    dbldbl z;
    double e, s, t;

    s = a.h * b.h;    
    t = fma (a.h, b.h, -s);
    t = fma (a.l, b.l, t);
    t = fma (a.h, b.l, t);
    t = fma (a.l, b.h, t);
    z.h = e = s + t;
    z.l = (s - e) + t;

    /* For result of zero or infinity, ensure that tail equals head */
    if (isinf (s)) {
        z.h = s;
        z.l = s;
    }
    if (z.h == 0) {
        z.l = z.h;
    }

    return z;
}

使用更准确的函数和导数评估计算的根是：

96.22963934659218 0
96.35706482518415 3 ± i * 0.06974975204 5672006

虽然实部现在精确到双精度范围内，但虚部仍然关闭。这样做的原因是，在这种情况下，二次方程对系数的微小差异很敏感。它们中的任何一个中的一个 ulp 错误都可能导致虚部的差异约为 10 ^-11。这可以通过将系数表示为高于双精度并在二次求解器中使用更高精度的计算来解决。

score 2 · Accepted Answer

如果您不想使用封闭的解决方案（或期望更大阶的多项式），最明显的方法是使用牛顿法计算近似根。

不幸的是，无法确定迭代时将获得哪些根，尽管它取决于起始值。

也见这里。

score 2 · Accepted Answer

2

请参阅D Herbison-Evans 在Graphics Gems V上发表的求解图形的四次和三次。

于 2011-02-05T13:22:08.663 回答

score 0 · Accepted Answer

/*******************************************************************************
 * FindCubicRoots solves:
 *      coeff[3] * x^3 + coeff[2] * x^2 + coeff[1] * x + coeff[0] = 0
 *  returns:
 *      3 - 3 real roots
 *      1 - 1 real root (2 complex conjugate)
 *******************************************************************************/

int FindCubicRoots(const FLOAT coeff[4], FLOAT x[3]);

http://www.realitypixels.com/turk/opensource/index.html#CubicRoots

c - 找到（三次）多项式实根的简单方法是什么？

5 回答 5

Related

Reference