Shared2018-08-29-202059.ipynbOpen in CoCalc
StatMech

SID: 440617329 ..

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

2a)

Please see attachment

2b)

Full working in attachment under 2a)

D(v)=(m2πkT)2πvemv2/2kTD(v) = (\frac{m}{2 \pi k T})2 \pi v e^{-mv^{2} / 2kT} \implies D(x)dx=2xex2dxD(x)dx = 2xe^{-x^{2}}dx, where x = vu\frac{v}{u}, u = (2kTm)1/2(\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')

quad(maxwell2d, 0, inf)
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.

2d)

Please see attachment

2e)

Please see attachment

Question 3

Let x = kT/ep

Noting that 0 is parity even

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

Please see attachment