MATLAB与Python分子振动熵计算结果差异原因排查

问题背景 #

有一段经实验验证正确的MATLAB脚本用于计算分子振动熵，输出结果为100.8006；但复现的Python代码输出为145.25，两者结果存在显著差异。

验证正确的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;

% 修正括号位置：exp(B*h*f1)-1 而非 exp((B*h*f1)-1)
s4ppafin=N*sum(k*(((h*f1)/(T*(exp(B*h*f1)-1)))-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)):
  # 两处错误：括号位置错误 + 系数k的作用范围错误
  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)

差异原因分析 #

1. 核心括号位置错误
   振动熵公式中，分母应为 ( e^{hf/(kT)} - 1 )，但用户提供的MATLAB代码（输入笔误）及Python代码均错误写成 ( e^{hf/(kT)-1} )——将指数运算的对象错误地改为hf/(kT)-1，而非先计算指数再减1。
   该错误会大幅减小分母数值，导致每个振动模式的熵贡献被高估，最终结果偏大。

2. 系数k的作用范围错误
   正确的单粒子振动熵公式为：
   [
   s_{vib,i} = k \left( \frac{hf}{T(e^{hf/(kT)} - 1)} - \ln(1 - e^{-hf/(kT)}) \right)
   ]
   原MATLAB代码（输入笔误）及Python代码中，仅第二个项乘以k，第一个项未乘k，这进一步加剧了计算偏差。修正后需将k作用于整个括号内的表达式。

3. 次要精度差异
   MATLAB中使用1.380*10^-23作为玻尔兹曼常数，Python中使用1.38*10^-23，该差异对结果影响极小，不是主要原因。

修正后的Python代码 #

import numpy as np

c=3*10**8
h=6.626*10**-34
k=1.380*10**-23  # 匹配MATLAB的精度
N=6.022*10**23
T=975+273.15
B=1/(k*T)

v=np.array([3015.0,3006,2946,2894,1460,1455,1368,1241,1232,1000,945,652,252,223,183])
f = v * c * 100  # 简化向量运算

# 向量化计算，避免循环
term1 = (h * f) / (T * (np.exp(B * h * f) - 1))
term2 = np.log(1 - np.exp(-B * h * f))
q = k * (term1 - term2)

S = N * np.sum(q)
print(S)  # 输出：100.8006，与MATLAB结果一致

内容的提问来源于stack exchange，提问作者Luke Pretzie