平方根( sqrt )の数値的解法②

引き続き平方根の数値的解法を見ていきたいと思います。ここではニュートン法での解法をまとめておきます。ほかの方法よりも収束が早く平方根の解法としてmath.hで定義されているのがニュートン法です。

＜ニュートン法＞

ニュートン法も二分法と同じく平方根の値を方程式の解として導き出します。

\begin{align*} f(x) = x^2 - \alpha \tag{1}\end{align*} 

\(\alpha\)を入力値として(1)の関数の\(x\)軸との交点が平方根の値となります。ここで任意の初期値 \(x_{0}\)を与え、\(x_{0}\)での式(1)の接線を求めます。

接線は一次関数\(g(x)\)として次のようにかけます。

\begin{align*} g(x) = f'(x_{0}) x  + \chi \tag{2}\end{align*}

切片\(\chi\)を求めておきましょう。\(g(x)\)は\(x=x_{0}\)の時、(1)より\(x_{0}^{2}-\alpha\)なので(2)に代入して

\begin{align*} x_{0}^{2} -\alpha = f'(x_{0})x_{0} + \chi \end{align*}

 ここで\( f'(x_{0}) = 2x_{0} \) なので、\(\chi = -(x_{0}^{2} + \alpha ) \) となります。(2)を書き直しておきます。

\begin{align*} g(x) = 2x_{0} x  - x_{0}^{2} - \alpha \tag{3} \end{align*}

 この接線とx軸との交点を\(x_1\)として\(g(x_1) = 0 \)から、

\begin{align*} 0 = 2x_{0} x_1  - x_{0}^{2} - \alpha \end{align*}

\begin{align*} x_1 = \frac{1}{2} (x_0 + \frac{\alpha} {x_{0}}) \tag{4}\end{align*}

 (4)のように漸化式の形で\(x_1\)が求まります。同様に\(x_2, x_3, \dots \) と求めていくことで解\(x_n = \sqrt{\alpha }\) が求まります。初期値\(x_0\)は二分法と同じく\(\frac{\alpha + 1}{2}\)にしておきます。ではプログラムで確かめてみましょう。

#include <iostream>

double sqrt(double alpha);

int _tmain(int argc, _TCHAR* argv[])
{
    std::cout << sqrt(2.0) << std::endl;

	return 0;
}

double sqrt(double alpha){

	double x0=0, x1=0.5*(1+alpha); 
	int n;

	for( n=0; x0!=x1; ++n ){
		x0 = x1;
		x1 = 0.5*( x0 + alpha / x0 );
	}

	std::cout << n << std::endl;

	return(x1);

}

このプログラムの結果は下記になります。

n = 5
sqrt(2) = 1.41421

 doubleで計算できるMAXまでiterationさせていますが5回しかかかっていません。二分法が54回かかりましたから、圧倒的に収束スピードが速いことが分かります。\(\sqrt{2}\)の最終桁の値が違いますが、二分法の時と同じくstd::coutの制約です。実際はもっと多くの桁まで正しい値になっています。補足ですが二次方程式の解でニュートン法を使う場合は問題になりませんが、接線を使いますので関数の形、初期値の選び方によっては収束しないので注意してください。