快速稳定的多项式求根器 - 第二部分
4.60/5 (3投票s)
具有实数系数的多项式的快速稳定求根器
目录
- 引言
- 简易方法
- 实数或复数根的陷阱
- 修改什么？
- 牛顿法的问题
- 求根器的合适起点
- 复数点处的多项式求值
- 合适的根停止准则
- 算术运算中的误差
- 多项式降阶策略
- K. Madsen 牛顿算法的实现
- C++ 代码
- 结论
- 参考文献
- 历史
引言
在第一篇论文（第一部分）中，我们开发了一种基于牛顿法的高效且鲁棒的多项式求根器，专门针对复数多项式系数。尽管许多工程问题涉及实数系数多项式，但核心方法仍然适用。将系数限制为实数确实增加了实现的复杂性，但也由于与复数算术相比，实数算术的计算开销减少而加快了计算速度。
简易方法
现在系数是实数，我们不再将多项式系数存储在 std::vector<std::complex<double>> 中，而是使用 std::vector<double>。一个直接的修改是将每个系数引用包装成 complex<double>(coefficients[i])。虽然这种方法可行，但它忽略了复数根总是成对出现这一特性，导致不必要的转换并降低了效率。
实数或复数根的陷阱
在我看来，为实数系数实现一个求根器比为复数系数设计的求根器更具挑战性。这是因为每次迭代可能产生一个实数根或一对复数根。复数根总是作为共轭对出现，需要格外小心以避免将实数根误分类为两个复数根。如果算法进入复平面但返回一个接近实轴的值，则需要进行测试以确认根的性质。找到一个根（指定为 z）后，我们检查 P(z.real()) <= P(z)。如果为 true，我们接受 z.real() 作为实数根并相应地对多项式进行降阶。否则，我们识别 z 及其复共轭为根，并使用二次公式对多项式进行降阶，将其次数减二。除此之外，第一篇论文中解决的多重根问题保持不变。
修改什么？
系数现在存储在 vector<double> 中，而不是 vector<complex<double>> 中。需要注意的是，复数根仍然以共轭对形式出现。
根据第一部分，步骤包括
- 找到一个初始起点
- 执行牛顿迭代，包括通过霍纳方法进行多项式求值
- 计算 P(z) 求值误差的最终上限
- 多项式降阶
- 求解二次方程
Ad 1) 初始点寻找算法保持不变，无需代码修改。
Ad 2) 尽管算法保持不变，但优化实数系数的霍纳方法会减少算术运算。需要进行代码调整以实现更高效的实现。
Ad 3) 实数系数简化了上限误差计算。我们从 Grant-Hutchins 算法切换到 Adam 算法以提高效率。
Ad 4) 与复数系数一次找到一个根不同，我们现在可能找到一个实数根或一对复数根。因此，多项式可以被一个实数根或两个复数根的二次形式除。
Ad 5) 求解二次或线性多项式的算法不同，因为它必须考虑根是实数还是复数。这在仅处理复数根时不是问题。
审查每个修改版本
- 霍纳方法
- 上限误差计算
- 实数或共轭复数降阶
- 线性或二次解
牛顿法的问题
在第一部分中，我们详细阐述了牛顿法的问题，并解释了多重根问题。为了克服这个问题，我们需要将牛顿步长乘以一个因子 m，从而得到修正的牛顿法。
其中 m 是根的重数，这都是众所周知的事情。挑战在于如何在实际生活中找到 m。
求根器的合适起点
我们使用与第一部分中相同的算法来寻找合适的起点，并且我们只展示了实数系数多项式的修改函数。
该算法为我们的根搜索计算了一个合理且合适的起点。
// 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;
};
复数点处的多项式求值
大多数求根方法都要求我们在某个点评估多项式。
为了评估多项式 P(z)，其中 P(z)=anzn+an-1zn-1,…,a1z+a0
我们使用霍纳 [4] 方法，由以下递推式给出
bn=an
bk=bk-1z+ak, k=n-1,…,0
P(z)=b0
然而，为了避免复数算术开销，您可以使用霍纳的优化方法，仅使用实数算术。
在具有实数系数的多项式 P(z) 在复数点 z 处求值的情况下，我们通常使用霍纳递推，但采用仅使用实数算术的特殊版本
z=x+iy
p=-2x
q=x2+y2
bn=an
bn-1=an-1-pbn
bk=ak-pbk+1-qbk+2, k=n-2,…,1
b0=a0-xb1-qb2
P(z)=b0+iyb1
此递推式的最后一项是 P(z) 的值。霍纳方法长期以来被认为是求多项式在给定点处最有效的方法。此算法适用于实数或复数点处的实数系数。
// Evaluate a polynomial with real coefficients a[] at a complex point z and
// return the result
// This is Horner's method of 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;
};
合适的根停止准则
在 [8] 中，他们介绍了许多不同的技术来计算合适的停止准则。另请参阅 [1]。许多求根器可以使用 Adams [5] 或 Hitching [6] 使用的方法来为具有实数或复数系数的多项式找到合适的停止准则。
在进行迭代方法时，您将需要在某个时候考虑要为您的求根器应用什么停止准则。由于大多数迭代求根器都使用多项式的求值来推进，因此自然而然地继续我们的搜索，直到 P(z) 的求值足够接近 0 才能在该点接受根。
算术运算中的误差
J.H.Wilkinson 在 [7] 中指出，执行代数运算的误差受以下限制
ε<½β1-t，其中 β 是基数，t 是精度（假设四舍五入到最接近）
注意 ½β1-t= β-t
对于 Intel 微处理器系列和 IEE754 浮点运算标准，64 位浮点算术的 β= 2 和 t=53，即 2-53
一个简单的上限
然后可以使用上述信息为 n 次多项式找到一个简单的上限。
| 多项式 | ||
| 操作数 | 实数系数 | 复数系数 |
| 实数点 | |ao|·2n·2-53 | |ao|·4n·2-53 |
| 复数点 | |ao|·4n·2-53 | |ao|·6n·2-53 |
一个更好的上限
在此类别中，我们有 Adams [5] 和 Grant & Hitchins [6] 等多项式停止准则。
多项式求根器通常可以处理实数或复数系数的多项式在实数或复数处求值。由于 Grant-Hitchin 的停止准则适用于复数系数多项式，我们将使用 Adams 的类似边界，但适用于实数系数多项式。
使用前面章节所示的霍纳方法，对在复数点 z 处求值的实数系数 an 的多项式使用算法。
Adams [5] 已经表明，可以使用以下递推式计算误差边界
en=|bn|7/9
ek=ek-1|z|+|bk|, k=n-1,...,0
e=(4.5e0-3.5(|b0|+|b1||z|)+|x||b1|)½β1-t
还存在其他有用的方法，请参阅 [1]。
// 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;
};
多项式降阶策略
找到一个根后，我们需要对当前多项式进行该根的综合除法，以降低多项式次数并准备寻找下一个根。问题随之而来，您是使用正向降阶还是反向降阶？
Wilkinson [7] 已经表明，要实现稳定的降阶过程，如果您按递增的量级查找多项式的根，则应选择正向降阶，并且始终首先使用量级最小的根来降阶多项式，当然，当按递减的量级查找根时，则选择相反的反向降阶。
为了进行正向降阶，我们尝试从最高系数 an 开始求解方程
anzn+an-1zn-1+⋯+a1z+a0=(bn-1zn-1+bn-2zn-2+⋯+b1z+b0)(z-R)
R 是根。
现在求解 b，您得到递推式
bn-1=an
bk=ak+1+R∙bk+1, k=n-2,…,0
这个简单的算法对于实数系数和实数根的多项式效果很好。一个特殊情况是实数系数和复数根。一个复数根及其复共轭根将与将多项式 P(Z) 除以两个复共轭根 (x+iy) 和 (x-iy) 的二次多项式（或 (z2-2xz+(x2+y2))）相同。令 r=-2x 且 u= x2+y2
你得到
P(z)=Q(z)(z2+rz+u)
其中 P(z)=anzn+an-1zn-1+⋯+a1z+a0
Q(z)=bn-2zn-2+bn-3zn-3+⋯+b1z+b0
递推式然后由下式给出
an=bn-2
an-1=bn-3+rbn-2
ak-2=bk-2+rbk-1+ubk, k=n-2,...,3
a2=b0+rb1+ub2
a1=rb0+ub1
a0=ub0
现在求解 b，您得到
bn-2=an
bn-3=an-1-rbn-2
bk=ak+2-rbk+1-ubk, k=n-4,…,0
// 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
};
K. Madsen 牛顿算法的实现
我们遵循与第一部分中相同的牛顿法实现，仅在必要时进行修改以适应适当的实数算术，并调用修改后的函数以在复数点评估多项式时查找误差上限，以及实数或复数根的多项式降阶。
C++ 代码
下面的 C++ 代码找到实数系数多项式的多项式根。复数系数多项式的相同版本可在第一部分中找到。详见 [1]。
/*
*******************************************************************************
*
* 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 : Newton.cpp
* Module ID Nbr :
* Description : Solve n degree polynomial using Newton's (Madsen) method
* --------------------------------------------------------------------------
* Change Record :
*
* Version Author/Date Description of changes
* ------- ------------- ----------------------
* 01.01 HVE/24Sep2023 Initial release
*
* End of Change Record
* --------------------------------------------------------------------------
*/
// define version string
static char _VNEWTON_[] = "@(#)Newton.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 Newton method
// 1) Eliminate all simple roots (roots equal to zero)
// 2) Find a suitable starting point
// 3) Find a root using Newton
// 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>> PolynomialRoots(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
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 of 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 a 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 Newton 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
vector<double> coeffprime;
// Calculate coefficients of p'(x)
for (int i = 0; i < n; i++)
coeffprime.push_back(coeff[i] * double(n - i));
// 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++)
{
// Calculate current P'(z)
p1z = horner(coeffprime, pz.z);
if (p1z.apz == 0.0) // P'(z)==0 then rotate and try again
dz *= rotation * complex<double>(Max_stepsize); // Multiply with 5
// to get away from saddlepoint
else
{
dz = pz.pz / p1z.pz; // next dz
// Check the Magnitude of Newton'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 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;
eval npz;
zn = pz.z;
for (div2 = pz.apz > pzprev.apz; steps <= n; ++steps)
{
if (div2 == true)
{ // Shorten steps
dz *= complex<double>(0.5);
zn = pzprev.z - dz;
}
else
zn -= dz; // try another step in the same direction
// Evaluate new try step
npz = horner(coeff, zn);
if (npz.apz >= pz.apz)
break; // Break if no improvement
// Improved => accept step and try another round of step
pz = npz;
if (div2 == true && steps == 2)
{ // 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 coefficienst
};
// 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
以下是上述源代码如何工作的一个示例。请注意详细说明阶段 1 和阶段 2 花费时间以及需要多少次旋转的统计数据。
For the real Polynomial:
+1x^4-10x^3+35x^2-50x+24
Start Newton Iteration for Polynomial=+1x^4-10x^3+35x^2-50x+24
Stage 1=>Stop Condition. |f(z)|<2.13e-14
Start : z[1]=0.2 dz=2.40e-1 |f(z)|=1.4e+1
Iteration: 1
Newton Step: z[1]=0.6 dz=-3.98e-1 |f(z)|=3.9e+0
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=1 dz=-3.98e-1 |f(z)|=2.0e-1
: Improved=>Continue stepping
Try Step: z[1]=1 dz=-3.98e-1 |f(z)|=9.9e-1
: No improvement=>Discard last try step
Iteration: 2
Newton Step: z[1]=1 dz=3.87e-2 |f(z)|=1.6e-2
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=1 dz=3.87e-2 |f(z)|=2.7e-1
: No improvement=>Discard last try step
Iteration: 3
Newton Step: z[1]=1 dz=-2.62e-3 |f(z)|=7.6e-5
In Stage 2=>New Stop Condition: |f(z)|<2.18e-14
Iteration: 4
Newton Step: z[1]=1 dz=-1.26e-5 |f(z)|=1.8e-9
In Stage 2=>New Stop Condition: |f(z)|<2.18e-14
Iteration: 5
Newton Step: z[1]=1 dz=-2.93e-10 |f(z)|=3.6e-15
In Stage 2=>New Stop Condition: |f(z)|<2.18e-14
Stop Criteria satisfied after 5 Iterations
Final Newton z[1]=1 dz=-2.93e-10 |f(z)|=3.6e-15
Alteration=0% Stage 1=40% Stage 2=60%
Deflate the real root z=1.0000000000000007
Start Newton Iteration for Polynomial=+1x^3-9x^2+25.999999999999993x-23.99999999999999
Stage 1=>Stop Condition. |f(z)|<1.60e-14
Start : z[1]=0.5 dz=4.62e-1 |f(z)|=1.4e+1
Iteration: 1
Newton Step: z[1]=1 dz=-7.54e-1 |f(z)|=3.9e+0
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=2 dz=-7.54e-1 |f(z)|=6.4e-2
: Improved=>Continue stepping
Try Step: z[1]=3 dz=-7.54e-1 |f(z)|=2.6e-1
: No improvement=>Discard last try step
Iteration: 2
Newton Step: z[1]=2 dz=-2.95e-2 |f(z)|=2.7e-3
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=2 dz=-2.95e-2 |f(z)|=5.4e-2
: No improvement=>Discard last try step
Iteration: 3
Newton Step: z[1]=2 dz=-1.32e-3 |f(z)|=5.3e-6
In Stage 2=>New Stop Condition: |f(z)|<1.42e-14
Iteration: 4
Newton Step: z[1]=2 dz=-2.63e-6 |f(z)|=2.1e-11
In Stage 2=>New Stop Condition: |f(z)|<1.42e-14
Iteration: 5
Newton Step: z[1]=2 dz=-1.04e-11 |f(z)|=3.6e-15
In Stage 2=>New Stop Condition: |f(z)|<1.42e-14
Stop Criteria satisfied after 5 Iterations
Final Newton z[1]=2 dz=-1.04e-11 |f(z)|=3.6e-15
Alteration=0% Stage 1=40% Stage 2=60%
Deflate the real root z=2.0000000000000036
Solve Polynomial=+1x^2-6.9999999999999964x+11.999999999999975 directly
Using the Newton Method, the Solutions are:
X1=1.0000000000000007
X2=2.0000000000000036
X3=4.000000000000011
X4=2.999999999999986
示例 2
相同的示例，只是在 x=1 处有一个重根。我们看到每一步都是双倍步长，与第一个根的重数 2 一致。请注意，对于多重根，它将所有搜索时间都花在了 stage1 中，而 stage2 中为 0%。
For the real Polynomial:
+1x^4-9x^3+27x^2-31x+12
Start Newton Iteration for Polynomial=+1x^4-9x^3+27x^2-31x+12
Stage 1=>Stop Condition. |f(z)|<1.07e-14
Start : z[1]=0.2 dz=1.94e-1 |f(z)|=6.9e+0
Iteration: 1
Newton Step: z[1]=0.5 dz=-3.23e-1 |f(z)|=2.0e+0
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=0.8 dz=-3.23e-1 |f(z)|=1.8e-1
: Improved=>Continue stepping
Try Step: z[1]=1 dz=-3.23e-1 |f(z)|=1.4e-1
: Improved=>Continue stepping
Try Step: z[1]=1 dz=-3.23e-1 |f(z)|=8.9e-1
: No improvement=>Discard last try step
Iteration: 2
Newton Step: z[1]=1 dz=8.71e-2 |f(z)|=3.1e-2
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=1 dz=8.71e-2 |f(z)|=9.6e-4
: Improved=>Continue stepping
Try Step: z[1]=0.9 dz=8.71e-2 |f(z)|=6.5e-2
: No improvement=>Discard last try step
Iteration: 3
Newton Step: z[1]=1 dz=-6.27e-3 |f(z)|=2.4e-4
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=1 dz=-6.27e-3 |f(z)|=2.6e-8
: Improved=>Continue stepping
Try Step: z[1]=1 dz=-6.27e-3 |f(z)|=2.3e-4
: No improvement=>Discard last try step
Iteration: 4
Newton Step: z[1]=1 dz=-3.28e-5 |f(z)|=6.5e-9
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=1 dz=-3.28e-5 |f(z)|=0
: Improved=>Continue stepping
Try Step: z[1]=1 dz=-3.28e-5 |f(z)|=6.5e-9
: No improvement=>Discard last try step
Stop Criteria satisfied after 4 Iterations
Final Newton z[1]=1 dz=-3.28e-5 |f(z)|=0
Alteration=0% Stage 1=100% Stage 2=0%
Deflate the real root z=0.9999999982094424
Start Newton Iteration for Polynomial=+1x^3-8.000000001790557x^2+19.000000012533903x-12.000000021486692
Stage 1=>Stop Condition. |f(z)|<7.99e-15
Start : z[1]=0.3 dz=3.16e-1 |f(z)|=6.8e+0
Iteration: 1
Newton Step: z[1]=0.8 dz=-4.75e-1 |f(z)|=1.5e+0
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=1 dz=-4.75e-1 |f(z)|=1.3e+0
: Improved=>Continue stepping
Try Step: z[1]=2 dz=-4.75e-1 |f(z)|=2.1e+0
: No improvement=>Discard last try step
Iteration: 2
Newton Step: z[1]=0.9 dz=3.54e-1 |f(z)|=5.7e-1
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=0.6 dz=3.54e-1 |f(z)|=3.7e+0
: No improvement=>Discard last try step
Iteration: 3
Newton Step: z[1]=1 dz=-8.28e-2 |f(z)|=3.6e-2
In Stage 2=>New Stop Condition: |f(z)|<6.64e-15
Iteration: 4
Newton Step: z[1]=1 dz=-5.87e-3 |f(z)|=1.7e-4
In Stage 2=>New Stop Condition: |f(z)|<6.66e-15
Iteration: 5
Newton Step: z[1]=1 dz=-2.88e-5 |f(z)|=4.1e-9
In Stage 2=>New Stop Condition: |f(z)|<6.66e-15
Iteration: 6
Newton Step: z[1]=1 dz=-6.90e-10 |f(z)|=8.9e-16
In Stage 2=>New Stop Condition: |f(z)|<6.66e-15
Stop Criteria satisfied after 6 Iterations
Final Newton z[1]=1 dz=-6.90e-10 |f(z)|=8.9e-16
Alteration=0% Stage 1=33% Stage 2=67%
Deflate the real root z=1.0000000017905577
Solve Polynomial=+1x^2-6.999999999999999x+12 directly
Using the Newton Method, the Solutions are:
X1=0.9999999982094424
X2=1.0000000017905577
X3=3.999999999999995
X4=3.0000000000000036
示例 3
一个同时具有实根和复共轭根的测试多项式。
For the real Polynomial:
+1x^4-8x^3-17x^2-26x-40
Start Newton Iteration for Polynomial=+1x^4-8x^3-17x^2-26x-40
Stage 1=>Stop Condition. |f(z)|<3.55e-14
Start : z[1]=-0.8 dz=-7.67e-1 |f(z)|=2.6e+1
Iteration: 1
Newton Step: z[1]=-2 dz=1.65e+0 |f(z)|=7.0e+1
Function value increase=>try shorten the step
Try Step: z[1]=-2 dz=8.24e-1 |f(z)|=3.1e+0
: Improved=>Continue stepping
Try Step: z[1]=-1 dz=4.12e-1 |f(z)|=1.8e+1
: No improvement=>Discard last try step
Iteration: 2
Newton Step: z[1]=-2 dz=6.27e-2 |f(z)|=1.5e-1
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=-2 dz=6.27e-2 |f(z)|=3.7e+0
: No improvement=>Discard last try step
Iteration: 3
Newton Step: z[1]=-2 dz=-2.74e-3 |f(z)|=2.9e-4
In Stage 2=>New Stop Condition: |f(z)|<1.91e-14
Iteration: 4
Newton Step: z[1]=-2 dz=-5.51e-6 |f(z)|=1.2e-9
In Stage 2=>New Stop Condition: |f(z)|<1.91e-14
Iteration: 5
Newton Step: z[1]=-2 dz=-2.22e-11 |f(z)|=0
In Stage 2=>New Stop Condition: |f(z)|<1.91e-14
Stop Criteria satisfied after 5 Iterations
Final Newton z[1]=-2 dz=-2.22e-11 |f(z)|=0
Alteration=0% Stage 1=40% Stage 2=60%
Deflate the real root z=-1.650629191439388
Start Newton Iteration for Polynomial=+1x^3-9.650629191439387x^2-1.0703897408530487x-24.233183447530717
Stage 1=>Stop Condition. |f(z)|<1.61e-14
Start : z[1]=-0.8 dz=-7.92e-1 |f(z)|=3.0e+1
Iteration: 1
Newton Step: z[1]=1 dz=-1.86e+0 |f(z)|=3.5e+1
Function value increase=>try shorten the step
Try Step: z[1]=0.1 dz=-9.30e-1 |f(z)|=2.5e+1
: Improved=>Continue stepping
Try Step: z[1]=-0.3 dz=-4.65e-1 |f(z)|=2.5e+1
: No improvement=>Discard last try step
Iteration: 2
Newton Step: z[1]=-7 dz=6.71e+0 |f(z)|=7.2e+2
Function value increase=>try shorten the step
Try Step: z[1]=-3 dz=3.35e+0 |f(z)|=1.5e+2
: Improved=>Continue stepping
Try Step: z[1]=-2 dz=1.68e+0 |f(z)|=4.9e+1
: Improved=>Continue stepping
: Probably local saddlepoint=>Alter Direction:
z[1]=(-0.4-i0.7) dz=(5.03e-1+i6.71e-1) |f(z)|=2.1e+1
Iteration: 3
Newton Step: z[1]=(0.3-i2) dz=(-6.86e-1+i1.17e+0) |f(z)|=1.9e+1
Function value decrease=>try multiple steps in that direction
Try Step: z[1]=(1-i3) dz=(-6.86e-1+i1.17e+0) |f(z)|=8.4e+1
: No improvement=>Discard last try step
Iteration: 4
Newton Step: z[1]=(-0.09-i1) dz=(4.11e-1-i3.44e-1) |f(z)|=3.0e+0
In Stage 2=>New Stop Condition: |f(z)|<1.29e-14
Iteration: 5
Newton Step: z[2]=(-0.18-i1.5) dz=(8.71e-2+i4.72e-2) |f(z)|=1.1e-1
In Stage 2=>New Stop Condition: |f(z)|<1.36e-14
Iteration: 6
Newton Step: z[3]=(-0.175-i1.55) dz=(-3.24e-3+i1.00e-3) |f(z)|=1.3e-4
In Stage 2=>New Stop Condition: |f(z)|<1.36e-14
Iteration: 7
Newton Step: z[6]=(-0.174685-i1.54687) dz=(1.28e-6+i3.82e-6) |f(z)|=1.8e-10
In Stage 2=>New Stop Condition: |f(z)|<1.36e-14
Iteration: 8
Newton Step: z[12]=(-0.174685404280-i1.54686888723) dz=(-2.07e-12-i5.33e-12) |f(z)|=3.6e-15
In Stage 2=>New Stop Condition: |f(z)|<1.36e-14
Stop Criteria satisfied after 8 Iterations
Final Newton z[12]=(-0.174685404280-i1.54686888723) dz=(-2.07e-12-i5.33e-12) |f(z)|=3.6e-15
Alteration=13% Stage 1=38% Stage 2=63%
Deflate the complex conjugated root z=(-0.17468540428030604-i1.5468688872313963)
Solve Polynomial=+1x-10 directly
Using the Newton Method, the Solutions are:
X1=-1.650629191439388
X2=(-0.17468540428030604-i1.5468688872313963)
X3=(-0.17468540428030604+i1.5468688872313963)
X4=10
结论
我们提出了一个改进的牛顿方法，基于第一部分建立的框架，以高效稳定地找到实数系数多项式的根。第一部分侧重于复数系数多项式——其中根仍然可以是实数——而第二部分深入探讨了实数系数多项式。第三部分将探讨高阶方法（如哈雷方法）所需的调整，而第四部分将展示将拉盖尔方法等替代方法轻松集成到现有框架中。一个展示这些各种方法的基于网络的多项式求解器可供进一步探索，可在 http://www.hvks.com/Numerical/websolver.php 上找到。多项式求根器演示了其中许多方法的实际应用。
参考文献
- H. Vestermark. A practical implementation of Polynomial root finders. Practical implementation of Polynomial root finders. www.hvks.com/Numerical/papers.html
- 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。
- 维基百科霍纳方法：https://en.wikipedia.org/wiki/Horner%27s_method
- Adams, D. 多项式求根的停止准则。
- Communications of the ACM 第 10 卷/第 10 期/ 1967 年 10 月 第 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，“多重根的修正牛顿法和高阶迭代。”，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。
历史
- 2023 年 10 月 15 日：第二部分首次发布