Einstein Solid

Considering two interacting Einstein solids:

Assuming that A and B are weakly coupled (just like ideal gas), the individual energies of the solids, q_A and q_B will change slowly. Under this assumption, the total number of energies q_total will be simply the sum of q_A and q_B.

For the first exercise, we seek to calculate the multiplicity as a function of q, when N_A = 300, N_B = 200, and q_total = 100.

import matplotlib.pyplot as plt
from math import log, exp

def factorial(z):
    if z == 0:
        return 1
    else:
        return z * factorial(z-1)
    good_z = False
    
def calc_omega1(NA,qA):
    return factorial(qA+NA-1)/(factorial(qA)*factorial(NA-1))

NA = 300
good_qA = False
while good_qA is False:        
    qA=int(input("type a value for qA energy units: "))
    if qA<0:
        print ("qA must be a positive integer!")
    else:
        good_qA = True
        
def calc_omega2(NB,qB):
    return factorial(qB+NB-1)/(factorial(qB)*factorial(NB-1))

NB = 200
qB = 100-qA

def omega_w(qA):
    return calc_omega1(NA,qA)*calc_omega2(NB,qB)
print ("NA=",NA, ", qA=",qA,",NB=",NB, ", qB=",qB, ' Multiplicity is', omega_w(qA), "")

Arbitrarily, we choose to make q to be a value of 40, which yields: type a value for qA energy units: 40 NA= 300 , qA= 40 ,NB= 200 , qB= 60 Multiplicity is 8.091930634801982e+111

For the second exercise, we seek to calculate the multiplicity as a function of q, when N_A = 3000, N_B = 2000, and q_total = 100. Using the same coding, "def factorial," will yield an error of "maximum recursion depth exceeded in comparison." To avoid such an error, we simply import "factorial" from "math" and call upon the function within the code. This method is provided in the following:

from math import factorial
def calc_omega1approx(NA,qA):
     return exp(log(factorial(qA+NA-1)) - log(factorial(qA)) - log(factorial(NA-1))) #approximation when N>>q

NA = 3000
good_qA = False
while good_qA is False:        
    qA=int(input("type a value for qA energy units: "))
    if qA<0:
        print ("qA must be a positive integer!")
    else:
        good_qA = True
        
def calc_omega2approx(NB,qB):
    return exp(log(factorial(qB+NB-1)) - log(factorial(qB)) - log(factorial(NB-1))) #approximation when N>>q

NB = 2000
qB = 100-qA

def omega_wapprox(qA):
    return calc_omega1approx(NA,qA)*calc_omega2approx(NB,qB)
print ("NA=",NA, ", qA=",qA,",NB=",NB, ", qB=",qB, ' Multiplicity is', omega_wapprox(qA), "")

Again, we choose to make q to be a value of 40, which yields: type a value for qA energy units: 40 NA= 3000 , qA= 40 ,NB= 2000 , qB= 60 Multiplicity is 6.424644088097059e+207

We now seek to provide a visual representation in the form of a plot, showing multiplicity as a function of q:

import matplotlib.pyplot as plt
from math import factorial

    
def calc_omega1approx(NA,qA):
     return exp(log(factorial(qA+NA-1)) - log(factorial(qA)) - log(factorial(NA-1))) #approximation when N>>q
        
def calc_omega2approx(NB,qB):
    return exp(log(factorial(qB+NB-1)) - log(factorial(qB)) - log(factorial(NB-1))) #approximation when N>>q

def omega_wapprox(qA):
    return calc_omega1approx(NA,qA)*calc_omega2approx(NB,qB)


NB = 2000
qB = 100-qA
NA = 3000
good_qA = False
while good_qA is False:        
    qA=int(input("type a value for qA energy units: "))
    if qA<0:
        print ("qA must be a positive integer!")
    else:
        good_qA = True  

    qA_series = range(qA+1)     
    omega_series = []

    for i in qA_series:
        omega_series.append(omega_wapprox(i))
    
    plt.plot(qA_series, omega_series)
    plt.title('qA='+str(qA))
    plt.xlabel('qA')
    plt.ylabel('$\Omega(qA)$')

    plt.show()

Acknowledgement

The original wiki page was created by J. Adams in my class