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

凹みTips

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

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

数値計算 C++

前のエントリ(d:id:hecomi:20100712)で,「boost::lambdaで積分計算」というソレらしい題名を付けたにも関わらず,sinやcosにバインドしただけ.というお粗末な内容になっていたので,今回はもっとガッチリとやってみます.
関係の無い人は聞き流して欲しいのですが,半導体,とりわけ化合物半導体の研究をしていると,「ボルツマン分布」の代わりに「フェルミ-ディラック分布」を使おう!というお話が出てきます.この結果,電子の統計を「フェルミ積分」というものを使って表そう!ということになります.このフェルミ積分は,以下の式で表されます.
f(E_f) = \int_0^\infty \frac{\sqrt{u}}{1+\exp \left( u - E_f \right)} du(厳密には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のお陰で,ここまで短く表現できるとは素晴らしい限りですね.