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

方程式の答え簡単に出ないかなー…，と考えてるあなた．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の実装

二分法は，Wikipedia（http://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))等）にしないのであれば，定義する必要ないし…．