読者です 読者をやめる 読者になる 読者になる

凹みTips

C++、JavaScript、Unity、ガジェット等の Tips について雑多に書いています。

boost::lambdaを使って2行で方程式の解を求める

数値計算 C++

方程式の答え簡単に出ないかなー…,と考えてるあなた.boost::lambdaを使えば「2行」で解けます.

boost::function<double(double)> f = _1*_1 - 3*_1 + 2; // f(x) = (x-1)(x-2)
std::cout << root(0.6, 1.3, f) << std::endl; // 1
/* std::cout << root(1.3, 2.4, f) << std::endl; // 2 */

ええ!?これだけ?

なぞと書いてみましたが,実際には二分法によって方程式の解を求める関数rootを作らなければならないので,もっと長くなります.しかしながら,一度この関数を作っておけば,後は利用するだけとなるので,C++数値計算をしている人にとって,コーディング時間が圧倒的に短くなるのではないでしょうか.

根を求める関数rootの実装

二分法は,Wikipediahttp://ja.wikipedia.org/wiki/%E4%BA%8C%E5%88%86%E6%B3%95)を見ると詳しいですが,簡単に説明しますと,

  1. 関数f(x)が0になりそうな前後のx(前:x1, 後:x2)を2つ選ぶ
  2. 中間値x3(=(x1+x2)/2)が解の位置からみて前後どちらにあるか調べる
  3. 前ならばx1 = x3,後ろならばx2 = x3とする
  4. 1.から繰り返す

という手順を繰り返すことによって,解の位置へ近づけていき,数値的に解を得る手法です.
これを実装した関数rootを以下に示します.

/*
@brief 与えられた2つの変数の範囲内にある根を二分法によって求める.
@param[in]  from 根の範囲(下)
@param[in]  to 根の範囲(上)
@param[in]  integrand 目的の方程式(0になる時の根)
@param[in]  accuracy 二分法の繰り返し回数
@param[in]  err 正確性
*/
template <class T, typename Functor>
T root(T from, T to, Functor equation, const int accuracy = 50, const double err = 1e-6)
{
	// accuracyをチェック
	if (accuracy <= 0) {
		std::cout << "Error! 'accuracy' must be greater than 1 (@root)." << std::endl;
		return ERROR_VALUE;
	}

	// errをチェック
	if (err <= 0) {
		std::cout << "Error! 'err' must be greater than 0 (@root)." << std::endl;
		return ERROR_VALUE;
	}

	// fromとtoが同じでないか
	if (from == to) {
		std::cout << "Error! 'from' is same as 'to'." << std::endl;
		return ERROR_VALUE;
	}

	// 計算の誤差
	const T maxErr = std::max(abs(equation(from)), abs(equation(to)))*err;

	// 間に解が無い場合は,範囲を広げて探索
	int cnt = 0;
	const int maxTry = 10;
	while (equation(from) * equation(to) >= 0) {
		// どちらかが解の場合,その値を返す
		if (abs(equation(from)) <= maxErr) {
			return from;
		}
		if (abs(equation(to)) <= maxErr) {
			return to;
		}
		if (cnt > maxTry) {
			std::cout << "Error! The solution to the equation cannot be found (@root)." << std::endl;
			return ERROR_VALUE;
		}
		T delta = to - from;
		from -= delta;
		to   += delta;
		cnt++;
	}

	T mid;
	for (int i=0; i<accuracy; i++) {
		// 中間値計算
		mid = (from + to)/2;

		// ピッタリの場合は値を返す
		if (abs(equation(mid)) <= maxErr) {
			break;
		}

		// 中間値を代入
		if (equation(from) * equation(mid) < 0) {
			to = mid;
		} else {
			from = mid;
		}
	}

	return mid;
}

丸め誤差を考えると,コードがややこしくなります…,が仕方ないですね.もっと巧いまとめ方があったら教えてください.
では,これで実際に方程式を解いてみましょう.

boost::function<double(double)> f = _1*_1 - 3*_1 + 2; // f(x) = (x-1)(x-2)
std::cout << root(2.2, 2.4, f) << std::endl; // 2
std::cout << root(10.0, 20.0, f) << std::endl; // ERROR
std::cout << root(10.0, 1.0, f) << std::endl; // 1
std::cout << root(1.0, 1.0, f) << std::endl; // ERROR
std::cout << root(2.0, 1.0, f) << std::endl; // 2

範囲外の時もまぁ上手く動いています.他にも,

boost::function<double(double)> f = bind(cos, _1/2); // f(x) = cos(x/2)
std::cout << root(0, 6, f) << std::endl; // 3.14159
boost::function<double(double)> f = bind(log, _1) - bind(sin, x); // f(x) = ln(x) - sin(x)
std::cout << root(0, 6, f) << std::endl; // 2.2191

なども計算できます.

追記

後で後悔したのですが,「boost::lambdaを使って1行で方程式を解く」の方がインパクトがあったかな,と思いました.

std::cout << root(0.6, 1.3, _1*_1 - 3*_1 + 2) << std::endl; // 1

入れ子(f(x) => f(sin(x))等)にしないのであれば,定義する必要ないし….