マクスウェル－ボルツマンの速度分布に沿った速度の与え方１：棄却法

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

で与えられます。この式はあくまで速度分布を表しているだけなので粒子の速度を与えることはできません。本項では、コンピュータで生成する疑似乱数を用いて、個々の粒子に速度の大きさを与えるアルゴリズムを考えます。

最も単純な方法は棄却法です。次のアルゴリズムで速度を決定します。
（１）乱数で速度の大きさを0から確率がほぼ０とみなせる大きさ（）までの間で決定。
（２）速度の大きさがvとなる確率を上記の式で算出。
（３）0から１の乱数を発生させて（２）の確率よりも小さければそのvを採択、大きければ不採択でもう一度（１）に戻る。

次の図は、として上記のアルゴリズムで無次元化したマクスウェル－ボルツマンの速度分布に従って粒子に速度を与えた場合の実際の速度分布です。 実線が無次元化したマクスウェル－ボルツマンの速度分布、散布図がN=10000個の粒子に速度分布です。

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

上記のアルゴリズムの場合、とdvの大きさにもよりますが、 （３）で不採択になる確率が高くなってしまうため、全ての粒子の運動量の決定までが非効率である点です。速度分布の区間分割数1000の場合、採択率は一番高いところでも0.8/1000=0.00008となり、10000個の粒子に速度を与えるのに数秒掛かってしまいます。。次回は効率よく速度を与えるアルゴリズムを考えます。

プログラムソース

//マクスウェル－ボルツマン速度分布
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 D = [];
for( var i=0; i<M; i++ ){
	D[ i ] = 0;
}

var start = new Date();
var n=0;
while( n<N ){

	//（２）速度を乱数で与える
	var v = v_max * Math.random();
	//速度分布に基づいた確率
	var p1 = f_bar( v )*dv;
	//採択を可否を決める値
	var p2 = Math.random();

	//採択の有無を決定
	if( p1 > p2 ){
		var nn = Math.floor( v/dv );
		D[ nn ]++;
		n++;
	}
}

参考図書

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

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

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

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