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.

fE,ET=2ππET2/3eEETE

Maxwell-Boltzmann-Verteilung bei drei Temperaturen.

Maxwell-Boltzmann-Verteilung bei drei Temperaturen.

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.

Test des Algorithmus für die Maxwell-Boltzmann-Verteilung

Viel Spaß damit! =)




Trackback URI | Kommentare als RSS

Einen Kommentar schreiben

XHTML: Du kannst folgende Tags verwenden: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong> <sub> <sup>

Hinweis: Ich behalte mir das Recht vor solche Kommentare, die Beleidigungen oder rechtswidrige Inhalte beinhalten erst nach einer Editierung freizugeben oder kommentarlos zu löschen. Ähnliches gilt auch für Kommentare die offensichtlich nur der Suchmaschinenoptimierung dienen.