Stirling's approximation and paramagnetism

Stirling's approximation

How large is \(n!\), and how accurate is Stirling's approximation? To refresh our memory about Matlab, and to get a feel for just how HUGE some of these numbers are, let's make a plot of \(\log(n!)\).

n = 100;
r = 10:n;
lognfac = log(factorial(r));
plot(r,lognfac);
xlabel('n'); ylabel('log(n!)');

How many digits are there in 100! ?

floor(log10(factorial(100)))+1

ans =

158

158 digits! By contrast, \(10^{100}\) only has 100 digits. So 100! is much, much bigger than a googol. The total number of atoms in the whole universe is estimated to be \(10^{80}\) or there about. Thus 100! is roughly the square of the number of particles in universe!

Stirling's approximation is the following:

\(\log(n!) = n \log(n) - n + \frac{1}{2} \log(2\pi n) + O(1/n)\,.\)

Often only the first two terms are good enough for applications. Occasionally the third term is also relevant. How accurate is this? Let's plot the residual error

a1 = r.*log(r);
a2 = a1 - r;
a3 = a2 + .5*log(2*pi*r);

approx1 = abs( lognfac - a1 );
approx2 = abs( lognfac - a2 );
approx3 = abs( lognfac - a3 );

plot(r,approx1,r,approx2,r,approx3);
title('Absolute error for k terms in Stirling approximation');
xlabel('n'); ylabel('|log(n!) - a_n|');

Looks like the first approximation is not so good. What is the percent error for each of these approximations?

relapprox1 = approx1./lognfac;
relapprox2 = approx2./lognfac;
relapprox3 = approx3./lognfac;

plot(r,relapprox1,r,relapprox2,r,relapprox3);
title('Relative error for k terms in Stirling approximation');
xlabel('n'); ylabel('|1 - a_k / log(n!)|');

As you can see, the percent error is already quite small for the two-term approximation. By the third, it's totally negligable.

Flipping coins

Suppose we have a coin that lands heads with p=.5 (a fair coin), and we flip it 50 times.

n=50;
p=.5;

k=25;
nchoosek(n,k) % number of ways to get k heads in n flips
ans*p^k*(1-p)^(n-k) % probability of k head in n flips

% Define a function for the probability of k heads
binom = @(n,k,p) nchoosek(n,k)*p^k*(1-p)^(n-k);

% Here is a combined plot for two values, p=0.1 and p=0.5.
for p=[.1,.5]
b=zeros(n+1,1);
for k = 0:50
    b(k+1) = binom(n,k,p);
end
figure(1)
bar(0:50,b)
hold on
end
title('Probability of heads in 50 weighted coin flips (p=0.1, 0.5)');
axis([0 50 0 0.25]);
xlabel('k'); ylabel('Pr(k)');
hold off

ans =

1.2641e+14


ans =

0.1123

Two-level Paramagnet

A paramagnet is a system which only has magnetization in the presence of an external magnetic field. This is in contrast with a ferromagnet which has strong interactions between the atoms (or molecules, etc.) that can cause it to retain a magnetization even in the presence of zero external field. A simple model for a paramagnet is a collection of \(n\) non-interacting particles with two energy levels each. We are going to consider just a single particle.

The energy levels of a paramagnet are described by the energy function

\(E_{\mu} = -\vec{\mu}\cdot \vec{B} = - \mu B \cos \theta\,.\)

The two-state paramagnet can only be either aligned or antialigned with the magnetic field, so the only allowed energy levels are \(\pm \mu B\). The aligned state has the lower (\(-\)) energy.

For a single particle, we just sum the exponentials of the energies:

\(Z = \e^{-\mu B/kT}+\e^{+\mu B/kT} = 2 \cosh(\mu B/kT)\,.\)

Now suppose that our collection of paramagnets is at \(T=300K\) (room temperature), so \(kT = .026\) eV. Suppose that \(\mu = 1\) eV/T.

We can derive the first formula as follows:

\(\begin{aligned} \avg{E}\, & =(-\mu B) \e^{+\mu B/kT}/Z + (+\mu B) \e^{-\mu B/kT}/Z \\ & = (-\mu B) \frac{\e^{+\mu B/kT} - \e^{-\mu B/kT}}{2 \cosh(\mu B/kT)} \\ & = -\mu B \tanh(\mu B/kT)\,. \end{aligned}\)

The second formula is given by

\(\begin{aligned} \avg{\vec{\mu}\cdot \hat{z}}\, & = (+\mu) \e^{+\mu B/kT}/Z + (-\mu) \e^{-\mu B/kT}/Z \\ & = \mu \tanh(\mu B/kT)\,. \end{aligned}\)

We can plot these as follows:

kT = .026;
mu = 1;
B = -.1:.001:.1;
aveenergy = -mu*B.*tanh(mu*B/kT);
avemag = mu*tanh(mu*B/kT);
plot(B,avemag,B,aveenergy);
title('Average mu and E of 1eV paramagnet vs. B at 300 K');
xlabel('B (Tesla)'); legend('average \mu (eV/T)', 'average E (eV)');

Note, \(.005 T\) is the strength of a typical refrigerator magnet, and \(0.15 T\) is the strength of a sunspot.