2015年03月01日

C#で立方根を求める(ニュートン法による近似解)

   このエントリーをはてなブックマークに追加 Clip to Evernote
どう書く?orgに感謝を込めて」シリーズ その73

■問題 (出題者:にしお さん)
xは0以上1000未満の実数です。 y * y * y = xになるような実数y(立方根)を小数点以下12桁以上の正確さで 求める関数cube_rootを作って下さい。
ただし、このお題の趣旨は実数区間での探索なので、 立方根関数があっても使ってはいけません。 指数関数と対数関数も禁止します。
Pythonで表現した入出力の例:
>>> cube_root(10.0)
2.1544346900318834
>>> _ ** 3
9.9999999999999947
>>> cube_root(100.0)
4.6415888336127793
>>> _ ** 3
100.00000000000003


まずは、立方根の前にイメージしやすい平方根について考ます。

※ ここでは、 x^n は、x の n乗を意味するものとします。

√aを求めることを考えて見ます。(√は、平方根の記号のつもり)
√aは、x^2 = a であるxの値ですので、
x^2 - a = 0 の方程式を解けば答えを求めることでができます。
y = f(x) = x^2 - a のグラフで考えてみると、y座標が0 である xの座標を求めればよいことになります。

では、どのようにして求めるかですが、コンピュータを使う場合は、ニュートン法を利用する方法が知られています。 こちら

考え方はこうです。
まず、解に近い値を x0 としたときに、座標 (x0,f(x0)) との接線を引き、その接線とx軸との交点をx1 とする。
同様に、この座標 (x1, f(x1))の接線を引き、x2を求める、続いて、(x2, f(x2))の接線を引き、x3を求める。
これを繰り返すことで, xnは、解に近づいてゆく、というものです。
xn と そのひとつ前の x の値の差が十分小さくなれば、解の近似値とみなせます。

図が無いので分かりにくいかもしれませんが、ここにある図3を見ていただければ、考え方は理解できると思います。


n 番目の x を xn 、n+1番目の x を xm とすれば、(x(n+1) だと紛らわしいので)

xm = xn - w

と表せます。このときの wは、xn の接線の傾きを g(xn) とすると、

w = f(xn) / g(xn)

です。

※ 傾き = 高さ/ 横幅 ですから、横幅 = 高さ / 傾き となります。

つまり、

xm = xn - (f(xn) / g(xn))

が導き出せます。

次に、グラフ f(x) の接線の傾き、g(xn) を求めます。

まず、y = f(x) のグラフの以下の2つの点を通る直線の傾きを考えてみます。
1: (x , f(x))
2: ((x+b) , f(x+b))

傾きは、

(f(x+b)-f(x)) / b

ですから、これを、平方根の式である f(x) = x^2 - a に当てはめてみると、

( f(x+b)-f(x) ) / b
= ( (x+b)^2 - a - (x^2 - a) ) / b
= ( x^2 + 2bx + b^2 - a - x^2 + a ) / b
= ( 2bx + b^2 ) / b
= ( 2x + b )


となります。ここで、b をどんどん小さくし、0に近づけていったときの値が、 x 地点の接線の傾きとなります。
つまり、

2x が傾きです。

よって、

xm = xn - (f(xn) / g(xn))
= xn - (xn^2 - a) / 2xn

※ xn でひとつの変数です。(x * nではありませんので注意してください)

となり、これを使って次々と新しい xm を求めていけばよいことになります。変化する値が十分に小さくなれば、 この xm が√a の 近似値となります。

これをコード化にしたものを以下に示します。

static double SquareRoot(double a, double xn = 1.0) {
    double e = (xn * xn - a) / (2 * xn);
    if (Abs(e) > 1.0e-14)
        return SquareRoot(a, xn - e);
    return xn - e;
}

static double Abs(double n) {
    return n < 0 ? -n : n;
}


√20 を求める場合は、以下のようにします。

Console.WriteLine(SquareRoot(20.0));

同様の考えで、立方根についても

g(x)
= ( f(x+b)-f(x) ) / b
= ( (x+b)^3 - a - (x^3 - a) ) / b
= ( x^3 + 3x^2b + 3xb^2 + b^3 - x^3 ) / b
= ( 3x^2b + 3xb^2 + b^3 ) / b
= ( 3x^2 + 3xb + b^2 )
となり、b を 0 にして、g(x) = 3x^2 となるので、

xm = xn - (xn^3 - a) / 3xn^2

が成り立ちます。
※ xn でひとつの変数です。(x * nではありませんので注意してください)

ここまで理解すれば。C#のコードにするのは難しくないですね。

static double CubeRoot(double a, double xn = 1.0) {
    double e = (xn * xn * xn - a) / (3 * xn * xn);
    if (Abs(e) > 1.0e-13)
        return CubeRoot(a, xn - e);
    return xn - e;
}

static double Abs(double n) {
    return n < 0 ? -n : n;
}

再帰使ってますが、ループ処理にして以下のように書いても答えは求まります。

static double CubeRoot(double a) {
    double x = 1.0;
    double e;
    do {
        e = (x * x * x - a) / (3 * x * x);
        x = x - e;
    } while (Abs(e) > 1.0e-13);
    return x;
}



 

この記事へのトラックバックURL

http://trackback.blogsys.jp/livedoor/gushwell/52407167