Skip to content

Newton's Method

Introduction

This article introduces Newton's method for finding roots, a technique proposed by Newton in the 17th century for approximating solutions to equations.

Specifically, the task is: for a continuous and monotonic function \(f(x)\) on \([a,b]\), find an approximate solution to the equation \(f(x)=0\).

Explanation

Initially, we start with a given \(f(x)\) and an approximate solution \(x_0\) (the choice of initial value relates to Newton fractals; refer to 3Blue1Brown's Newton Fractal).

Suppose our current approximate solution is \(x_i\). We draw the line \(l\) tangent to \(f(x)\) at the point \((x_i,f(x_i))\), and denote the x-coordinate of its intersection with the x-axis as \(x_{i+1}\). This gives a better approximate solution. Repeat this iterative process. Based on the geometric meaning of the derivative, we have:

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

Rearranging yields the recurrence relation:

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

Intuitively, if \(f(x)\) is sufficiently smooth, then as the number of iterations increases, \(x_i\) will converge to the true solution.

Newton's method has quadratic convergence, meaning the number of accurate digits approximately doubles with each iteration. For a proof of convergence of Newton's method, refer to citizendium - Newton method Convergence analysis

Of course, Newton's method also has limitations; see Xiaolin Wu - Roots of Equations pages 18 - 20 analysis for details.

Computing Square Roots

We attempt to use Newton's method to compute square roots. Let \(f(x)=x^2-n\), then the approximate solution to this equation is an approximation of \(\sqrt{n}\). Thus we obtain:

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

When implementing, pay attention to setting an appropriate precision. Code is as follows:

Implementation

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
double sqrt_newton(double n) {
  constexpr static 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

Computing Integer Square Root

Although we can call the sqrt() function to obtain the square root value, here we discuss a variant of Newton's method for finding the maximum integer solution to the inequality \(x^2\le n\). We still consider a process similar to Newton iteration but with slight modifications to the boundary conditions. If during iteration the approximate solution decreases from the previous step and then increases in the current step, we skip this iteration and exit the loop.

Implementation

 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

High-Precision Square Root

Finally, consider Newton's method for high-precision computation. The iterative method remains unchanged, but we need to pay attention to the initial approximate solution, i.e., the value of \(x_0\). Since numbers requiring high precision are typically very large, the choice of initial value significantly affects algorithm efficiency. A natural idea is to set \(x_0=2^{\left\lfloor\frac{1}{2}\log_2n\right\rfloor}\), which allows quick computation of \(x_0\) while being relatively close to the square root approximation.

Implementation

Java code implementation:

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

Practical performance: For \(n=10^{1000}\), the runtime of this algorithm is 60 ms. If we do not optimize the initial value \(x_0\) and start iterating directly from \(x_0=1\), the runtime increases to 120 ms.

Exercises

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用