二分法和牛頓迭代法

先說一個面試題：問 1.2 - 0.2  ==  1 ？

　　答案是False！ 為什么？

其原因在於十進制和二進制的轉換上，計算機先要把十進制的數轉化為二進制，然后再計算。
但是，在轉化中，浮點數轉化為二進制，就出問題了，例如：
十進制的 0.1，轉化為二進制是：0.0001100110011001100110011001100110011001100110011…（不能精確）
也就是說，轉化為二進制后，不會精確等於十進制的 0.1。
同時，計算機存儲的位數是有限制的，所以，就出現上述現象了。

這種問題不僅僅是 Python 中有，所有支持浮點數運算的編程語言都會遇到，它不是 Python 的 bug 

 

 

 

 

求一個數的平方根函數sqrt(int num)，在大多數語言中都提供實現。那么要求一個數的平方根，是怎么實現的哪？

　　實際上求平方根的算法方法主要有兩種：二分法( binary search)和牛頓迭代法( Newton iteration):

說到這兩種方法，心里就感慨重重，記得大學開了一門叫《計算方法》的課程，后來就忘了，直到有一天老師講到

這里的時候才想起這種計算方法自己曾幾何時學過！

 

1、二分法：（夾逼法）　　

　　以求根號5為例：

　　a:折半：       5/2=2.5

　　b:平方校驗:  2.5*2.5=6.25>5，並且得到當前上限2.5

　　c:再次向下折半:2.5/2=1.25

　　d:平方校驗：1.25*1.25=1.5625<5,得到當前下限1.25

　　e:再次折半:2.5-(2.5-1.25)/2=1.875

　　f:平方校驗：1.875*1.875=3.515625<5,得到當前下限1.875

 

每次得到當前和5進行比較，並且記下下限和上線，依次迭代，逐漸逼近我們設好的誤差范圍：

import math
from math import sqrt
 
def sqrt_binary(num):
    x=sqrt(num)    # 求出系統給我們的准確值
    y=num/2.0　　　　#中分
    low=0.0
    up=num*1.0
    count=1
    while abs(y-x)>0.00000001:
        print count,y          
        count+=1           
        if (y*y>num):
            up=y           #二分后中間值比預期大了
            y=low+(y-low)/2
        else:
            low=y
            y=up-(up-y)/2   #二分后中間值比預期小了
    return y
 
print(sqrt_binary(5))
print(sqrt(5))

 

運行結果:

　　1 2.5
　　2 1.25
　　3 1.875
　　4 2.1875
　　5 2.34375
　　6 2.265625
　　7 2.2265625
　　8 2.24609375
　　9 2.236328125
　　10 2.2314453125
　　11 2.23388671875
　　12 2.23510742188
　　13 2.23571777344
　　14 2.23602294922
　　15 2.23617553711
　　16 2.23609924316
　　17 2.23606109619
　　18 2.23608016968
　　19 2.23607063293
　　20 2.23606586456
　　21 2.23606824875
　　22 2.23606705666
　　23 2.2360676527
　　24 2.23606795073
　　25 2.23606809974
　　26 2.23606802523
　　27 2.23606798798
2.23606796935
2.2360679775

經過27次二分迭代，得到的值和系統sqrt()差別在0.0000001，

 

2、牛頓迭代法：

仔細思考一下就能發現，我們需要解決的問題可以簡單化理解。

從函數意義上理解：我們是要求函數f(x) = x²，使f(x) = num的近似解，即x² - num = 0的近似解。

從幾何意義上理解：我們是要求拋物線g(x) = x² - num與x軸交點（g(x) = 0）最接近的點。

 

我們假設g(x0)=0，即x0是正解，那么我們要做的就是讓近似解x不斷逼近x0，這是函數導數的定義

 

 可以由此得到： (可以寫出f'(x)的導數形式）   

 

從幾何圖形上看，因為導數是切線，通過不斷迭代，導數與x軸的交點會不斷逼近x0。

所以對於一般情況：

將m=2帶入得：

代碼如下：

def sqrt_newton(num):
    x=sqrt(num)
    y=num/2.0
    count=1
    while abs(y-x)>0.00000001:
        print count,y
        count+=1
        y=((y*1.0)+(1.0*num)/y)/2.0000
    return y
 
print(sqrt_newton(5))
print(sqrt(5))

 

運算結果： 

 1 2.5
 2 2.25
 3 2.23611111111

 

對於牛頓迭代法來說：還可以求立方，甚至等多（只需要修改m得值即可）

def cube_newton(num):
    x=num/3.0
    y=0
    count=1
    while abs(x-y)>0.00000001:
        print count,x
        count+=1
        y=x
        x=(2.0/3.0)*x+(num*1.0)/(x*x*3.0)
    return x
 
print(cube_newton(27))