多项式算法高级 C++

快速稳定的多项式求根器 - 第四部分

Henrik Vestermark

0/5 (0投票)

2023年11月29日

CPOL

15分钟阅读

3483

Laguerre 方法用于稳定高效的多项式求根。

引言

在第一篇文章（第一部分）中，我们开发了一种基于 Newton 方法的高效鲁棒的多项式求根器，专门针对复系数多项式。在第二部分中，我们详细阐述了如何处理实系数多项式。在第三部分中，我们探讨了如何使用相同的框架实现高阶方法（Halley），并讨论了与标准 Newton 方法相比，使用高阶方法是否能获得优势。在第四部分中，我们将介绍 Laguerre 方法，在第五部分中，我们将深入研究其中一种同时求根方法，例如 Durand-Kerner 方法，这是一种完全不同的方法，但我们将看到如何将其相对容易地映射到第一、二、三部分建立的现有框架中。

为什么选择 Laguerre 方法？

在撰写本文之前，我仔细考虑了选择一种与系列前几部分不同的方法。我选择 Laguerre 方法有几个令人信服的理由，特别是与其他求根方法相比，它被低估的效率。

Laguerre 方法在处理高次多项式方面表现出色，使其成为一个绝佳的选择。与许多在多项式次数增加时在精度和收敛性方面遇到困难的其他方法不同，Laguerre 方法即使对于复杂的高次多项式，也能保持鲁棒和有效。

其卓越的收敛性尤其值得称道，特别是当初始猜测值 reasonably close to a true root 时。这种快速收敛部分归因于该方法在其迭代公式中同时使用了多项式的一阶和二阶导数。

此外，Laguerre 方法在处理复数根方面表现出色，这一点尤为突出。这是一个关键特性，因为许多高次多项式都具有复数根。它能巧妙地找出实根和复根。

更重要的是，与其他方法在初始猜测值不精确时可能失败或收敛到错误根不同，Laguerre 方法具有强大的全局收敛特性。这意味着它更有可能从更广泛的起始点收敛到某个根。

在实际应用中，尤其是在多项式方程普遍存在的工程和物理科学领域，Laguerre 方法已成为解决复杂问题的可靠且高效的工具。

Laguerre 方法：一种高效的多项式根近似技术

Laguerre 方法是一种强大的数值方法，用于近似多项式方程的根。它表现出卓越的有效性，尤其是在处理具有复数根的多项式时，通常优于包括著名的 Newton 方法在内的其他求根方法。该方法最初由 Laguerre 在 1898 年提出，是数值分析的重要贡献，如 McNamee [8] 和 [11] 的参考文献所示。

Laguerre 方法的一个显著特点是它同时使用了多项式函数 P(x) 的一阶和二阶导数。这种导数信息的整合导致了三阶收敛，这与第三部分中描述的 Halley 方法的收敛速率一致。多项式效率指数，表示为 3^(1/3)=1.44，由 Laguerre 方法和 Halley 方法共享。Laguerre 方法是

选择符号 ± 以最大化分母的绝对值，并且

其中 P(x) 是待求根的多项式，P’(x) 是一阶导数，P’’(x) 是多项式的二阶导数。我们注意到 P(x_n) 位于 G 和 H 的分母中，并且永不为零，否则我们将根据此标准停止搜索。

该方法具有对大多数函数进行全局收敛的优点，但由于需要计算一阶和二阶导数，因此计算成本比其他方法更高。

Laguerre 方法对于简单根具有三阶收敛速率。这意味着每次迭代的准确数字大约会增加两倍，使其成为寻找多项式根的非常高效的方法。然而，对于重根，收敛速率可能不同。三阶收敛使得 Laguerre 方法比许多其他求根技术（如二阶收敛速率的 Newton 方法）更快。

这是一个例子：

第一个根是一个复共轭根 x=(0.80-i1.22)。第三个根是一个简单的实根，约为 1.5。其余两个根直接从降阶的多项式中找到。

上图显示了接近根时的收敛速率，符合预期。

对于我们到目前为止介绍的所有求根方法，当存在重根时（重数 > 1），所有方法都会简化为线性收敛。同样，对于 Newton 和 Halley 方法，存在 Laguerre 方法的修改版本，即使对于重根也能保持三阶收敛速率。

修改后的版本是

其中 m 是根的重数。

大多数情况下，您事先不知道 m，但可以使用与第一部分至第三部分中介绍的相同技术（请参阅 Newton 方法的详细描述），其中我们继续使用 m=2 到 n 的以下公式，只要对于每个 m 都是。

修改后的 Laguerre 方法对实系数或复系数多项式都效果很好，是寻找多项式零点的一种非常稳定的方法。

比较 Laguerre 和 Newton 方法？

为了将不同的方法与其他方法进行比较，我们可以使用一个著名的效率指标来查看它与其他基于导数的方法相比如何。

效率指数为，其中 q 是方法的收敛阶数，p 是方法的多项式求值次数。对于 Newton 方法，p 为 2，因为我们需要每次迭代求值 P(z) 和 P'(z)，而 Newton 方法的收敛阶数为 q=2，因此效率指数为 2^(1/2)=1.42。

对于 Laguerre 方法，我们每次迭代需要求值 P(x)、P'(x) 和 P''(x)，因此效率指数为 3^(1/3)=1.44。

略大于 Newton 方法，但不足以让我们期望使用 Laguerre 方法能有任何可衡量的收益。

修改什么？

与 Newton 方法（第二部分）相比，我们可以很幸运地重用 Newton 方法中已有的绝大多数代码。

第二部分中的步骤包括：

寻找一个初始点
执行 Laguerre 迭代，包括通过 Horner 方法进行多项式求值
计算最终上限
多项式降阶
求解二次方程

第 1, 3, 4, 5 项）与 Newton 方法相同，无需修改。

第 2 项）我们可以不变地使用 Horner 方法求值 P(x)、P'(x) 和 P''(x)。但是，我们需要添加另一个向量来存储二阶导数 P''(x)。处理重数的变步长可以从m更改为前面列出的公式。否则，我们可以再次重用变步长，或者减小步长并在第一和第二部分中都显示它。

Laguerre 方法的实现

这是与第二部分和第三部分相同的源，除了需要评估 Laguerre 迭代步而不是 Newton 或 Halley 步的更改。

此求根器的实现遵循第一部分中概述的方法。

首先，我们消除零根（等于零的根）
然后我们找到一个合适的起始点开始 Laguerre 迭代，这还包括基于 P(x) 可接受值的终止条件，我们将在此停止当前迭代。
开始 Laguerre 迭代
1. 第一步是找到 dx_n=，当然，还要决定分母为零时应该发生什么。分母中的符号选择是为了最大化分母的绝对值。当这种情况出现时，通常是由于局部最小值，最佳做法是改变方向，乘以一个因子 dx_n=dx_n(0.6+i0.8)k。这相当于将方向旋转 53 度，并将方向乘以因子 k。k=5 是一个合理的值。
2. 此外，很容易意识到，如果 P'(x_n)~0，您可能会得到一个不合理的步长 dx_n，因此引入一个限制因子，如果 abs(dx_n)>5·abs(dx_n-1)（即当前步长大于前一次迭代步长的 5 倍），则减小当前步长。同样，如果发生这种情况，您将通过 dx_n=dx_n(0.6+i0.8)5(abs(dx_n-1)/abs(dx_n)) 来改变方向。
3. 这两个修改（a 和 b）使他的方法非常有弹性，并确保它总是收敛到根。
4. 下一个问题是如何处理重数 > 1 的问题，这会将三阶收敛速率降低到线性收敛速率。在找到合适的 dx_n 并得到新的 x_n+1=x_n-dx_n 后，我们检查 P(x_n+1)>P(x_n)：如果是，我们考虑修正后的 x_n+1=x_n-0.5dx_n，如果 P(x_n+1)≥P(x_n)，则我们使用原始的 x_n+1 作为下一次迭代的新起点。如果不是，那么我们接受 x_n+1 作为更好的选择，并继续考虑一个新修正的 x_n+1=x_n-0.25dx_n。另一方面，如果新的 P(x_n+1)≥P(x_n)，我们使用前一个 x_n+1 作为下一次迭代的新起点。如果不是，那么我们假设我们正在接近一个新的鞍点，并且方向通过 dx_n=dx_n(0.6+i0.8) 进行改变，我们使用 x_n+1=x_n-dx_n 作为下一次迭代的新起点。另一方面，如果 P(x_n+1)≤P(x_n)：那么我们正在朝着正确的方向前进，然后我们继续沿该方向前进，使用 , m=2,..,n，只要 P(x_n+1)≤P(x_n)，并为下一次迭代选择最佳的 m。这个过程的好处是，如果存在重数为 m 的根，那么 m 也是步长大小的最佳选择，这将保持三阶收敛速率，即使对于重根也是如此。
  进程 a-d 继续进行，直到达到停止条件，此时根 x_n 被接受并降阶到多项式中。然后对降阶后的多项式启动新的根搜索。

我们将迭代分为两个阶段。阶段 1 和阶段 2。在阶段 1 中，我们试图进入收敛圆，在那里我们确信 Laguerre 方法将收敛到一个根。由于这是不同的方法，我们不能使用与 Newton 方法相同的计算。但是，我们仍然会这样做，因为在实践中，它效果很好。当我们进入该圆后，我们自动切换到阶段 2。在阶段 2 中，我们跳过步骤 d)，只使用未修改的 Laguerre 步：，直到满足停止条件。万一我们移出了收敛圆，我们就切换回阶段 1 并继续迭代。

我们使用与 Newton 和 Halley 方法相同的标准来切换到阶段 2。

C++ 代码

下面的 C++ 代码使用 Laguerre 方法查找具有实系数多项式的多项式根。为了实现 Laguerre 方法，对 Newton 版本仅做了很少的更改。

/*
 *******************************************************************************
 *
 *                       Copyright (c) 2023
 *                       Henrik Vestermark
 *                       Denmark, USA
 *
 *                       All Rights Reserved
 *
 *   This source file is subject to the terms and conditions of
 *   Henrik Vestermark Software License Agreement which restricts the manner
 *   in which it may be used.
 *
 *******************************************************************************
*/
/*
 *******************************************************************************
 *
 * Module name     :   Laguerre.cpp
 * Module ID Nbr   :
 * Description     :   Solve n degree polynomial using Laguerre's method
 * --------------------------------------------------------------------------
 * Change Record   :
 *
 * Version Author/Date Description of changes
 * -------  ------------- ----------------------
 * 01.01 HVE/24Sep2023 Initial release
 *
 * End of Change Record
 * --------------------------------------------------------------------------
*/

// define version string 
static char _VLaguerre_[] = "@(#)testLaguerre.cpp 01.01 -- Copyright (C) Henrik Vestermark";

#include <algorithm>
#include <vector>
#include <complex>
#include <iostream>
#include <functional>

using namespace std;

constexpr int       MAX_ITER = 50;

// Find all polynomial zeros using a modified Laguerre method
// 1) Eliminate all simple roots (roots equal to zero)
// 2) Find a suitable starting point
// 3) Find a root using the Laguerre method
// 4) Divide the root up in the polynomial reducing its order with one
// 5) Repeat steps 2 to 4 until the polynomial is of the order of two 
//    whereafter the remaining one or two roots are found by the direct formula
// Notice:
//      The coefficients for p(x) is stored in descending order. coefficients[0] is a(n)x^n, 
//      coefficients[1] is a(n-1)x^(n-1),...,  coefficients[n-1] is a(1)x, coefficients[n] is a(0)
//
static vector<complex<double>> PolynomialRootsLaguerre(const vector<double>& coefficients)
{
    struct eval { complex<double> z{}; complex<double> pz{}; double apz{}; };
    const complex<double> complexzero(0.0);  // Complex zero (0+i0)
    size_t n;       // Size of Polynomial p(x)  
    eval pz;        // P(z)
    eval pzprev;    // P(zprev)
    eval p1z;       // P'(z)
    eval p1zprev;   // P'(zprev)
    complex<double> z;      // Use as temporary variable
    complex<double> dz;     // The current stepsize dz
    complex<double> newtondz;
    int itercnt;    // Hold the number of iterations per root
    vector<complex<double>> roots;             // Holds the roots of the Polynomial
    vector<double> coeff(coefficients.size()); // Holds the current coefficients of P(z)

    copy(coefficients.begin(), coefficients.end(), coeff.begin());
    // Step 1 eliminate all simple roots
    for (n = coeff.size() - 1; n > 0 && coeff.back() == 0.0; --n)
        roots.push_back(complexzero);  // Store zero as the root

    // Compute the next starting point based on the polynomial coefficients
    // A root will always be outside the circle from the origin and radius min
    auto startpoint = [&](const vector<double>& a)
    {
        const size_t n = a.size() - 1;
        double a0 = log(abs(a.back()));
        double min = exp((a0 - log(abs(a.front()))) / static_cast<double>(n));

        for (size_t i = 1; i < n; i++)
            if (a[i] != 0.0)
            {
                double tmp = exp((a0 - log(abs(a[i]))) / static_cast<double>(n - i));
                if (tmp < min)
                    min = tmp;
            }

        return min * 0.5;
    };

    // Evaluate a polynomial with real coefficients a[] at a complex point z and
    // return the result 
    // This is Horner's method, avoiding complex arithmetic
    auto horner = [](const vector<double>& a, const complex<double> z)
    {
        const size_t n = a.size() - 1;
        double p = -2.0 * z.real();
        double q = norm(z);
        double s = 0.0;
        double r = a[0];
        eval e;

        for (size_t i = 1; i < n; i++)
        {
            double t = a[i] - p * r - q * s;
            s = r;
            r = t;
        }

        e.z = z;
        e.pz = complex<double>(a[n] + z.real() * r - q * s, z.imag() * r);
        e.apz = abs(e.pz);
        return e;
    };

    // Calculate an upper bound for the rounding errors performed in a
    // polynomial with real coefficient a[] at a complex point z. 
    // (Adam's test)

    auto upperbound = [](const vector<double>& a, const complex<double> z)
    {
        const size_t n = a.size() - 1;
        double p = -2.0 * z.real();
        double q = norm(z);
        double u = sqrt(q);
        double s = 0.0;
        double r = a[0];
        double e = fabs(r) * (3.5 / 4.5);
        double t;

        for (size_t i = 1; i < n; i++)
        {
            t = a[i] - p * r - q * s;
            s = r;
            r = t;
            e = u * e + fabs(t);
        }
        t = a[n] + z.real() * r - q * s;
        e = u * e + fabs(t);
        e = (4.5 * e - 3.5 * (fabs(t) + fabs(r) * u) +
            fabs(z.real()) * fabs(r)) * 0.5 * pow((double)_DBL_RADIX, -DBL_MANT_DIG + 1);

        return e;
    };

    // Do Laguerre iteration for polynomial order higher than 2
    for (; n > 2; --n)
    {
        const double Max_stepsize = 5.0; // Allow the next step size to be 
                                         // up to 5 times larger than the 
                                         // previous step size
        const complex<double> rotation = complex<double>(0.6, 0.8);  // Rotation amount
        double r;              // Current radius
        double rprev;          // Previous radius
        double eps;            // The iteration termination value
        bool stage1 = true;    // By default it start the iteration in stage1
        int steps = 1;         // Multisteps if > 1
        eval p2z;              // P''(z)
        vector<double> coeffprime;   // vector holding the prime coefficients
        vector<double> coeffprime2;  // Laguerre vector holding both the 
                                     // prime and double prime coefficients

        // Calculate coefficients of p'(x)
        for (int i = 0; i < n; i++)
            coeffprime.push_back(coeff[i] * double(n - i));

        // Calculate coefficients of p''(x)
        for (int i = 0; i < n - 1; i++)       // Laguerre
            coeffprime2.push_back(coeffprime[i] * double(n - i - 1));   // Laguerre 

        // Step 2 find a suitable starting point z
        rprev = startpoint(coeff);      // Computed startpoint
        z = coeff[n - 1] == 0.0 ? complex<double>(1.0) : 
                                  complex<double>(-coeff[n] / coeff[n - 1]);
        z *= complex<double>(rprev) / abs(z);

        // Setup the iteration
        // Current P(z)
        pz = horner(coeff, z);

        // pzprev which is the previous z or P(0)
        pzprev.z = complex<double>(0);
        pzprev.pz = coeff[n];
        pzprev.apz = abs(pzprev.pz);

        // p1zprev P'(0) is the P'(0)
        p1zprev.z = pzprev.z;
        p1zprev.pz = coeff[n - 1];       // P'(0)
        p1zprev.apz = abs(p1zprev.pz);

        // Set previous dz and calculate the radius of operations.
        dz = pz.z;                 // dz=z-zprev=z since zprev==0
        r = rprev *= Max_stepsize; // Make a reasonable radius of the 
                                   // maximum step size allowed
        // Preliminary eps computed at point P(0) by a crude estimation
        eps = 2 * n * pzprev.apz * pow((double)_DBL_RADIX, -DBL_MANT_DIG);

        // Start iteration and stop if z doesnt change or apz <= eps
        // we do z+dz!=z instead of dz!=0. if dz does not change z 
        // then we accept z as a root
        for (itercnt = 0; pz.z + dz != pz.z && pz.apz > 
             eps && itercnt < MAX_ITER; itercnt++)
        {
            complex<double> G, H, gp, gm, u;

            // Calculate current P'(z) and P''(z)
            p1z = horner(coeffprime, pz.z);
            p2z = horner(coeffprime2, pz.z);
            // Compute G and H
            G = p1z.pz / pz.pz;
            H = G * G - p2z.pz / pz.pz;
            H = (complex<double>(static_cast<int>(n)) * H - G * G); // Save H for 
                                                                    // later=nH-G^2
            u = sqrt(complex<double>(static_cast<int>(n) - 1) * H);
            gp = G + u;
            gm = G - u;
            if (abs(gp) < abs(gm))
                gp = gm;
            // Calculate dz, change directions if zero
            if (abs(gp) == 0.0)                 // If Laguerre denominator is zero, 
                                                // then rotate previous dz direction
                dz *= rotation * complex<double>(Max_stepsize);
            else
                dz = complex<double>(static_cast<int>(n)) / gp;

            // Check the Magnitude of Laguerre's step
             r = abs(dz);
             if (r > rprev) // Large than 5 times the previous step size
             {   // then rotate and adjust step size to prevent 
                 // wild step size near P'(z) close to zero
                 dz *= rotation * complex<double>(rprev / r);
                 r = abs(dz);
             }
             rprev = r * Max_stepsize;  // Save 5 times the current step size 
                                        // for the next iteration check of 
                                        // reasonable step size

             // Calculate if stage1 is true or false. Stage1 is false 
             // if the Newton/Laguerre converge otherwise true
             z = (p1zprev.pz - p1z.pz) / (pzprev.z - pz.z);
             stage1 = (abs(z) / p1z.apz > p1z.apz / pz.apz / 4) || (steps != 1);
      
            // Step accepted. Save pz in pzprev

            pzprev = pz;
            z = pzprev.z - dz;      // Next z
            pz = horner(coeff, z);  //ff = pz.apz;
            steps = 1;

            if (stage1)
            {  // Try multiple steps or shorten steps depending if P(z) 
               // is an improvement or not P(z)<P(zprev)
                bool div2;
                complex<double> zn, dzn=dz;
                eval npz;

                steps++;
                for (div2 = pz.apz > pzprev.apz; steps <= n; ++steps)
                {
                    if (div2 == true)
                    {  // Shorten steps
                        dzn *= complex<double>(0.5);
                        zn = pzprev.z - dz;
                    }
                    else
                    {   // Compute new dz
                        dzn = sqrt(complex<double>(double(n)/double(steps) - 1) * H);
                        gp = G + dzn;
                        gm = G - dzn;
                        if (abs(gp) < abs(gm))
                            gp = gm;
                        dzn = complex<double>(static_cast<int>(n)) / gp; 
                        zn = pzprev.pz - dzn;
                    }
                    // Evaluate new try step
                    npz = horner(coeff, zn);
                    if (npz.apz >= pz.apz)
                    {
                        --steps; break; // Break if no improvement
                    }

                    // Improved => accept step and try another round of step
                    pz = npz;
                    dz = dzn;

                    if (div2 == true && steps == 3)
                    {   // To many shorten steps => try another direction and break
                        dz *= rotation;
                        z = pzprev.z - dz;
                        pz = horner(coeff, z);
                        break;
                    }
                }
            }
            else
            {   // calculate the upper bound of error using Grant & 
                // Hitchins's test for complex coefficients
                // Now that we are within the convergence circle.
                eps = upperbound(coeff, pz.z);
            }
        }

        // Real root forward deflation.
        //
        auto realdeflation = [&](vector<double>& a, const double x)
        {
            const size_t n = a.size() - 1;
            double r = 0.0;

            for (size_t i = 0; i < n; i++)
                a[i] = r = r * x + a[i];

            a.resize(n);    // Remove the highest degree coefficients.
        };

        // Complex root forward deflation for real coefficients
        //
        auto complexdeflation = [&](vector<double>& a, const complex<double> z)
        {
            const size_t n = a.size() - 1;
            double r = -2.0 * z.real();
            double u = norm(z);

            a[1] -= r * a[0];
            for (int i = 2; i < n - 1; i++)
                a[i] = a[i] - r * a[i - 1] - u * a[i - 2];

            a.resize(n - 1); // Remove top 2 highest degree coefficients
        };

        // Check if there is a very small residue in the imaginary part by trying
        // to evaluate P(z.real). if that is less than P(z). 
        // We take that z.real() is a better root than z.
        z = complex<double>(pz.z.real(), 0.0);
        pzprev = horner(coeff, z);
        if (pzprev.apz <= pz.apz)
        { // real root 
            pz = pzprev;
            // Save the root
            roots.push_back(pz.z);
            realdeflation(coeff, pz.z.real());
        }
        else
        {   // Complex root
            // Save the root
            roots.push_back(pz.z);
            roots.push_back(conj(pz.z));
            complexdeflation(coeff, pz.z);
            --n;
        }
    }   // End Iteration

    // Solve any remaining linear or quadratic polynomial
    // For Polynomial with real coefficients a[], 
    // The complex solutions are stored in the back of the roots
    auto quadratic = [&](const std::vector<double>& a)
    {
        const size_t n = a.size() - 1;
        complex<double> v;
        double r;

        // Notice that a[0] is !=0 since roots=zero has already been handle
        if (n == 1)
            roots.push_back(complex<double>(-a[1] / a[0], 0));
        else
        {
            if (a[1] == 0.0)
            {
                r = -a[2] / a[0];
                if (r < 0)
                {
                    r = sqrt(-r);
                    v = complex<double>(0, r);
                    roots.push_back(v);
                    roots.push_back(conj(v));
                }
                else
                {
                    r = sqrt(r);
                    roots.push_back(complex<double>(r));
                    roots.push_back(complex<double>(-r));
                }
            }
            else
            {
                r = 1.0 - 4.0 * a[0] * a[2] / (a[1] * a[1]);
                if (r < 0)
                {
                    v = complex<double>(-a[1] / (2.0 * a[0]), 
                                a[1] * sqrt(-r) / (2.0 * a[0]));
                    roots.push_back(v);
                    roots.push_back(conj(v));
                }
                else
                {
                    v = complex<double>((-1.0 - sqrt(r)) * a[1] / (2.0 * a[0]));
                    roots.push_back(v);
                    v = complex<double>(a[2] / (a[0] * v.real()));
                    roots.push_back(v);
                }
            }
        }
        return;
    };

    if (n > 0)
        quadratic(coeff);

    return roots;
}

示例 1

以下是上述源代码工作方式的示例。

对于实数多项式

+1x^4-10x^3+35x^2-50x+24

开始 Laguerre 迭代，多项式为 +1x^4-10x^3+35x^2-50x+24
阶段 1 => 停止条件。|f(z)|<2.13e-14
开始: z[1]=0.2 dz=2.40e-1 |f(z)|=1.4e+1

迭代：1

Laguerre 步: z[1]=1 dz=-7.46e-1 |f(z)|=8.9e-2
函数值减小 => 尝试在该方向进行多次步进
尝试步进: z[1]=1 dz=-7.46e-1 |f(z)|=8.1e-1
: 没有改进。放弃上一步

迭代：2

Laguerre 步: z[1]=1 dz=-1.45e-2 |f(z)|=2.1e-6
进入阶段 2。新的停止条件: |f(z)|<2.18e-14

迭代：3

Laguerre 步: z[1]=1 dz=-3.53e-7 |f(z)|=1.1e-14
进入阶段 2。新的停止条件: |f(z)|<2.18e-14
经过 3 次迭代满足停止条件
最终 Laguerre z[1]=1 dz=-3.53e-7 |f(z)|=1.1e-14
调整=0% 阶段 1=33% 阶段 2=67%
降阶实根 z=1.0000000000000009
开始 Laguerre 迭代，多项式为 +1x^3-9x^2+25.999999999999993x-23.999999999999986
阶段 1 => 停止条件。|f(z)|<1.60e-14
开始 : z[1]=0.5 dz=4.62e-1 |f(z)|=1.4e+1

迭代：1

Laguerre 步: z[1]=2 dz=-1.52e+0 |f(z)|=4.7e-2
函数值减小 => 尝试在该方向进行多次步进
尝试步进: z[1]=2 dz=-1.52e+0 |f(z)|=3.4e-1
: 没有改进。放弃上一步

迭代：2

Laguerre 步: z[1]=2 dz=-2.27e-2 |f(z)|=1.4e-6
进入阶段 2。新的停止条件: |f(z)|<1.42e-14

迭代：3

Laguerre 步: z[1]=2 dz=-6.91e-7 |f(z)|=0
进入阶段 2。新的停止条件: |f(z)|<1.42e-14
经过 3 次迭代满足停止条件
最终 Laguerre z[1]=2 dz=-6.91e-7 |f(z)|=0
调整=0% 阶段 1=33% 阶段 2=67%
降阶实根 z=2
直接求解多项式 +1x^2-7x+11.999999999999993

使用 Laguerre 方法，解为

X1=1.0000000000000009
X2=2
X3=4.000000000000007
X4=2.999999999999993

示例 2

同一个例子，只是有一个双根 x=1。我们看到每一步都是双步，与第一个根的重数为 2 相符。

对于实数多项式

+1x^4-9x^3+27x^2-31x+12

开始 Laguerre 迭代，多项式为 +1x^4-9x^3+27x^2-31x+12
阶段 1 => 停止条件。|f(z)|<1.07e-14
开始 : z[1]=0.2 dz=1.94e-1 |f(z)|=6.9e+0

迭代：1

Laguerre 步: z[1]=0.8 dz=-6.32e-1 |f(z)|=2.1e-1
函数值减小 => 尝试在该方向进行多次步进
尝试步进: z[1]=1 dz=-6.32e-1 |f(z)|=3.5e-6
: 改进。继续步进
尝试步进: z[1]=1 dz=-6.32e-1 |f(z)|=1.2e-1
: 没有改进。放弃上一步

迭代：2

Laguerre 步: z[1]=1 dz=-5.59e-4 |f(z)|=2.5e-7
函数值减小 => 尝试在该方向进行多次步进
尝试步进: z[1]=1 dz=-5.59e-4 |f(z)|=0
: 改进。继续步进
尝试步进: z[1]=1 dz=-5.59e-4 |f(z)|=2.5e-7
: 没有改进。放弃上一步
经过 2 次迭代满足停止条件
最终 Laguerre z[1]=1 dz=-5.59e-4 |f(z)|=0
调整=0% 阶段 1=100% 阶段 2=0%
降阶实根 z=0.9999999999981045
开始 Laguerre 迭代，多项式为 +1x^3-8.000000000001895x^2+19.00000000001327x-12.000000000022744
阶段 1 => 停止条件。|f(z)|<7.99e-15
开始 : z[1]=0.3 dz=3.16e-1 |f(z)|=6.8e+0

迭代：1

Laguerre 步: z[1]=1 dz=-6.83e-1 |f(z)|=6.3e-3
函数值减小 => 尝试在该方向进行多次步进
尝试步进: z[1]=1 dz=-6.83e-1 |f(z)|=1.2e+0
: 没有改进。放弃上一步

迭代：2

Laguerre 步: z[1]=1 dz=-1.05e-3 |f(z)|=4.8e-11
进入阶段 2。新的停止条件: |f(z)|<6.66e-15

迭代：3

Laguerre 步: z[1]=1 dz=-7.96e-12 |f(z)|=8.9e-16
进入阶段 2。新的停止条件: |f(z)|<6.66e-15
经过 3 次迭代满足停止条件
最终 Laguerre z[1]=1 dz=-7.96e-12 |f(z)|=8.9e-16
调整=0% 阶段 1=33% 阶段 2=67%
降阶实根 z=1.0000000000018952
直接求解多项式 +1x^2-7x+12.000000000000004

使用 Laguerre 方法，解为

X1=0.9999999999981045
X2=1.0000000000018952
X3=3.9999999999999964
X4=3.0000000000000036

示例 3

一个同时具有实根和复共轭根的测试多项式。

对于实数多项式

+1x^4-8x^3-17x^2-26x-40

开始 Laguerre 迭代，多项式为 +1x^4-8x^3-17x^2-26x-40
阶段 1 => 停止条件。|f(z)|<3.55e-14
开始 : z[1]=-0.8 dz=-7.67e-1 |f(z)|=2.6e+1

迭代：1

Laguerre 步: z[1]=-2 dz=1.17e+0 |f(z)|=1.9e+1
函数值减小 => 尝试在该方向进行多次步进
尝试步进: z[1]=-3 dz=1.17e+0 |f(z)|=9.3e+1
: 没有改进。放弃上一步

迭代：2

Laguerre 步: z[1]=-2 dz=-2.91e-1 |f(z)|=8.4e-2
进入阶段 2。新的停止条件: |f(z)|<1.91e-14

迭代：3

Laguerre 步: z[1]=-2 dz=1.58e-3 |f(z)|=2.0e-8
进入阶段 2。新的停止条件: |f(z)|<1.91e-14

迭代：4

Laguerre 步: z[1]=-2 dz=-3.81e-10 |f(z)|=2.8e-14
进入阶段 2。新的停止条件: |f(z)|<1.91e-14

迭代：5

Laguerre 步: z[1]=-2 dz=-5.34e-16 |f(z)|=0
进入阶段 2。新的停止条件: |f(z)|<1.91e-14
经过 5 次迭代满足停止条件
最终 Laguerre z[1]=-2 dz=-5.34e-16 |f(z)|=0
调整=0% 阶段 1=20% 阶段 2=80%
降阶实根 z=-1.650629191439388
开始 Laguerre 迭代，多项式为 +1x^3-9.650629191439387x^2-1.0703897408530487x-24.233183447530717
阶段 1 => 停止条件。|f(z)|<1.61e-14
开始 : z[1]=-0.8 dz=-7.92e-1 |f(z)|=3.0e+1

迭代：1

Laguerre 步: z[1]=(-0.4+i1) dz=(-4.08e-1-i1.45e+0) |f(z)|=7.2e+0
函数值减小 => 尝试在该方向进行多次步进
尝试步进: z[1]=(0.5+i2) dz=(-4.08e-1-i1.45e+0) |f(z)|=4.3e+1
: 没有改进。放弃上一步

迭代：2

Laguerre 步: z[2]=(-0.17+i1.5) dz=(-2.09e-1-i9.45e-2) |f(z)|=1.0e-2
进入阶段 2。新的停止条件: |f(z)|<1.36e-14

迭代：3

Laguerre 步: z[5]=(-0.17469+i1.5469) dz=(-5.79e-5+i3.16e-4) |f(z)|=2.7e-11
进入阶段 2。新的停止条件: |f(z)|<1.36e-14

迭代：4

Laguerre 步: z[13]=(-0.1746854042803+i1.546868887231) dz=(-7.87e-13+i3.54e-13) |f(z)|=3.6e-15
进入阶段 2。新的停止条件: |f(z)|<1.36e-14
经过 4 次迭代满足停止条件
最终 Laguerre z[13]=(-0.1746854042803+i1.546868887231) dz=(-7.87e-13+i3.54e-13) |f(z)|=3.6e-15
调整=0% 阶段 1=25% 阶段 2=75%
降阶复共轭根 z=(-0.17468540428030604+i1.5468688872313963)
直接求解多项式 +1x-10

使用 Laguerre 方法，解为

X1=-1.650629191439388
X2=(-0.17468540428030604+i1.5468688872313963)
X3=(-0.17468540428030604-i1.5468688872313963)
X4=10

迭代路径趋向于前两个根。

建议

由于效率指数与 Laguerre 方法相当，因此使用 Laguerre 方法而不是 Halley 或 Newton 方法没有优缺点。我建议为了简单起见，继续使用第一部分和第二部分中介绍的 Newton 方法。但是，如果您选择 Laguerre 方法，您不会失望的。

结论

我们介绍了一种改进的 Laguerre 方法，其收敛速率为 3（与 Halley 方法相当），并在第一、二、三部分建立的框架基础上，高效稳定地求解实系数多项式的根。我认为，是品味和偏好问题，决定是否使用 Laguerre 方法而不是例如 Newton 或 Halley 方法。一个展示这些各种方法的基于 Web 的多项式求解器可供进一步探索，您可以在多项式根查找器中找到它，该求解器演示了许多这些方法的实际应用。在第五部分，我们将考察一种所谓的同步方法（Durand-Kerner，又名 Weierstrass 方法）。

参考文献

H. Vestermark. A practical implementation of Polynomial root finders. Practical implementation of Polynomial root finders vs 7.docx (hvks.com)
Madsen. A root-finding algorithm based on Newton Method, Bit 13 (1973) 71-75。
A. Ostrowski, Solution of equations and systems of equations, Academic Press, 1966。
Wikipedia Horner’s Method: https://en.wikipedia.org/wiki/Horner%27s_method
Adams, D A stopping criterion for polynomial root finding. Communication of the ACM Volume 10/Number 10/ October 1967 Page 655-658
Grant, J. A. & Hitchins, G D. Two algorithms for the solution of polynomial equations to limiting machine precision. The Computer Journal Volume 18 Number 3, pages 258-264
Wilkinson, J H, Algebraic Processes中的舍入误差, Prentice-Hall Inc, Englewood Cliffs, NJ 1963
McNamee, J.M., Numerical Methods for Roots of Polynomials, Part I & II, Elsevier, Kidlington, Oxford 2009
H. Vestermark, “A Modified Newton and higher orders Iteration for multiple roots.”, www.hvks.com/Numerical/papers.html
M.A. Jenkins & J.F. Traub, ”A three-stage Algorithm for Real Polynomials using Quadratic iteration”, SIAM J Numerical Analysis, Vol. 7, No.4, December 1970。
Wikipedia Laguerre's method, https://en.wikipedia.org/wiki/Laguerre%27s_method

历史

2023年11月29日：原始提交中缺少两张图片

快速稳定的多项式求根器 - 第四部分

目录

引言

为什么选择 Laguerre 方法？

Laguerre 方法：一种高效的多项式根近似技术

比较 Laguerre 和 Newton 方法？

修改什么？

Laguerre 方法的实现

C++ 代码

示例 1

迭代：1

迭代：2

迭代：3

迭代：1

迭代：2

迭代：3

示例 2

迭代：1

迭代：2

迭代：1

迭代：2

迭代：3

示例 3

迭代：1

迭代：2

迭代：3

迭代：4

迭代：5

迭代：1

迭代：2

迭代：3

迭代：4

建议

结论

参考文献

历史