當計算1.0E21和1.0E23的平方根時，爲什麼Newton-Raphson方法不會收斂？

Question

我有以下函數用於計算使用Newton-Raphson方法的數字的平方根。

double get_dx(double y, double x) 
{ 
    return (x*x - y)/(2*x); 
} 

double my_sqrt(double y, double initial) 
{ 
    double tolerance = 1.0E-6; 
    double x = initial; 
    double dx; 
    int iteration_count = 0; 
    while (iteration_count < 100 && 
      fabs(dx = get_dx(y, x)) > tolerance) 
    { 
     x -= dx; 
     ++iteration_count; 
     if (iteration_count > 90) 
     { 
     printf("Iteration number: %d, dx: %lf\n", iteration_count, dx); 
     } 
    } 

    if (iteration_count < 100) 
    { 
     printf("Got the result to converge in %d iterations.\n", iteration_count); 
    } 
    else 
    { 
     printf("Could not get the result to converge.\n"); 
    } 

    return x; 
} 

它們適用於大多數數字。但是，當y爲1.0E21和1.0E23時，它們不收斂。可能還有其他數字，我還不知道，但功能不會收斂。開始使用的東西 -

1. y*0.5：

   我的初始值進行測試。

2. 1.0E10 - 最終結果附近的數字。
3. sqrt(y)+1.0 - 顯然，與最終答案非常接近。

在所有這些情況下，1.0E21和1.0E23的函數均未收斂。我嘗試了一個較低的號碼，1.0E19和更高的號碼，1.0E25。它們都不是問題。

這裏有一個完整的程序：

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

double get_dx(double y, double x) 
{ 
    return (x*x - y)/(2*x); 
} 

double my_sqrt(double y, double initial) 
{ 
    double tolerance = 1.0E-6; 
    double x = initial; 
    double dx; 
    int iteration_count = 0; 
    while (iteration_count < 100 && 
      fabs(dx = get_dx(y, x)) > tolerance) 
    { 
     x -= dx; 
     ++iteration_count; 
     if (iteration_count > 90) 
     { 
     printf("Iteration number: %d, dx: %lf\n", iteration_count, dx); 
     } 
    } 

    if (iteration_count < 100) 
    { 
     printf("Got the result to converge in %d iterations.\n", iteration_count); 
    } 
    else 
    { 
     printf("Could not get the result to converge.\n"); 
    } 

    return x; 
} 

void test(double y) 
{ 
    double ans = my_sqrt(y, 0.5*y); 
    printf("sqrt of %lg: %lg\n\n", y, ans); 

    ans = my_sqrt(y, 1.0E10); 
    printf("sqrt of %lg: %lg\n\n", y, ans); 

    ans = my_sqrt(y, sqrt(y) + 1.0); 
    printf("sqrt of %lg: %lg\n\n", y, ans); 
} 

int main() 
{ 
    test(1.0E21); 
    test(1.0E23); 
    test(1.0E19); 
    test(1.0E25); 
} 

和這裏的輸出（建立在64位Linux，用gcc 4.8.4）：

Iteration number: 91, dx: -0.000002 
Iteration number: 92, dx: 0.000002 
Iteration number: 93, dx: -0.000002 
Iteration number: 94, dx: 0.000002 
Iteration number: 95, dx: -0.000002 
Iteration number: 96, dx: 0.000002 
Iteration number: 97, dx: -0.000002 
Iteration number: 98, dx: 0.000002 
Iteration number: 99, dx: -0.000002 
Iteration number: 100, dx: 0.000002 
Could not get the result to converge. 
sqrt of 1e+21: 3.16228e+10 

Iteration number: 91, dx: 0.000002 
Iteration number: 92, dx: -0.000002 
Iteration number: 93, dx: 0.000002 
Iteration number: 94, dx: -0.000002 
Iteration number: 95, dx: 0.000002 
Iteration number: 96, dx: -0.000002 
Iteration number: 97, dx: 0.000002 
Iteration number: 98, dx: -0.000002 
Iteration number: 99, dx: 0.000002 
Iteration number: 100, dx: -0.000002 
Could not get the result to converge. 
sqrt of 1e+21: 3.16228e+10 

Iteration number: 91, dx: 0.000002 
Iteration number: 92, dx: -0.000002 
Iteration number: 93, dx: 0.000002 
Iteration number: 94, dx: -0.000002 
Iteration number: 95, dx: 0.000002 
Iteration number: 96, dx: -0.000002 
Iteration number: 97, dx: 0.000002 
Iteration number: 98, dx: -0.000002 
Iteration number: 99, dx: 0.000002 
Iteration number: 100, dx: -0.000002 
Could not get the result to converge. 
sqrt of 1e+21: 3.16228e+10 

Iteration number: 91, dx: 0.000027 
Iteration number: 92, dx: 0.000027 
Iteration number: 93, dx: 0.000027 
Iteration number: 94, dx: 0.000027 
Iteration number: 95, dx: 0.000027 
Iteration number: 96, dx: 0.000027 
Iteration number: 97, dx: 0.000027 
Iteration number: 98, dx: 0.000027 
Iteration number: 99, dx: 0.000027 
Iteration number: 100, dx: 0.000027 
Could not get the result to converge. 
sqrt of 1e+23: 3.16228e+11 

Iteration number: 91, dx: -0.000027 
Iteration number: 92, dx: -0.000027 
Iteration number: 93, dx: -0.000027 
Iteration number: 94, dx: -0.000027 
Iteration number: 95, dx: -0.000027 
Iteration number: 96, dx: -0.000027 
Iteration number: 97, dx: -0.000027 
Iteration number: 98, dx: -0.000027 
Iteration number: 99, dx: -0.000027 
Iteration number: 100, dx: -0.000027 
Could not get the result to converge. 
sqrt of 1e+23: 3.16228e+11 

Iteration number: 91, dx: 0.000027 
Iteration number: 92, dx: 0.000027 
Iteration number: 93, dx: 0.000027 
Iteration number: 94, dx: 0.000027 
Iteration number: 95, dx: 0.000027 
Iteration number: 96, dx: 0.000027 
Iteration number: 97, dx: 0.000027 
Iteration number: 98, dx: 0.000027 
Iteration number: 99, dx: 0.000027 
Iteration number: 100, dx: 0.000027 
Could not get the result to converge. 
sqrt of 1e+23: 3.16228e+11 

Got the result to converge in 35 iterations. 
sqrt of 1e+19: 3.16228e+09 

Got the result to converge in 6 iterations. 
sqrt of 1e+19: 3.16228e+09 

Got the result to converge in 1 iterations. 
sqrt of 1e+19: 3.16228e+09 

Got the result to converge in 45 iterations. 
sqrt of 1e+25: 3.16228e+12 

Got the result to converge in 13 iterations. 
sqrt of 1e+25: 3.16228e+12 

Got the result to converge in 1 iterations. 
sqrt of 1e+25: 3.16228e+12 

任何人都可以解釋爲什麼上面的功能不要」 t收斂爲1.0E21和1.0E23？

Answer 1

該答案具體解釋了爲什麼算法無法收斂1e23的輸入。

您面臨的問題被稱爲「大量小差異」。具體而言，您計算的是(x*x - y)/(2*x)，其中y是1e23，而x大約是3.16e11。

1e23的IEEE-754 encoding是0x44b52d02c7e14af6。因此指數是0x44b，它是1099的十進制數。由於指數編碼爲偏移二進制，因此必須減少1023。然後你必須減去另一個52才能找到LSB的權重。

0x44b ==> 1099 - 1023 - 52 = 24 ==> 1 LSB = 2^24 

所以一個LSB​​的值是16777216

結果從這個代碼

double y = 1e23; 
double x = sqrt(y); 
double dx = x*x - y; 
printf("%lf\n", dx); 

其實-16777216。

其結果是，當你計算x*x - y，結果是要麼將是零，或者它會通過這是2*3.16e11是的16777216劃分16777216多給0.000027的錯誤，這是不是在您公差。

兩個最接近sqrt(y)值

0x4252682931c035a0 
0x4252682931c035a1 

，當你方這些數字你

0x44b52d02c7e14af5 
0x44b52d02c7e14af7 

所以沒有一個符合預期的結果，這是

0x44b52d02c7e14af6 

，因此算法永遠不會收斂。

---

該算法適用於1e25的原因是由於運氣。爲1e25編碼是

0x45208b2a2c280291 

爲sqrt(1e25)編碼是

0x428702337e304309 

，當你方這個數字，你會得到

0x45208b2a2c280291 

因此x*x - y正好是0，和算法幸運。

Answer 2

這是一個浮點精度問題，連同使用固定的tolerance。 （1）一個「假」的零差異，這不是一個問題，或（2）一個錯誤的零差異，這是不成問題的，或者（2）由於差異持續大於容差而收斂失敗。

使用固定容差1e-6的另一個問題是它不適用於小數字。例如，假設你採用1e-16的平方根。平方根將爲1e-8。收斂性測試將錯誤地決定在第一次迭代中找到解決方案。在這種情況下，需要更小的容差。

更復雜的收斂測試可能會有意義。例如，檢查當前估計值是否比先前估計值更接近。