C++でニュートン法を美しく書く その１

誰もが一度は通る道かもしれませんが，数値計算の例として，ニュートン法を書いてみたいと思います．なぜc++かというと，私がc++でいろいろと数値計算をしているからです．数値計算というとやっぱりFortranと言う気もしますが，c++はそこそこ実行速度も速いし，いろいろライブラリも整ってるし，オブジェクト指向はそこそこ便利と言うこともあって，私はC++は数値計算には結構向いていると思います．あともちろん，自分が所属している組織のルールで言語指定っていうのもあるしね．

それで，一応例題としては，

y=x^2-5

で，y=0となるときのxを求めてみます．つまり$\sqrt{5}$（2.236･･･）を求めるってことですね．
（$-\sqrt{5}$については取りあえずおいておきます．）

早速ですが，ニュートン法の第1バージョン

main.cpp

#include <iostream>

void main()
{
	double x;
	double y;

	x = 5.0;//初期値

	double y1;
	double y2;
	double x1;
	double x2;
	double diff_x;
	double error = 0.00001;
	double delta = 0.00001;

	double next_x;

	x1 = x;

	while (1) {

		x2 = x1 + delta;

		//ここがy=x^2-5
		y1 = x1*x1 - 5.0;
		y2 = x2*x2 - 5.0;

		diff_x = (y2 - y1) / (x2 - x1);
		next_x = x1 - y1 / diff_x;

		if (fabs(y1 - 0.0) < error) {
			break;
		}
		x1 = next_x;
	}

	x = x1;
	y = y1;

	std::cout << "x = " << x << std::endl;
	std::cout << "y = " << y << std::endl;

}

とりあえずで良いならこれで何とか解は求まる．
でも，いろいろ改善したいことはあって，

• 式を直打ちしているので，ニュートン法の部分が使い回せない．しかも同じ式を二回も打っていて無駄だしバグになりやすい．
• 今回はたまたま解に収束してくれたから良いけど，何回もループを回して解にたどりつかなかったらエラーを返して欲しい．

とか，挙げだしたら切りが無い．
これを少しずつクールにして行こう．
まずは，解きたい関数を直打ちするのは良くない．直そう．同じ式を何度も打つのはバグの温床になる．

main.cpp

#include <iostream>

int main()
{
	double x;
	double y;

	x = 5.0;

	auto func = [](double x) {
		return x * x - 5.0;
	};

	double y1;
	double y2;
	double x1;
	double x2;
	double diff_x;
	double error = 0.00001;
	double delta = 0.00001;
	double next_x;

	x1 = x;

	while (1) {

		x2 = x1 + delta;

		y1 = func(x1);
		y2 = func(x2);

		diff_x = (y2 - y1) / (x2 - x1);
		next_x = x1 - y1 / diff_x;

		if (fabs( y1 - 0 ) < error) {
			break;
		}
		x1 = next_x;
	}

	x = x1;
	y = y1;

	std::cout << "x = " << x << std::endl;
	std::cout << "y = " << y << std::endl;
	return 0;

}

よし，これで関数部分に$y=x^2-5$を書き出したから，汎用性は（少し）出てきた．（ちなみにラムダ式を使っているのは私の勉強も兼ねている．いろいろと新しいものを取り入れてみたいのである）
でもこのままだと，while文とかを直打ちしないといけないから，まだまだ汎用性は無い．
ということで，while文のあたりをまるまる関数にしてしまおう．

main.cpp

#include <iostream>

double solve_by_newton_method(double x ) {

	auto func = [](double x) {
		return x * x - 5.0;
	};

	double y1;
	double y2;
	double x1;
	double x2;
	double diff_x;
	double error = 0.00001;
	double delta = 0.00001;
	double next_x;

	x1 = x;

	while (1) {

		x2 = x1 + delta;

		y1 = func(x1);
		y2 = func(x2);

		diff_x = (y2 - y1) / (x2 - x1);
		next_x = x1 - y1 / diff_x;

		if (fabs(y1 - 0) < error) {
			break;
		}
		x1 = next_x;
	}
	
	return x1;

}

void main()
{
	double x;
	double y;
	double xini;

	xini = 5.0;//初期値

	x = solve_by_newton_method(xini);
	y = x * x - 5.0;

	std::cout << "x = " << x << std::endl;
	std::cout << "y = " << y << std::endl;

	return 0;
}

うーっん．main文はすっきりしたけど，今度はsolve_by_newton_method関数の中に$y=x^2-5$が直打ちされてしまった．こりゃだめだ．しかもmain文の最後で直打ちが発生してしまっている．
ってことで，なんとかして$y=x^2-5$を外に持ってこなくてはならない．もういちどfuncの定義をmain文に持ってきたい．
と言うことで次のバージョンは以下の通り．

main.cpp

#include <iostream>

auto solve_by_newton_method = []( auto func , double x ) {

	double y1;
	double y2;
	double x1;
	double x2;
	double diff_x;
	double error = 0.00001;
	double delta = 0.00001;
	double next_x;

	x1 = x;

	while (1) {

		x2 = x1 + delta;

		y1 = func(x1);
		y2 = func(x2);

		diff_x = (y2 - y1) / (x2 - x1);
		next_x = x1 - y1 / diff_x;

		if (fabs(y1 - 0) < error) {
			break;
		}
		x1 = next_x;
	}
	
	return x1;
};

int main()
{
	double x;
	double y;
	double xini;

	xini = 5.0;//初期値

	auto func = [](double x) {
		return x * x - 5.0;
	};

	x = solve_by_newton_method( func ,xini );
	y = func(x);

	std::cout << "x = " << x << std::endl;
	std::cout << "y = " << y << std::endl;
	return 0;
}

関数だと上手くかけなかったので，solve_by_newton_methodはラムダ式にしてあります．

コメントで教えてもらったんですが，

template<typename Func>
double solve_by_newton_method(Func&& func, double x, double error, double delta){

これでも大丈夫でした．

そしてmain文にラムダ式で$y=x^2-5$を書いてしまい，ラムダ式をそのままsolve_by_newton_methodに渡すことで，ニュートン法の解きたい関数をそのまま渡せるのです．

だから，ラムダ式の形で解きたい方程式を定義すれば，そのままsolve_by_newton_methodへ渡せば解ける！あとは，かゆいところに手を届かせたいので，ニュートン法のパラメータとなるdeltaとerrorも外部から入力できるようにしてしまいましょう．

main.cpp

#include <iostream>

auto solve_by_newton_method = []( auto func , double x ,double error , double delta ) {

	double y1;
	double y2;
	double x1;
	double x2;
	double diff_x;
	double next_x;

	x1 = x;

	while (1) {

		x2 = x1 + delta;

		y1 = func(x1);
		y2 = func(x2);

		diff_x = (y2 - y1) / (x2 - x1);
		next_x = x1 - y1 / diff_x;

		if (fabs(y1 - 0) < error) {
			break;
		}
		x1 = next_x;
	}
	
	return x1;
};

int main()
{
	double x;
	double y;
	double xini;

	xini = 5.0;//初期値

	auto func = [](double x) {
		return x * x - 5.0;
	};

	x = solve_by_newton_method( func ,xini , 0.00001 , 0.00001);
	y = func(x);

	std::cout << "x = " << x << std::endl;
	std::cout << "y = " << y << std::endl;
	


	auto func2 = [](double x) {
		return x * x * x - 5.0;
	};

	x = solve_by_newton_method(func2, xini, 0.00001, 0.00001);
	y = func2(x);

	std::cout << "x = " << x << std::endl;
	std::cout << "y = " << y << std::endl;

	return 0;
}

と，こんな風に，func2を定義して，solve_by_newton_methodへ渡せば解いてくれますし，ニュートン法の収束条件なども，solve_by_newton_methodの引数として与えました．
これで，そこそこ使い回せるニュートン法ができあがりました．