Maxwell-Boltzmann-Verteilung in C++

C++ bietet einige statistische Verteilung an, aber eine wichtige Verteilung der statistischen Physik – die Maxwell-Boltzmann-Verteilung wird nicht angeboten (zumindest konnte ich sie nicht finden), also habe ich mir selbst einen Generator für diese Verteilung geschrieben. Der Algorithmus ist nur wenige Zeilen lang: man würfelt ein Koordinatenpaar (Energie und Verteilungsdichte) und schaut ob der Wert der gewürfelten Verteilungsdichte kleiner als der Funktionswert der Maxwell-Boltzmann-Funktion ist. Wenn ja, dann wird der Energiewert akzeptiert, ansonsten wird weiter gewürfelt.

$$f\left(E,E_T\right)=\frac{2\mathrm{\pi }}{{\left(\mathrm{\pi }{E}_T\right)}^{2/3}}{\mathrm{e}}^{–\frac{\mathrm{E}}{{\mathrm{E}}_{\mathrm{T}}}}\sqrt{\mathrm{E}}$$

Der Algorithmus wurde in eine Klasse verpackt, die zur Instantiierung die Temperatur(Energie) und den maximalen Energiewert benötigt.

/*
	Erzeugt eine Maxwell-Boltzmann-(Energie)Verteilung

	f(Ex,Etemp)=((2*Pi)/(Pi*Etemp)^(3/2))*e^(-Ex/Etemp)*Sqrt(Ex)


	Autor:	Maxim Singer 
	WWW:	www.virtual-maxim.de
	Datum:	03 Januar 2013
*/


#include <random>

#pragma once

class MaxwellBoltzmannDistribution
{
public:
	/**
	*	Konstruktor
	*
	*	energyTemp - Energiewert, an dem sich das Maximum der Verteilung befindet. Auch als kb*T (k-Boltzmann)*(Temperatur) bekannt.
	*	energyCut  - Maximaler Energiewert, der verwendet werden soll. 
	*/
	MaxwellBoltzmannDistribution(double energyTemp, double energyCut);

	/**
	*	Gibt einen Maxwell-Boltzmann verteilten Energiewert im Intervall [0,energyCut)
	*/
	double getEnergyValue();

private:

	double energyTemp;
	double energyCut;

	std::mt19937 randEng;	// random engine	
};

/*
	Erzeugt eine Maxwell-Boltzmann-(Energie)Verteilung


	Autor:	Maxim Singer 
	WWW:	www.virtual-maxim.de
	Datum:	03 Januar 2013
*/

#include "MaxwellBoltzmannDistribution.h"

MaxwellBoltzmannDistribution::MaxwellBoltzmannDistribution(double _energyTemp, double _energyCut)
	: energyCut(_energyCut), energyTemp(_energyTemp)
{}
	

double MaxwellBoltzmannDistribution::getEnergyValue()
{
	const double PI = 3.14159265359;
	std::uniform_real<double> randUniform(0.0, 1.0);	// zufallsgenerator für reelle Zahlen

	// Normierungsfaktor der Verteilung 
	const double factor = (2*PI) / std::pow(PI*energyTemp, 1.5);

	// Würfelt so lange, bis ein Energiewert in der Verteilung liegt
	double energy = 0.0;
	double density = 0.0;
	do
	{
		energy = randUniform(randEng)*energyCut;
		density = factor * std::sqrt(energy) * std::exp(-energy/energyTemp);
	}
	while(randUniform(randEng) > density);

	return energy;
}

Einen visuellen Test kann mit folgendem Code durchführen.

/*
	Testet die Klasse MaxwellBoltzmannDistribution


	Autor:	Maxim Singer 
	WWW:	www.virtual-maxim.de
	Datum:	03 Januar 2013
*/


#include <iostream>
#include <fstream>
#include <algorithm>

#include "MaxwellBoltzmannDistribution.h"



int main()
{	
	// erstelle und zeichne eine Maxwell-Boltzmann-Verteilung

	size_t energyCut = 75;		// Kästchen auf der x-Achse
	size_t amplitude = 1000;	// größerer Wert verbessert die Statistik

	// Instanz der Verteilung erstellen
	MaxwellBoltzmannDistribution mbDist(20.0, energyCut);

	// ein Histogramm erzeugen
	std::vector<size_t> histogram(energyCut);
	for(size_t i = 0; i < energyCut; ++i)
		histogram.push_back(0);

	// Histogramm mit werten füllen
	for(size_t i = 0; i < energyCut*amplitude; ++i)
		histogram.at(static_cast<size_t>(mbDist.getEnergyValue()))++;

	// Werte im Histogramm normieren, damit sie auf den Bildschirm passen
	for_each(histogram.begin(), histogram.end(), [amplitude](size_t& v){ v/=((double)amplitude/10.0);});

	// Das Maximum im Histogramm bestimmen
	size_t maxValue = *std::max_element(histogram.begin(), histogram.end());

	// Histogramm zeichnen
	for(size_t j = 0; j < maxValue; ++j)
	{
		for(size_t i = 0; i < energyCut; ++i)
		{
			if(maxValue - histogram.at(i) < j)
				std::cout << "#";
			else 
				std::cout << " ";
		}
		std::cout << std::endl;
	}


	std::cin.get();

	return 0;
}

Führt man das Programm aus, so erscheint folgende Ausgabe.

Viel Spaß damit! =)