很牛的牛頓迭代法

在mit公開課《電腦科學與程式設計導論》的第五講中，講到編寫求解平方根的函式sqrt時，提到了牛頓迭代法。今天仔細一查，發現這是乙個用途很廣、很牛的計算方法。

首先，考慮如何編寫乙個開平方根的函式sqrt(float num, float e)。引數num是要求開平方根的實數，引數e是計算結果可以達到多大誤差。這是乙個無法得到精確解，只能求出近似解的問題。該如何編寫呢？

1. 傳統的二分法

我們可以先猜測乙個值作為解，看這個值是否在誤差範圍內。如果沒有達到誤差要求，就構造乙個更好的猜測值，繼續迭代。猜測值可以簡單地初始化為num/2，但怎樣在下一輪迭代前構造乙個更好的猜測值呢？我們不妨參照二分查詢演算法的思路，解的初始範圍是[0, num]，用二分法逐漸縮小範圍。

private

static

float

sqrt(

float

num,

float

e)else

count++;

system.

out.printf(

"try %f, e: %f\n"

, guess, e0); }

while

(e0 > e);

system.

out.printf(

"try %d times, result: %f\n"

, count, guess);

return

guess; }

在區間[low, high]內取中點(low+high)/2作為猜測值。如果guess*guess大於num，說明猜測值偏大，則在區間[low, guess]進行下一輪迭代，否則在區間[guess, high]中繼續。當誤差值e0小於能夠接受的誤差值e時停止迭代，返回結果。

取num=2, e=0.01進行測試，執行結果如下：

try 1.000000, e: 1.000000

try 1.500000, e: 0.250000

try 1.250000, e: 0.437500

try 1.375000, e: 0.109375

try 1.437500, e: 0.066406

try 1.406250, e: 0.022461

try 1.421875, e: 0.021729

try 1.414063, e: 0.000427

try 8 times, result: 1.414063

可見嘗試了八次才達到0.01的誤差。

2. 神奇的牛頓法

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

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

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

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

可以由此得到

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

3. 牛頓法的實現與測試

public

static

void

main(string args)

private

static

float

sqrtnewton(

float

num,

float

e)while

(e0 > e);

system.

out.printf(

"try %d times, result: %f\n"

, count, guess);

return

guess; }

與二分法的對比測試結果：

try 1.000000, e: 1.000000

try 1.500000, e: 0.250000

try 1.250000, e: 0.437500

try 1.375000, e: 0.109375

try 1.437500, e: 0.066406

try 1.406250, e: 0.022461

try 1.421875, e: 0.021729

try 1.414063, e: 0.000427

try 8 times, result: 1.414063

try 1.500000, e: 0.250000

try 1.416667, e: 0.006945

try 2 times, result: 1.416667

try 1.000000, e: 1.000000

try 1.500000, e: 0.250000

try 1.250000, e: 0.437500

try 1.375000, e: 0.109375

try 1.437500, e: 0.066406

try 1.406250, e: 0.022461

try 1.421875, e: 0.021729

try 1.414063, e: 0.000427

try 1.417969, e: 0.010635

try 1.416016, e: 0.005100

try 1.415039, e: 0.002336

try 1.414551, e: 0.000954

try 1.414307, e: 0.000263

try 1.414185, e: 0.000082

try 14 times, result: 1.414185

try 1.500000, e: 0.250000

try 1.416667, e: 0.006945

try 1.414216, e: 0.000006

try 3 times, result: 1.414216

try 1.000000, e: 1.000000

try 1.500000, e: 0.250000

try 1.250000, e: 0.437500

try 1.375000, e: 0.109375

try 1.437500, e: 0.066406

try 1.406250, e: 0.022461

try 1.421875, e: 0.021729

try 1.414063, e: 0.000427

try 1.417969, e: 0.010635

try 1.416016, e: 0.005100

try 1.415039, e: 0.002336

try 1.414551, e: 0.000954

try 1.414307, e: 0.000263

try 1.414185, e: 0.000082

try 1.414246, e: 0.000091

try 1.414215, e: 0.000004

try 16 times, result: 1.414215

try 1.500000, e: 0.250000

try 1.416667, e: 0.006945

try 1.414216, e: 0.000006

try 3 times, result: 1.414216

可以看到隨著對誤差要求的更加精確，二分法的效率很低下，而牛頓法的確非常高效。

可在兩三次內得到結果。

如果搞不清牛頓法的具體原理，可能就要像我一樣複習下數學知識了。極限、導數、泰勒展開式、單變數微分等。

4. 更快的方法

在quake原始碼中有段求sqrt的方法，大概思路是只進行一次牛頓迭代，得到能夠接受誤差範圍內的結果。

因此肯定是更快的。

float q_rsqrt( float number )

參考文章

quake3原始碼中的sqrt

牛頓迭代方程的近似解

很牛的牛頓迭代法

牛頓迭代法

牛頓迭代法

牛頓迭代法

相關推薦