boost::lambdaで積分計算 其の二

前のエントリ（d:id:hecomi:20100712)で，「boost::lambdaで積分計算」というソレらしい題名を付けたにも関わらず，sinやcosにバインドしただけ．というお粗末な内容になっていたので，今回はもっとガッチリとやってみます．
関係の無い人は聞き流して欲しいのですが，半導体，とりわけ化合物半導体の研究をしていると，「ボルツマン分布」の代わりに「フェルミ-ディラック分布」を使おう！というお話が出てきます．この結果，電子の統計を「フェルミ積分」というものを使って表そう！ということになります．このフェルミ積分は，以下の式で表されます．
（厳密には1/2次のフェルミ積分）
本式を，これまでのエントリで出てきたCGnuplot（gnuplotでグラフを描画するためのクラス：d:id:hecomi:20100709），CLerpクラス（線形補間をするクラス: d:id:hecomi:20100710），integral関数（積分を実行する関数:d:id:hecomi:20100712）を用いて，連続関数に仕立て上げましょう．

フェルミ積分計算

コードは以下になります．

コード：

#include <vector>
#include <functional>
#include <cmath>
#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>
#include <boost/function.hpp>
#include <algorithm>
#include "calc.h" // makeVector, setVector, CLerp, integralがcalc名前空間に宣言されている
#include "gnuplot.h"

using namespace calc;
using namespace std;
using namespace boost::lambda;

int main()
{
	// サンプルポイント数
	const int N = 80;

	// 計算に使うベクトル
	vector<double> dFi, Ef;

	// サンプル座標セット（-20〜60をNこ分割 → size = N+1）
	makeVector(Ef, -20.0, 60.0, N);

	// 被フェルミ積分関数（_1: 積分変数， _2: フェルミレベル）
	boost::function<double(double, double)> f = bind(sqrt, _1) / (1 + bind(exp, _1 - _2));

	// フェルミ積分格納
	for (int i=0; i<=N; i++) {
		// フェルミレベルを代入
		boost::function<double(double)> integrand = bind2nd(f, Ef[i]);
		// 積分実行・格納
		dFi.push_back(integral(0.0, Ef[i]+200.0, integrand));
	}

	// 線形補間
	CLerp<double> fi(Ef, dFi);

	// プロット用ベクトル
	vector<double> x, sFi;
	makeVector(x, -20.0, 100.0, 1200);
	setVector(sFi, x, fi);

	// プロット
	CGnuplot gp;
	gp.SetLogPlotY(); // 片対数グラフプロットに変更
	gp.Plot(x, sFi);
	gp.SetTitle("Fermi Integral with Linear Interpolation");
	gp.SetLabel("Ef (*k_BT/q [eV])", "Fermi Integral");
	gp.DumpToEps("fermi_integral");
	int num;
	cin >> num;

	return 0;
}

出力：

プログラムの味噌

悩んだ箇所は2箇所あります．
まず1箇所目が，2変数関数を定義する部分です．

boost::function<double(double, double)> f = bind(sqrt, _1) / (1 + bind(exp, _1 - _2));

boost::lambdaを用いて，2変数を定義しました．_1が積分変数，_2が関数の変数です．bindの部分だけちょっと式が読みづらくなりますが，C++の速度で計算が実行出来ると思えば我慢できるレベルです．というかシンプルで素晴らしいです．

2箇所目は，これを1変数関数にする部分です．

boost::function<double(double)> integrand = bind2nd(f, Ef[i]);

while文の中に入っているのですが，やっていることは要は_2をEf[i]とバインドさせている（代入している）訳です．これによって被積分関数integrandを作成します．

後は，0〜無限大（Ef[i]+200で代用）まで積分することによって，あるEfの時の積分結果が得られるので，それをpush_backしてdFi(discrete:離散 fermi integralのつもり)に格納しておきます．

dFi.push_back(integral(0.0, Ef[i]+200.0, integrand));

そして，得られた結果をCLerpで線形補間すれば，オペレータ()がオーバーロードされているので，連続関数の完成となります．より細かいプロットベクトルxを用意して，CGnuplotクラスで結果を吐き出せば完成となります．

boostのお陰で，ここまで短く表現できるとは素晴らしい限りですね．