Discrepancy Between MATLAB and Python

Question

I have a MATLAB script for calculating the entropy of a molecule according to its vibrational frequencies, and which I know to be correct by comparing it with experimental results. The MATLAB script is as follows:

h=6.626*10^-34;
k=1.380*10^-23;
N=6.022*10^23;
n=N/1000;
R=N*k;
P0=101325; % Abhijit's operating pressure
m=16.042*1.67*10^-27;
T=975+273.15; % Abhijit's operating temperature
B=1/(k*T);
c=3*10^8;

v1=[3015,3006,2946,2894,1460,1455,1368,1241,1232,1000,945,652,252,223,183];
f1=v1*c*100;

s4ppafin=N*sum(((h*f1)/(T*(exp((B*h*f1)-1))))-(k*log(1-exp(-h*B*f1))))

This MATLAB code produces an output of 100.8006. My attempt to reproduce this code in Python is as follows:

import numpy as np
# Do harmonic oscillator entropy approx. for the 4thppaFin system
c=3*10**8 # Approx. speed of light, m/s
h=6.626*10**-34 # Planck's constant, J/s
k=1.38*10**-23  # Boltzmann's constant, J/K
N=6.022*10**23  # Avogardo's number, particles/mole
R=N*k           # Ideal Gas Constant, J/(mol*k)
T=975+273.15    # Refernce temperature, Abhijit's operating condition
B=1/(k*T)       # Thermodynamic Beta

v=np.array([3015.0,3006,2946,2894,1460,1455,1368,1241,1232,1000,945,652,252,223,183],dtype=float)
f=np.array([c*100*i for i in v],dtype=float)
print(f)

q=[0]*len(f) # Part. func. init. for HO approx.
# q=np.array(q)

for i in range (len(f)):
  q[i]=((h*f[i])/(T*(np.exp((B*h*f[i])-1))))-(k*np.log(1-np.exp(-h*B*f[i])));

print(q)

S=N*sum(q)
print(S)

Here, the output is 145.25.

What is the cause for the Python script to produce a different number?

Answer 1

Note that your MATLAB script, in the last line:

s4ppafin=N*sum(((h*f1)/(T*(exp((B*h*f1)-1))))-(k*log(1-exp(-h*B*f1))))

has a matrix division (/). I simplify your code for clarity:

m1 = h*f1;
m2 = T*(exp((B*h*f1)-1));
m3 = k*log(1-exp(-h*B*f1));
s4ppafin = N * sum(m1/m2 - m3)

We see m1/m2. Both are 1x15 vectors. The division of the two vectors here is the least squares solution to an overdetermined system of equations:

x * m2 = m1   =>   x = m1 / m2

Because x is scalar, this is actually easy to find:

x = sum(m1 .* m2) ./ sum(m2 .^ 2)

You can compute this in Python in the same way. Note that MATLAB's ./ (element-wise division) is / in Python, MATLAB's .* is * in Python, and MATLAB's * is @ in Python.

So in Python you can write:

import numpy as np

# Do harmonic oscillator entropy approx. for the 4thppaFin system
c = 3e8 # Approx. speed of light, m/s
h = 6.626e-34     # Planck's constant, J/s
k = 1.38e-23      # Boltzmann's constant, J/K
N = 6.022e23      # Avogardo's number, particles/mole
R = N*k           # Ideal Gas Constant, J/(mol*k)
T = 975+273.15    # Refernce temperature, Abhijit's operating condition
B = 1/(k*T)       # Thermodynamic Beta

v1 = np.array([3015.0,3006,2946,2894,1460,1455,1368,1241,1232,1000,945,652,252,223,183])
f1 = v1*c*100

m1 = h*f1
m2 = T*(np.exp((B*h*f1)-1))
m3 = k*np.log(1-np.exp(-h*B*f1))
x = np.sum(m1 * m2) / sum(m2**2)
s4ppafin = N * np.sum(x - m3)
print(s4ppafin)

This produces the same output as your MATLAB code (100.8006).

Now, whether you intended to include this least squares solution in your original MATLAB code or not I cannot say, it seems to me that you didn't realize what / does in MATLAB with matrix inputs (because otherwise you would not have translated it this way in Python). Using ./ in MATLAB instead will produce the same result as your Python code (145.2506):

s4ppafin = N * sum(m1./m2 - m3)
%                    ^^