1.0E21と1.0E23の平方根を計算するとき、ニュートンラフソン法が収束しないのはなぜですか？

Question

私は、ニュートン・ラフソン法を使って数値の平方根を計算するために以下の関数を持っています。

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); 
} 

、ここでは（GCC 4.8.4を使用して、64ビットのLinux上に構築された）出力です：

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です。

IEEE-754 encoding（1e23）は、0x44b52d02c7e14af6です。したがって、指数は0x44bで1099です。指数はオフセットバイナリとして符号化されているため、これを1023だけ減らす必要があります。次に、LSBの重みを見つけるために別の52を引く必要があります。 LSB

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

だから、単一のコード

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

から値を16777216

結果を持っているが、実際にある-16777216。その結果、

、あなたがx*x - yを計算する場合、結果はどちらかでゼロになるだろうか、それが2*3.16e11あるで16777216を割る16777216の倍数になるだろうことは、あなたの中にされていない、0.000027の誤差を与えます耐性。 sqrt(y)へ

2つの最も近い値は

0x4252682931c035a0 
0x4252682931c035a1 

であり、あなたは、これらの数字を二乗するとき、あなたは

0x44b52d02c7e14af5 
0x44b52d02c7e14af7 

を取得するので、どちらもが

0x44b52d02c7e14af6 

で所望の結果が、一致しました

であり、アルゴリズムは収束することができません。

---

アルゴリズムが1e25で動作する理由は運が原因です。 1e25のエンコーディングが

0x45208b2a2c280291 

あるsqrt(1e25)のエンコーディングは

0x428702337e304309 

であり、あなたがその数を二乗するとき、あなたはx*x - yしたがって

0x45208b2a2c280291 

を取得し、正確に0、およびアルゴリズムであり、幸運になる。

Answer 2

これは、固定されたtoleranceの使用に関連して、浮動小数点精度の問題です。

減算される数値の大きさが十分に大きければ、（1）「偽」のゼロ差（これは問題ではない）、または（2）その差が公差よりもずっと大きいため収束に失敗した。

固定公差1e-6を使用する際のもう1つの問題は、小さな番号に対しては機能しないことです。たとえば、平方根が1e-16であるとします。平方根は1e-8になります。コンバージェンステストは最初の反復で解決策が見つかったと誤って判断します。この場合、はるかに小さい公差が必要となる。

より洗練されたコンバージェンステストはおそらく意味をなさないでしょう。たとえば、現在の見積もりが以前の見積もりより近いかどうかを確認します。