昨天查mathomatic的MAN时看到一个很有意思的计算平方根的方法(当然不是叫你在Excel里n^1/2啦!):

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

x_(k+1)=1/2(x_k+n/(x_k))

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

1-> x_new = ( x_old + y/x_old )/2

                       y
            (x_old + -----)
                     x_old
#1: x_new = ---------------
                   2

1-> calculate x_old 1
Enter y: 2
Enter initial x_old: 1

 x_new = 1.5

1-> calculate x_old 2
Enter y: 2
Enter initial x_old: 1

 x_new = 1.4166666666667

1-> calculate x_old 3
Enter y: 2
Enter initial x_old: 1

 x_new = 1.4142156862745

1-> calculate x_old 10
Enter y: 2
Enter initial x_old: 1
Convergence 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})}\,\!

简化等式得到:

x_(n+1)=x_n-(f(x_n))/(f^'(x_n))

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

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

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

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

f(X) = X2 - n

求f(X)的一阶导为:

f'(X) = 2X

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

Xk+1 = Xk - (Xk2 - n)/2Xk

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

x_(k+1)=1/2(x_k+n/(x_k))

用泰勒公式推导

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

f(x_0+epsilon)=f(x_0)+f^'(x_0)epsilon+1/2f^('')(x_0)epsilon^2+....

仅保留等式右边前两项:

f(x_0+epsilon) approx f(x_0)+f^'(x_0)epsilon.

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

epsilon_0=-(f(x_0))/(f^'(x_0))

再令X1 = X0 + ε0，得到ε1…依此类推可知:

epsilon_n=-(f(x_n))/(f^'(x_n))

转化为:

x_(n+1)=x_n-(f(x_n))/(f^'(x_n))

引申

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

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

#include <stdio.h>
#include <math.h>

#define N 768

main() {
        float x=1;
        int i;
        for (i=1;i<=1000;i++) {  // recursion times : 1000
                x = (x + N/x)/2;
        }
        printf("The square root of %d is %f\n",N,x);
}

14 Reponses So Far ↓

  1. hermy:

    丫我一周上了8节数学课我orz你还给我看数学你有没有良心啊!!!

  2. punkid:

    这不正是教你怎么用你那8节课所学的东西吗! 不然就课堂上老师教你的那点，除了能考试做题，还能干嘛。

    你看，我这多好的一免费老师啊!

  3. Tony:

    我对代数没啥兴趣

  4. cyndi:

    我是不是总留下与文章主题无关的评论啊。
    这一次我想说，我想要个好看的模板呀。
    。。。

  5. punkid:

    上面留言的都是O.T者

    第二句话暂时忽略…做个blogbus的theme还好说，wordpress的工程量太大了，而且我向来都是心血来潮的话主动给别人做theme，而不是别人要我做theme :D

  6. fcicq:

    求平方根立方根都不应该需要除法… 自己继续找资料去吧 :D

  7. punkid:

    给个思路啊? 啥算法来着?

  8. feng:

    mathomatic??????

  9. punkid:

    Mathomatic is a portable, general purpose CAS (Computer Algebra System) written entirely in the C programming language. It is free and open source software, published under the GNU Lesser General Public License (LGPL version 2.1). This is a console mode application that does symbolic math and quick calculations. It compiles and runs under any operating system with a C compiler, especially GNU/Linux.

  11. punkid:

    对于算法我可是个门外汉啊…写这种一知半解的东西纯粹给自己提个醒的 :D
    正在研究那篇文章里的数理

  12. Shan:

    “再令X1 = X0 + ε0，得到ε1…依此类推可知:”
    请问εn 表示什么？
    结尾不明白诶，能不能在解释一下？

  14. Shan:

    哦，得到εn的步骤其实就是迭代的过程~~赞~
    我是看过The Art…再来看的，谢拉~

Leave a Reply ↓