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

math.h で定義されている初等関数の数値計算に関してまとめておこうかと思います。(必ずしもmath.hで定義されているものと同じとは限りませんので注意ください)まずは一番初歩的な平方根( sqrt )の解法を幾つか紹介します。

<二分法>

平方根を方程式の問題として扱い\(f(x) = 0 \) となる\(x\)を挟み込んで段々と探す方法です。計算量が多いのですが、数値計算手法の中で最も基本的(力技で値を探す方法として基本的)なものとなります。それでは二分法で平方根を求めてみましょう。

 \begin{align*} x = \sqrt{\alpha} \tag{1} \end{align*}

 \(x\)が知りたい値です。両辺二乗して右辺を左辺に持ってきて\( x^2 - \alpha = 0 \)を得ます。この二次方程式の解が平方根の値になります。関数の形に直しておきましょう。

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

(2) は下記のような形をしています。x軸と曲線との交点を数値的に求めます。

ここからは手続きの話しになります。交点 となる\(x_{\alpha}\)に対して\(a\lt x_{\alpha}\lt b\)となる\(a, b\)を初期値として与えます。

手続き①

\(a\)と\(b\)の和の２分の１となる数値を式(2)で評価してみます。

 \begin{align*} x_{1}=\frac{a+b}{2}\end{align*}

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

手続き②

 ここで\(f(x_{1}) \lt 0\)であるならば、解 \(x_{\alpha}\)は \(x_{1} \lt x_{\alpha} \lt b \)の領域にいるはずなので

\begin{align*} x_{2}=\frac{x_{1}+b}{2}\end{align*}

 として再び式(2)で評価します。

ここで\( f(x_{1}) \gt 0\)であるならば、解 \(x_{\alpha}\)は \(a \lt x_{\alpha} \lt x_{1} \)の領域にいるはずなので 

\begin{align*} x_{2}=\frac{a+x_{1}}{2}\end{align*}

 として再び式(2)で評価します。

手続き③

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

上記のように式(2)で評価して手続き②に戻ります。

このように、手続き②と③を繰り返す事で解がどんどん近似されていきます。許容誤差 \(\beta\) を最初に決めておき \( x_{n} - x_{n-1} \lt \beta \) となった時点でiterationをやめにします。さてここで問題になるのが初期値\(a, b\)の決め方です。平方根においては次の関係を使用します。

\begin{align*} 0 \le \sqrt{x} \le \frac{x+1}{2} \end{align*} 

 平方根においてよく知られた関係です。高校数学でよく出る証明問題かと思います。では実際にプログラムを見てみましょう。

#include <iostream>

double sqrt(double alpha);

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

	return 0;
}

double sqrt(double alpha)
{
	int n;
	double a = 0, b = (alpha+1)/2;
	double x1, x2, y;

	x2 = b;
	x1 = (a + b) / 2 ;

	for(n=0;x2-x1 != 0;++n){
		x2 = x1;
		y = x1*x1-alpha; 
		if( y < 0 ){
			a = x1;
			x1 = ( x1 + b )/2;
		}else if( y > 0){
			b = x1;
			x1 = ( x1 + a )/2;
		}else{
			break;
		}
	}

	std::cout << "n = " << n << std::endl; 

	return(x1);
}

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

n = 54
sqrt(2) = 1.41421

 doubleで計算できるMAXまでiterationさせていますが54回かかっています。\(\sqrt{2}\)の最終桁の値が違いますが、これは単なるstd::coutの制約です。実際はもっと多くの桁まで正しい値になっています。

まとめ

・二分法は力技で解を探すので計算量が多いが計算機での解法の基礎が凝縮されている。