Question 1¶

Kb = 1.38064852e-23;
%Kb = 8.6173303e-5 % ev
Z = @(E,T) sum(exp(-E./(Kb*T)));
P = @(Q,E,T) (1/Z(E,T))*exp(-E(Q)/(Kb*T));

f = 4.8e13;
h = 6.62607004e-34;
E = (0.5*f*h)*[1:0.1:101];


1a)¶

T = 300
P(1,E,T) %First state
P(2,E,T) %Second state
P(3,E,T) %Third state

T = 300 ans = 0.31883 ans = 0.21718 ans = 0.14793

1b)¶

T = 700
P(1,E,T) %First state
P(2,E,T) %Second state
P(3,E,T) %Third state

T = 700 ans = 0.15172 ans = 0.12870 ans = 0.10917
T1 = 200
T2 = 500
T3 = 800
states = 1:5;
Prob1 = [];
Prob2 = [];
Prob3 = [];
for s = states
Prob1 = [Prob1, P(s,E,T1)];
Prob2 = [Prob2, P(s,E,T2)];
Prob3 = [Prob3, P(s,E,T3)];
end

plot(states, Prob1);
hold on;
plot(states, Prob2);
plot(states, Prob3);

T1 = 200 T2 = 500 T3 = 800 Question 2¶

2b)¶

Full working in attachment under 2a)

$D(v) = (\frac{m}{2 \pi k T})2 \pi v e^{-mv^{2} / 2kT}$ $\implies$ $D(x)dx = 2xe^{-x^{2}}dx$, where x = $\frac{v}{u}$, u = $(\frac{2kT}{m})^{1/2}$

%% 2D case
maxwell2d = @(x) 2*x.*exp(-x.^2);
%% 3D case
maxwell3d = @(x) 4/sqrt(pi)*(x.^2).*exp(-x.^2);

figure(1)
ezplot(maxwell2d, [0,4])
hold on;
ezplot(maxwell3d, [0,4])
title('Maxwell speed distrobution 2D vs 3D')
legend('2D', '3D')
y = 0:0.1:4;
TwoDCharSpeed = 1/sqrt(2); % show v = u/sqrt(2) when x = 1/sqrt(2)
ThreeDCharSpeed = 1; % show v = u at x = 1
plot(0*y + TwoDCharSpeed, y, 'b')
plot(0*y + ThreeDCharSpeed, y, 'r')


ans = 1.0000 2c)¶

We see that the 2D distrobution has a lower most likely speed at v = u/sqrt(2) compared to v = u for the 3D case.
Both have the same most likely velocity vector of 0.
Both distrobution amplitudes increase as we increase x, indicating the larger multiplicites of vectors available at higher speeds.
However, after the peak likely speed, the probability amplitudes begin to fall (with skewed right tails) as the increasing multiplicity of available vector states for a given speed is compeonsated by fact that the probability of obtaining any such invididual vector decays exponentially with increasing speed.

Question 3¶

Let x = kT/ep

3a)¶

% define a vectorized version of the Boltzmann factor
J = [0:2:1000]; %even
boltz = @(x) exp( -(J.*(J+1))'./x ); % matrix / vector division
Z = @(x) (2*J+1)*boltz(x); % partition function
p = @(x) boltz(x)./Z(x); % probability of a *fixed* J, matrix / vector division
aveJ = @(x) (J.*(2*J+1))*p(x); % average J, including degeneracy term

aveE = @(x) (J.*(J+1).*(2*J+1))*p(x);
aveE_sq = @(x) (((J.*(J+1)).^2).*(2*J+1))*p(x);
sigma_E_sq = @(x) aveE_sq(x)-aveE(x).^2;
heatC = @(x) sigma_E_sq(x)./(x.^2);

figure(1)
ezplot(@(x) aveJ(x),[0,7]);
title('Average rotation level vs temperature for parahydrogen')
xlabel('kT/ep');
ylabel('ave j');
hold on;

figure(2)
ezplot(@(x) heatC(x),[0,7]);
title('Low-temperature heat capacity for parahydrogen')
xlabel('kT/ep');
ylabel('C/k');  3b)¶

% define a vectorized version of the Boltzmann factor%
J = [1:2:1001];
boltz = @(x) exp( -(J.*(J+1))'./x ); % matrix / vector division
Z = @(x) (2*J+1)*boltz(x); % partition function
p = @(x) boltz(x)./Z(x); % probability of a *fixed* J, matrix / vector division
aveJ = @(x) (J.*(2*J+1))*p(x); % average J, including degeneracy term

aveE = @(x) (J.*(J+1).*(2*J+1))*p(x);
aveE_sq = @(x) (((J.*(J+1)).^2).*(2*J+1))*p(x);
sigma_E_sq = @(x) aveE_sq(x)-aveE(x).^2;
heatC = @(x) sigma_E_sq(x)./(x.^2);

figure(1)
ezplot(@(x) aveJ(x),[0,7]);
title('Average rotation level vs temperature for orthohydrogen')
xlabel('kT/ep');
ylabel('ave j');
hold on;

figure(2)
ezplot(@(x) heatC(x),[0,7]);
title('Low-temperature heat capacity for orthohydrogen')
xlabel('kT/ep');
ylabel('C/k');  3c)¶

J = [0:2:250];
boltz = @(x) exp( -(J.*(J+1))'./x ); % matrix / vector division
Z = @(x) (2*J+1)*boltz(x); % partition function
p = @(x) boltz(x)./Z(x); % probability of a *fixed* J, matrix / vector division
aveJ = @(x) (J.*(2*J+1))*p(x); % average J, including degeneracy term
aveE = @(x) (J.*(J+1).*(2*J+1))*p(x);
aveE_sq = @(x) (((J.*(J+1)).^2).*(2*J+1))*p(x);
sigma_E_sq = @(x) aveE_sq(x)-aveE(x).^2;
heatC = @(x) sigma_E_sq(x)./(x.^2);

J = [1:2:251];
boltz = @(x) exp( -(J.*(J+1))'./x ); % matrix / vector division
Z = @(x) (2*J+1)*boltz(x); % partition function
p = @(x) boltz(x)./Z(x); % probability of a *fixed* J, matrix / vector division
aveJ = @(x) (J.*(2*J+1))*p(x); % average J, including degeneracy term
aveE = @(x) (J.*(J+1).*(2*J+1))*p(x);
aveE_sq = @(x) (((J.*(J+1)).^2).*(2*J+1))*p(x);
sigma_E_sq = @(x) aveE_sq(x)-aveE(x).^2;
heatC2 = @(x) sigma_E_sq(x)./(x.^2);

figure(1)
ezplot(@(x) (1/4)*heatC(x)+(3/4)*heatC2(x),[0,7]);
title('Low-temperature heat capacity for normal hydrogen')
xlabel('kT/ep');
ylabel('C/k');
ylim([0 1.5])
hold on;
x = 0:7;
y = 0:1;
plot(x, 0*x+1)
plot(x, 0*x+0.5)
plot(0*y + 1.65, y) %% Guess the x is around 1.6

Kb = 1.38064852e-23;
ep = 0.0076;
epj = ep*1.602e-19;
T = (1.65*epj)/Kb %Kelvin

T = 145.50

3d)¶

jodd = [1:2:1001];
jeve = [0:2:1000];
J = [jodd, jodd, jodd, jeve];
boltz = @(x) exp( -(J.*(J+1))'./x ); % matrix / vector division
Z = @(x) (2*J+1)*boltz(x); % partition function
p = @(x) boltz(x)./Z(x); % probability of a *fixed* J, matrix / vector division
aveJ = @(x) (J.*(2*J+1))*p(x); % average J, including degeneracy term

aveE = @(x) (J.*(J+1).*(2*J+1))*p(x);
aveE_sq = @(x) (((J.*(J+1)).^2).*(2*J+1))*p(x);
sigma_E_sq = @(x) aveE_sq(x)-aveE(x).^2;
heatC = @(x) sigma_E_sq(x)./(x.^2);

figure(1)
ezplot(@(x) aveJ(x),[0,7]);
title('Average rotation level vs temperature for normal hydrogen - all terms')
xlabel('kT/ep');
ylabel('ave j');
hold on;

figure(2)
ezplot(@(x) heatC(x),[0,7]);
title('Low-temperature heat capacity for normal hydrogen - all terms')
xlabel('kT/ep');
ylabel('C/k');
ylim([0 2.1])  Just for a nice figure¶

J = [0:2:1000];
boltz = @(x) exp( -(J.*(J+1))'./x ); % matrix / vector division
Z = @(x) (2*J+1)*boltz(x); % partition function
p = @(x) boltz(x)./Z(x); % probability of a *fixed* J, matrix / vector division
aveJ = @(x) (J.*(2*J+1))*p(x); % average J, including degeneracy term
aveE = @(x) (J.*(J+1).*(2*J+1))*p(x);
aveE_sq = @(x) (((J.*(J+1)).^2).*(2*J+1))*p(x);
sigma_E_sq = @(x) aveE_sq(x)-aveE(x).^2;
heatC = @(x) sigma_E_sq(x)./(x.^2);
ezplot(@(x) heatC(x),[0,7]);
hold on;

J = [1:2:1001];
boltz = @(x) exp( -(J.*(J+1))'./x ); % matrix / vector division
Z = @(x) (2*J+1)*boltz(x); % partition function
p = @(x) boltz(x)./Z(x); % probability of a *fixed* J, matrix / vector division
aveJ = @(x) (J.*(2*J+1))*p(x); % average J, including degeneracy term
aveE = @(x) (J.*(J+1).*(2*J+1))*p(x);
aveE_sq = @(x) (((J.*(J+1)).^2).*(2*J+1))*p(x);
sigma_E_sq = @(x) aveE_sq(x)-aveE(x).^2;
heatC = @(x) sigma_E_sq(x)./(x.^2);
ezplot(@(x) heatC(x),[0,7]);

J = [0:2:250];
boltz = @(x) exp( -(J.*(J+1))'./x ); % matrix / vector division
Z = @(x) (2*J+1)*boltz(x); % partition function
p = @(x) boltz(x)./Z(x); % probability of a *fixed* J, matrix / vector division
aveJ = @(x) (J.*(2*J+1))*p(x); % average J, including degeneracy term
aveE = @(x) (J.*(J+1).*(2*J+1))*p(x);
aveE_sq = @(x) (((J.*(J+1)).^2).*(2*J+1))*p(x);
sigma_E_sq = @(x) aveE_sq(x)-aveE(x).^2;
heatC = @(x) sigma_E_sq(x)./(x.^2);
J = [1:2:251];
boltz = @(x) exp( -(J.*(J+1))'./x ); % matrix / vector division
Z = @(x) (2*J+1)*boltz(x); % partition function
p = @(x) boltz(x)./Z(x); % probability of a *fixed* J, matrix / vector division
aveJ = @(x) (J.*(2*J+1))*p(x); % average J, including degeneracy term
aveE = @(x) (J.*(J+1).*(2*J+1))*p(x);
aveE_sq = @(x) (((J.*(J+1)).^2).*(2*J+1))*p(x);
sigma_E_sq = @(x) aveE_sq(x)-aveE(x).^2;
heatC2 = @(x) sigma_E_sq(x)./(x.^2);
ezplot(@(x) (1/4)*heatC(x)+(3/4)*heatC2(x),[0,7]);

jodd = [1:2:1001];
jeve = [0:2:1000];
J = [jodd, jodd, jodd, jeve];
boltz = @(x) exp( -(J.*(J+1))'./x ); % matrix / vector division
Z = @(x) (2*J+1)*boltz(x); % partition function
p = @(x) boltz(x)./Z(x); % probability of a *fixed* J, matrix / vector division
aveJ = @(x) (J.*(2*J+1))*p(x); % average J, including degeneracy term
aveE = @(x) (J.*(J+1).*(2*J+1))*p(x);
aveE_sq = @(x) (((J.*(J+1)).^2).*(2*J+1))*p(x);
sigma_E_sq = @(x) aveE_sq(x)-aveE(x).^2;
heatC = @(x) sigma_E_sq(x)./(x.^2);
ezplot(@(x) heatC(x),[0,7]);

title('Low-temperature heat capacity for hydrogen')
xlabel('kT/ep');
ylabel('C/k');
legend('Parahydrogen', 'Orthohydrogen', 'Normal hydrogen', 'Equilibrium')
ylim([0 2.1]) Question 4¶