マクスウェル－ボルツマンの速度分布に沿った速度の与え方２：累積分布関数法

N個の粒子にマクスウェル－ボルツマン速度分布に従った速度を与えるアルゴリズムを解説します。 前回解説した棄却法では試行回数のに対して確率的にほとんどが不採択となってしまうため、必要個数の速度を与えるまでに時間がかかってしまいます。そこで本項ではより効率的なアルゴリズムを考えます。

考え方

１個の粒子の速度の大きさがvからv+dvの間に存在する確率は前回の「無次元化したマクスウェル－ボルツマン速度分布」で示したとおり

で与えられます。これは確率密度を表しているのでv=0から無限大まで積分すると１になります。 そこで、マクスウェル－ボルツマン速度分布に対する累積分布関数を

と定義し、擬似乱数で得られた値に対応するを決定する方法で粒子の速度を与える方法を考えます。次の図の青線がマクスウェル－ボルツマン速度分布に対する累積分布関数です。

図：無次元化したマクスウェル－ボルツマンの速度分布に従った速度を与えた例

上記の累積分布関数を用いて速度を与える手順は以下の通りです。
（１）あらかじめ、累積分布関数の値をから まで数値的計算して配列に格納しておく。
（今回は速度分布の分割数を1000とするので、配列要素数は0から1000までの1001個となります）
（２）乱数で0から1までの小数を決定。
（３）この小数に対応する累積分布関数値が格納されている配列の要素番号を取得し、そこから速度の大きさを見積る。

この方法の場合、あらかじめ累積分布関数の値を用意する必要がありますが、棄却法のような無駄となるプロセスが無くなるため（乱数の発生回数は粒子数の個数分だけ）、効率的に速度を決定することができます。実際、棄却法よりも1/30程度の時間で速度を得ることができました。なお、数値積分はシンプソン法を用いています。

プログラムソース

//マクスウェル?ボルツマン速度分布
function f_bar( v_bar ){
	return  4.0 / Math.sqrt(Math.PI) * v_bar * v_bar * Math.exp( -v_bar*v_bar );
}

//粒子数
var N = 100000;
//速度分布分割数
var M = 1000;
var v_min = 0;
var v_max = 4;
var dv = (v_max - v_min)/M;

//データ配列
var data1 = [];
var I = []; //累積積分値を格納する配列
var D = [];	//散布図作成用配列
for( var j=0; j<M; j++ ){
	var v = v_min + (v_max - v_min)/M * j;
	D[ j ] = 0;
	I[ j ] = SimpsonIntegral( f_bar, v_min, v, M ); //累積積分値の計算
	data1.push([ v, I[ j ] ]);
}
plot2D.pushData( data1 ); //データ１

//累積積分値の最後を１とする
I[ M ] = 1;

var start = new Date();
for( var i=0; i<N; i++ ){

	var nv = Math.random();

	for( var j=0; j<M; j++){

		if( I[ j ] < nv && nv < I[ j+1 ] ) {

			D[ j ]++;
			//v = v_min + (v_max - v_min)/M * j;
			break;

		}

	}

}

参考図書

３次元グラフィックス関連

・three.jsによるHTML5グラフィックス上　【改定版】
・three.jsによるHTML5グラフィックス下　【改定版】
・three.jsによるHTML5 3Dグラフィックス 【新機能と応用】

物理シミュレーション関連

・HTML5による物理シミュレーション
・HTML5による物理シミュレーション【拡散・波動編】
・HTML5による物理シミュレーション　【剛体編１】
・HTML5による物理シミュレーション　【剛体編２】