利用牛顿迭代法求平方根

来源：互联网发布：人工智能第一股编辑：程序博客网时间：2024/04/28 16:35

求n的平方根，先假设一猜测值X₀ = 1，然后根据以下公式求出X₁，再将X₁代入公式右边，继续求出X₂…通过有效次迭代后即可求出n的平方根，X_k+1

先让我们来验证下这个巧妙的方法准确性，来算下2的平方根 (Computed by Mathomatic)

1-> x_new = ( x_old + y/x_old )/2y(x_old + -----)x_old#1: x_new = ---------------21-> calculate x_old 1Enter y: 2Enter initial x_old: 1 x_new = 1.51-> calculate x_old 2Enter y: 2Enter initial x_old: 1 x_new = 1.41666666666671-> calculate x_old 3Enter y: 2Enter initial x_old: 1 x_new = 1.41421568627451-> calculate x_old 10Enter y: 2Enter initial x_old: 1Convergence reached after 6 iterations. x_new = 1.4142135623731...

可见，随着迭代次数的增加，运算值会愈发接近真实值。很神奇的算法，可是怎么来的呢? 查了下wikipedia和wolfram，原来算法的名字叫Newton’s Iteration (牛顿迭代法)。

下面是极其つまらない(boring)的数理介绍，不喜欢数学的言下之意也就是绝大部分人可以略过了。

简单推导

假设f(x)是关于X的函数:

An illustration of one iteration of Newton's method

求出f(x)的一阶导，即斜率:

$f'(x_{n}) = /frac{ /mathrm{rise} }{ /mathrm{run} } = /frac{ /mathrm{/Delta y} }{ /mathrm{/Delta x} } = /frac{ f( x_{n} ) - 0 }{ x_{n} - x_{n+1} } = /frac{0 - f(x_{n})}{(x_{n+1} - x_{n})}/,/!$

简化等式得到:

然后利用得到的最终式进行迭代运算直至求到一个比较精确的满意值，为什么可以用迭代法呢?理由是中值定理(Intermediate Value Theorem):

如果f函数在闭区间[a,b]内连续，必存在一点x使得f(x) = c，c是函数f在闭区间[a,b]内的一点

我们先猜测一X初始值，例如1，当然地球人都知道除了1本身之外任何数的平方根都不会是1。然后代入初始值，通过迭代运算不断推进，逐步靠近精确值，直到得到我们主观认为比较满意的值为止。例如要求768的平方根，因为25² = 625，而30² = 900，我们可先代入一猜测值26，然后迭代运算，得到较精确值:27.7128。

回到我们最开始的那个”莫名其妙”的公式，我们要求的是N的平方根，令x² = n，假设一关于X的函数f(x)为:

f(X) = X² - n

求f(X)的一阶导为:

f'(X) = 2X

代入前面求到的最终式中:

X_k+1 = X_k - (X_k² - n)/2X_k

化简即得到我们最初提到的那个求平方根的神奇公式了:

用泰勒公式推导

我之前介绍过在The Art and Science of C一书中有用到泰勒公式求平方根的算法，其实牛顿迭代法也可以看作是泰勒公式(Taylor Series)的简化，先回顾下泰勒公式:

仅保留等式右边前两项:

令f(X₀+ε) = 0，得到:

再令X₁ = X₀ + ε₀，得到ε₁…依此类推可知:

转化为:

引申

从推导来看，其实牛顿迭代法不仅可以用来求平方根，还可以求立方根，甚至更复杂的运算。

同样，我们还可以利用C语言来实现下那个最简单的求平方根的公式(尽管我们可以直接用sqrt()完成)

#include <stdio.h>#include <math.h>#define N 768main() {float x=1;int i;for (i=1;i<=1000;i++) {  // recursion times : 1000x = (x + N/x)/2;}printf("The square root of %d is %f/n",N,x);}