跳转至

牛頓迭代法

引入

本文介紹如何用牛頓迭代法（Newton's method for finding roots）求方程的近似解，該方法於 17 世紀由牛頓提出。

具體的任務是，對於在 \([a,b]\) 上連續且單調的函數 \(f(x)\)，求方程 \(f(x)=0\) 的近似解。

解釋

初始時我們從給定的 \(f(x)\) 和一個近似解 \(x_0\) 開始（初值的問題與 Newton 分形有關，可參考 3Blue1Brown 的 牛頓分形）。

假設我們目前的近似解是 \(x_i\)，我們畫出與 \(f(x)\) 切於點 \((x_i,f(x_i))\) 的直線 \(l\)，將 \(l\) 與 \(x\) 軸的交點橫座標記為 \(x_{i+1}\)，那麼這就是一個更優的近似解。重複這個迭代的過程。 根據導數的幾何意義，可以得到如下關係：

\[ f'(x_i) = \frac{f(x_i)}{x_{i} - x_{i+1}} \]

整理後得到如下遞推式：

\[ x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)} \]

直觀地説，如果 \(f(x)\) 比較平滑，那麼隨着迭代次數的增加，\(x_i\) 會越來越逼近方程的解。

牛頓迭代法的收斂率是平方級別的，這意味着每次迭代後近似解的精確數位會翻倍。 關於牛頓迭代法的收斂性證明可參考 citizendium - Newton method Convergence analysis

當然牛頓迭代法也同樣存在着缺陷，詳情參考 Xiaolin Wu - Roots of Equations 第 18 - 20 頁分析

求解平方根

我們嘗試用牛頓迭代法求解平方根。設 \(f(x)=x^2-n\)，這個方程的近似解就是 \(\sqrt{n}\) 的近似值。於是我們得到

\[ x_{i+1}=x_i-\frac{x_i^2-n}{2x_i}=\frac{x_i+\frac{n}{x_i}}{2} \]

在實現的時候注意設置合適的精度。代碼如下

實現

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
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
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() 函數來獲取平方根的值，但這裏還是講一下牛頓迭代法的變種算法，用於求不等式 \(x^2\le n\) 的最大整數解。我們仍然考慮一個類似於牛頓迭代的過程，但需要在邊界條件上稍作修改。如果 \(x\) 在迭代的過程中上一次迭代值得近似解變小，而這一次迭代使得近似解變大，那麼我們就不進行這次迭代，退出循環。

實現

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
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
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

高精度平方根

最後考慮高精度的牛頓迭代法。迭代的方法是不變的，但這次我們需要關注初始時近似解的設置，即 \(x_0\) 的值。由於需要應用高精度的數一般都非常大，因此不同的初始值對於算法效率的影響也很大。一個自然的想法就是考慮 \(x_0=2^{\left\lfloor\frac{1}{2}\log_2n\right\rfloor}\)，這樣既可以快速計算出 \(x_0\)，又可以較為接近平方根的近似解。

實現

給出 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;
}

實踐效果：在 \(n=10^{1000}\) 的時候該算法的運行時間是 60 ms，如果我們不優化 \(x_0\) 的值，直接從 \(x_0=1\) 開始迭代，那麼運行時間將增加到 120 ms。

習題

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：Marcythm, iamtwz, nutshellfool, sshwy, allenanswerzq, countercurrent-time, Enter-tainer, H-J-Granger, hly1204, Ir1d, Menci, NachtgeistW, SukkaW, Tiphereth-A, Xeonacid
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用