boost::lambdaで積分計算

boost::lambdaについて

boost::lambdaはC++でlambda expressionを可能にしてくれる，なんというか非常に気持ち悪いながら，強力且つ便利なものです．
記法については，boostのページ（http://www.boost.org/doc/libs/1_35_0/doc/html/lambda.html（自分は余り読んでませんが…））に詳しく書いてありますが，次のように説明されることが多いようです．

#include <iostream>
#include <algorithm>
#include <iterator>
#include <boost/lambda/lambda.hpp>

using namespace boost::lambda;

int main()
{
	int a[5] = {1, 2, 3, 4, 5};
	std::for_each(a, a+5, _1*=_1);
	std::copy(a, a+5, std::ostream_iterator<int>(std::cout, " "));
	return 0;
}
// 出力： 1 4 9 16 25

…，何を言ってるのかわからねーと（以下略）．って感じですよね．
中身としてはlambda functorを返す巧妙な仕組みとなっています．つまり各演算子がオーバーロードされているので，_1*_1などを行えば返り値としてlambda functorが返ってくるわけです．対して，次のようなコードはエラーとなります．

double a[5] = {1, 2, 3, 4, 5};
std::for_each(a, a+5, _1=sin(_1));
// 'sin' : 3 オーバーロードのどれも、すべての引数の型を変換できませんでした

これは，boost::lambda::bindを用いると解決できます．

double (*pSin)(double) = &sin;
std::for_each(a, a+5, _1 = bind(pSin, _1));

sinを直接書くとエラーになってしまうのはちょっと面倒ですが，bindすれば関数合成も出来ると思うと，ちょっと数値計算活用へ向けた期待が湧いてきます（※1）．

#include <iostream>
#include <cmath>
#include <algorithm>
#include <iterator>
#include <boost/function.hpp>
#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>

using namespace boost::lambda;

int main()
{
	double a[5] = {1.0, 2.0, 3.0, 4.0, 5.0};
	double (*pSin)(double) = &sin;
	double (*pCos)(double) = &cos;
	boost::function<double(double)> f = bind(pSin, _1) * bind(pCos, _1);
	std::transform(a, a+5, a, f);
	std::copy(a, a+5, std::ostream_iterator<double>(std::cout, " "));
	return 0;
}

boost::lambdaを用いた積分計算

前置きが長くなりましたが，C++でよりシンプルに積分計算を行えないかと思い，実装を行ってみました．積分は，高校のときにも習ったかもしれませんが，次の式で近似することが出来ます．

$\int_{a}^{b} f(x)dx \simeq \sum_{k=0}^{N-1} f(k \Delta x) \Delta x$
（ただし $\Delta x = \frac{b-a}{N}$ ）

しかしながら本式は誤差が大きく，例えばf(x) = x, a = 0, b = 5, N = 5としたとき，本来は 5 * 5 / 2 = 12.5 になって欲しいところが，10にしかなりません．N = 10にしても，11.25であり，1/Nでしか誤差が小さくなりません．Nを10000とか100000へ持っていければ良いのですが，コストを考えるとそこまで大きくは出来ません．
そこで，「シンプソンの積分公式」（※2）を用います．これは，kで参照している要素とk+1/2，k+1をそれぞれ重み付けして足すものです．

$\int_{a}^{b} f(x)dx \simeq \sum_{k=0}^{N-1} \left( f(k \Delta x) + 4 f(k \Delta x + \Delta/2) + f((k+1) \Delta x) \right) \frac{\Delta x}{6}$

では，これを実装した積分関数integralを実装してみます．

/*
@brief 等間隔ベクトルをセット
@param[out] y 代入結果を格納する変数
@param[in]  first 下限
@param[in]  end 上限
@param[in]  div セットする総データ数（下限〜上限までの分割数）
*/
template <class T, template <class A, class Allocator = std::allocator<A> > class Container>
void makeVector(Container<T> &x, const T first, const T end, const int div)
{
	if (div <= 0) {
		std::cout << "Error! div <= 0 (@setVector)" << std::endl;
	}
	for (int i=0; i<=div; i++) {
		x.push_back(first + ((end - first)*(T)i)/(T)div);
	}
}

/*
@brief ファンクタを用いた代入
@param[out] y 代入結果を格納する変数
@param[in]  x 代入に用いる引数
@param[in]  func y=func(x)
*/
template <class T, template <class A, class Allocator = std::allocator<A> > class Container, typename Functor>
void setVector(Container<T> &y, Container<T> &x, Functor func)
{
	y.clear();
	std::transform(x.begin(), x.end(), std::back_inserter(y), func);
}

/*
@brief 積分計算
@param[in]  from 積分範囲（下）
@param[in]  to 積分範囲（上）
@param[in]  integrand 被積分関数ファンクタ
@param[in]  N 分割数
*/
template <class T, typename Functor>
T integral(const T from, const T to, Functor integrand, const int div = 100)
{
	// N
	const int N = div*2;
	// ?x
	const T dx = (to - from) * 2 / N;
	
	// x, y
	std::vector<T> x, y;
	makeVector(x, from, to, N);
	setVector(y, x, integrand);

	// 結果
	T res = 0;
	int i = 0;
	while (i < N) {
		res += dx/6 * (y[i] + 4*y[i+1] + y[i+2]);
		i+=2;
	}

	return res;
}

では，これを実行してみます．

double res;
res = integral(0.0, pi, bind(sin, _1));
std::cout << res << std::endl; // 2
res = integral(0.0, pi, bind(cos, _1));
std::cout << res << std::endl; // 1.45717e-016

なかなか良さそうですね．
ひとつ気になるのが，なぜかsin, cosがそのままでコンパイル通るということですが…．

追記

cmathをインクルードしたときに，std名前空間の中には，オーバーロードされたstd::sinが3種類あります．
しかしながら

using namespace std;

すると，なぜか

double sin(double X)

のみになります．…なぜ？ちなみにVisual Studio 2008のインテリセンスで確認．

参考サイト

※1 http://d.hatena.ne.jp/Cryolite/20040907
boost::functionのコードを参考にさせて頂きました．というか，そもそもboost::functionをこのエントリーで知りました．感謝です．
※2 http://letsphysics.blog17.fc2.com/blog-entry-241.html
長方形近似の誤差が大きく困っていた際に「シンプソンの積分公式」を参考にさせて頂きました