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

凹みTips

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

boost::lambdaで積分計算

数値計算 C++

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
長方形近似の誤差が大きく困っていた際に「シンプソンの積分公式」を参考にさせて頂きました