牛顿迭代法

引入

本文介绍如何用牛顿迭代法(Newton's method for finding roots)求方程的近似解,该方法于 17 世纪由牛顿提出。

具体的任务是,对于在 上连续且单调的函数 ,求方程 的近似解。

解释

初始时我们从给定的 和一个近似解 开始(初值的问题与 Newton 分形有关,可参考 3Blue1Brown 的 牛顿分形)。

假设我们目前的近似解是 ,我们画出与 切于点 的直线 ,将 轴的交点横坐标记为 ,那么这就是一个更优的近似解。重复这个迭代的过程。 根据导数的几何意义,可以得到如下关系:

整理后得到如下递推式:

直观地说,如果 比较平滑,那么随着迭代次数的增加, 会越来越逼近方程的解。

牛顿迭代法的收敛率是平方级别的,这意味着每次迭代后近似解的精确数位会翻倍。 关于牛顿迭代法的收敛性证明可参考 citizendium - Newton method Convergence analysis

当然牛顿迭代法也同样存在着缺陷,详情参考 Xiaolin Wu - Roots of Equations 第 18 - 20 页分析

求解平方根

我们尝试用牛顿迭代法求解平方根。设 ,这个方程的近似解就是 的近似值。于是我们得到

在实现的时候注意设置合适的精度。代码如下

实现

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
// C++ Version
double sqrt_newton(double n) {
  const double eps = 1E-15;
  double x = 1;
  while (true) {
    double nx = (x + n / x) / 2;
    if (abs(x - nx) < eps) break;
    x = nx;
  }
  return x;
}
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# Python Version
def sqrt_newton(n):
    eps = 1e-15
    x = 1
    while True:
        nx = (x + n / x) / 2
        if abs(x - nx) < eps:
            break
        x = nx
    return x

求解整数平方根

尽管我们可以调用 sqrt() 函数来获取平方根的值,但这里还是讲一下牛顿迭代法的变种算法,用于求不等式 的最大整数解。我们仍然考虑一个类似于牛顿迭代的过程,但需要在边界条件上稍作修改。如果 在迭代的过程中上一次迭代值得近似解变小,而这一次迭代使得近似解变大,那么我们就不进行这次迭代,退出循环。

实现

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
// C++ Version
int isqrt_newton(int n) {
  int x = 1;
  bool decreased = false;
  for (;;) {
    int nx = (x + n / x) >> 1;
    if (x == nx || (nx > x && decreased)) break;
    decreased = nx < x;
    x = nx;
  }
  return x;
}
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
# Python Version
def isqrt_newton(n):
    x = 1
    decreased = False
    while True:
        nx = (x + n // x) // 2
        if x == nx or (nx > x and decreased):
            break
        decreased = nx < x
        x = nx
    return x

高精度平方根

最后考虑高精度的牛顿迭代法。迭代的方法是不变的,但这次我们需要关注初始时近似解的设置,即 的值。由于需要应用高精度的数一般都非常大,因此不同的初始值对于算法效率的影响也很大。一个自然的想法就是考虑 ,这样既可以快速计算出 ,又可以较为接近平方根的近似解。

实现

给出 Java 代码的实现:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
public static BigInteger isqrtNewton(BigInteger n) {
  BigInteger a = BigInteger.ONE.shiftLeft(n.bitLength() / 2);
  boolean p_dec = false;
  for (;;) {
    BigInteger b = n.divide(a).add(a).shiftRight(1);
    if (a.compareTo(b) == 0 || a.compareTo(b) < 0 && p_dec)
      break;
    p_dec = a.compareTo(b) > 0;
    a = b;
  }
  return a;
}

实践效果:在 的时候该算法的运行时间是 60 ms,如果我们不优化 的值,直接从 开始迭代,那么运行时间将增加到 120 ms。

习题


评论