MATLAB 和 Python 之间的差异

Question

我有一个 MATLAB 脚本，用于根据振动频率计算分子的熵，通过与实验结果进行比较，我知道它是正确的。 MATLAB脚本如下：

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))))

此 MATLAB 代码产生 100.8006 的输出。我尝试在 Python 中重现这段代码如下：

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)

在这里，输出是 145.25.

Python脚本产生不同数字的原因是什么？

Answer 1

注意你的 MATLAB 脚本，在最后一行：

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

有一个矩阵除法（

/

）。为了清楚起见，我简化了您的代码：

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

我们看到

m1/m2

。两者都是 1x15 向量。这里两个向量的除法是超定方程组的最小二乘解：

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

因为

x

是标量，这个其实很好找：

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

您可以用相同的方式在 Python 中计算它。注意MATLAB的

./

（按元素划分）在Python中是

/

，MATLAB的

.*

在Python中是

*

，MATLAB的

*

在Python中是

@

所以在 Python 中你可以这样写：

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)

这会产生与您的 MATLAB 代码 (100.8006) 相同的输出。

现在，无论您是否打算在原始 MATLAB 代码中包含这个最小二乘解，我都不能说，在我看来，您似乎没有意识到

/

在 MATLAB 中用矩阵输入做了什么（因为否则你不会有在 Python 中以这种方式翻译它）。在 MATLAB 中使用

./

将产生与您的 Python 代码 (145.2506) 相同的结果：

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