Brian Warehouse

Some birds aren`t meant to be caged, their feathers are just too bright... ...
posts - 40, comments - 16, trackbacks - 0, articles - 1

牛顿法求平方根【转】

Posted on 2010-08-17 13:49 Brian 阅读(710) 评论(0)  编辑 收藏 引用 所属分类: 概念和技术

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

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
...

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

下面是数理介绍,不喜欢数学的言下之意也就是绝大部分人可以略过了。

简单推导

假设DE>f(x)DE>是关于DE>XDE>的函数:

An illustration of on<wbr>e iteration of Newton's method

求出DE>f(x)DE>的一阶导,即斜率:

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):

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

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

回到我们最开始的那个”莫名其妙”的公式,我们要求的是DE>NDE>的平方根,令DE>x2 = nDE>,假设一关于DE>XDE>的函数DE>f(x)DE>为:

DE>f(X) = X2 - nDE>

DE>f(X)DE>的一阶导为:

DE>f'(X) = 2XDE>

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

DE>Xk+1 = Xk - (Xk2 - n)/2XkDE>

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

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.

DE>f(X0+ε) = 0DE>,得到:

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

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

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

转化为:

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


只有注册用户登录后才能发表评论。
网站导航: 博客园   IT新闻   BlogJava   博问   Chat2DB   管理