# Plotting the maxwell-Boltzmann velocity distribution in matlab

I'm trying to plot a maxwell-boltzman velocity distribution in matlab.

I have also asked this question at cross validated without much luck

The PDF is

f(v)=sqrt(m/2*pi*k*T) * exp(-m*v^2/2*k*T)

See this link under Distribution for the velocity vector

I have defined the following function:

(This simple function is split over so many lines for debugging purposes.)

function [ output_args ] = mypdf( V_i )
m=2.18e-25;
k=1.38e-23;
T=500;

intm = m/(2*pi*k*T);
intmsqr = sqrt(intm);
exponent = -m*V_i .*V_i;
exponent = exponent/(2*k*T);
exponent = exp(exponent);
output_args = intmsqr*exponent;
end


I then ran this function for a variety of input velocities, saved them and plotted them like this:

>> velocities=-1000:1000;
>> results=mypdf(velocities);
>> plot(results)


To my horror I simply get a perfectly symmetrical distribution rather than classic maxwell-boltzman shape.

I also tried simply plotting the speeds of particles using another matlab function:

function [ output_args ] = mb_speed( V_i )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

m=2.18e-25;
k=1.38e-23;
T=500;

term1 = m/(2*pi*k*T);

term1 = term1* term1 * term1;

term1 = sqrt(term1);

term2 = 4*pi*V_i .* V_i;

term3 = -m*V_i.*V_i;

term3= term3/2*k*T;

term3 = exp(term3);

output_args = term1 * term2 .* term3;

end


I then plotted it using

>> speeds=0:1000;
>> r2=mb_speed(speeds);
>> plot(r2)


This just produces an exponential curve and not the classic shape.

What am I doing wrong in both of these cases?

Thanks.

-

First, this question has nothing to do with Matlab. Besides, if you think there is an error in your code implementation, this is most certainly the wrong stackexchange.

That said, the error is a physical one. You are confusing the distribution of velocity vectors $\vec{v}$, which is a 3D Gaussian (by that I mean a function from $\mathbb{R}^3$ to $\mathbb{R}$ that is the product of Gaussians in three Cartesian coordinates) with the 1D distribution of those vectors' magnitudes. The former had better be symmetric about 0 - do you expect that gas in thermal equilibrium is moving in one particular direction? This is in fact exactly what you plot (well, actually you plot a 1D slice of the 3D domain, essentially the distribution of $x$-velocities).

You need a $v^2$ factor in the distribution (and a few modified constants) to get the asymmetric thing you are looking for. This is the same $v^2$ (or $r^2$) that arises, for instance, when converting a spherically symmetric integrand from Cartesian to spherical coordinates. When you do that, by the way, don't try to plug in negative values - they don't make physical sense, since again it is a distribution of magnitudes.

In terms of the notation in that Wikipedia article, the distribution of vectors is $f_\mathbf{v}$, which is $f_v^3$, and you are plotting $f_v$. What you want to plot is what the article calls unsubscripted $f$ (well, one of the many distinct things it calls $f$): $$f(v) = \sqrt{\left(\frac{m}{2\pi kT}\right)^3} 4\pi v^2 \exp\left(\frac{-mv^2}{2kT}\right).$$

-
Thanks; What I eventually want to do is take the 1D Velocity (the symmetric distribution I plotted), turn it into a CDF and use it to assign velocities to particles in a MD simulation. Will it be enough to sample from it three times for each particle (and in so doing get V_x,V_y,V_z for each)? – RRs_Ghost Feb 15 '13 at 9:48
Yes - unless you have some anisotropy (e.g. a magnetic field permeating the region), $x_x$, $x_y$, and $x_z$ are all on equal footing and are thus IID variables. – Chris White Feb 15 '13 at 17:40

The first distribution (velocity) you describe should be symmetric, as detailed above.

In the second (speed), I believe you are simply mistaking your distribution to be exponential because of an order of operations mistake in your code. You have:

term3= term3/2*k*T;


What you should have is:

term3= term3/(2*k*T);


As you have it now, the already very small value of m is being multiplied by k, which is also very small. You'll notice that if you plot $y=x^2 e^{-Cx^2}$ where $C$ is very close to zero, it will look like an exponential plot until $x$ becomes sufficiently large.

-
Thanks for the code correction :) (I don't have enough rep to upvote) – RRs_Ghost Feb 15 '13 at 9:50