牛顿迭代法(Newton-Raphson)实现sqrt

Posted by qingdujun on 2019-07-09

牛顿法最初由艾萨克·牛顿在《流数法》(Method of Fluxions,1671年完成,在牛顿去世后的1736年公开发表)中提出。约瑟夫·鲍易也曾于1690年在Analysis Aequationum中提出此方法。

首先,选择一个接近函数 $f(x)$零点的$x{0}$,计算相应的$f(x_0)$和切线斜率 $f’(x_0)$(这里$f’$表示函数$f$的导数)。然后我们计算穿过点$(x{0},f(x_{0})$并且斜率为$f’(x_0)$的直线和$x$轴的交点的$x$坐标,也就是求如下方程的解:

  • 我们将新求得的点的$x$坐标命名为$x1$,通常$x_1$会比$x{0}$更接近方程$f(x)=0$的解。因此我们现在可以利用$x_1$开始下一轮迭代。迭代公式可化简为如下所示:

对于$sqrt$问题,

  • 则存在 $f(x) = x^2 - a$ 要求解$x$的值
  • 此时 $f’(x) = 2 * x$

以下为牛顿迭代法求开根号的C++版本实现… (实现了$f(x)$、$f’(x)$函数,将上述迭代公式带入即可)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
double origin(double x, int a) {
return pow(x, 2) - a; //f(x) = x^2 - a
}

double derivative(double x) {
return 2 * x; //f'(x) = 2 * x
}

double newton_raphson(double a, double precision = 1e-5, int max_iter_num = 50) {
if (fabs(a) < precision) {
return 0.0;
}

double x0 = a, x1 = 0.0;

for (int k = 0; k < max_iter_num; ++k) {
x1 = x0 - origin(x0, a) / derivative(x0); //newton-raphson
if (fabs(x1-x0) < precision) {
break;
}
x0 = x1;
}

return x1;
}

References:
[1] 维基百科,牛顿法, https://zh.wikipedia.org/wiki/%E7%89%9B%E9%A1%BF%E6%B3%95
[2] 果壳问答,求牛顿开方法的算法及其原理,此算法能开任意次方吗?