關於matlab/python轉換的建議

Question

我想把下面的matlab代碼翻譯成python代碼。該代碼計算氘分子的對數狀態的數值，然後繪製結果。當我嘗試將它翻譯爲python時，似乎我陷入了一個嵌套的for循環中，該循環計算出一筆總和。過去幾天我一直在網上搜索，但沒有成功。

因爲這是一個物理代碼，我會在代碼中提到一些方面。所以首先我們計算分區函數（Z）。之後，有一個計算能量，它是ln（Z）到β的偏導數。由此我們可以計算比熱（大約）作爲能量與溫度的導數。

所以MATLAB代碼如下所示：

epsilon = 0.0038*1.60217662*10^-19; 
k = 1.38*10^-23; 
T = 1:.1:2000; 
beta = 1./(k*T); 

%partitionfunction 
clear Z Zodd; 
for i = 1:length(T) 
    clear p;  
    for s = 1:2:31; 
     a = 2*s+1; 
     b = s^2+s; 
     p(s) = 3*a*exp(-b*epsilon*beta(i)); 
    end 
    Zodd(i) = sum(p); 
end 

%energy 
ln_Zodd = log(Zodd); 

for i = 1 : (length(T)-1) 
    Epara(i) = -(ln_Zodd(i+1)-ln_Zodd(i))/(beta(i+1)-beta(i)); 
end 

%heat capacity 
for i = 1 : (length(T)-2) 
    Cpara(i) = (Epara(i+1)-Epara(i))/(T(i+1)-T(i)); 
end 

%plot 
x = k*T/epsilon; 
plot(x(1:6000),Cpara(1:6000)/k, 'r'); 
axis([0 7 0 1.5]); 
ylabel('C_v/k'); 
xlabel('kT/eps'); 

相應的Python代碼：

import numpy as np 
import matplotlib.pyplot as plt 
import math 

epsilon=0.0038*1.60217662*10**-19 
k = 1.38*10**-23 
T = np.arange(1,2000,0.1) 
beta = 1/(k*T) 

#partitionfunction 
for i in np.arange(1,len(T)): 
    for s in np.arange(1,31,2): 
    p[s] = 3*(2*s+1)*math.exp(-(s**2+s)*epsilon*beta(i)) 
    Zodd[i] = sum(p) 

#energy 
ln_Zodd = math.log(Zodd) 

for i in np.arange(1,(len(T) - 1)): 
    Epara[i]=- (ln_Zodd(i + 1) - ln_Zodd(i))/(beta(i + 1) - beta(i)) 

#heat capacity 
for i in np.arange(1,(len(T) - 2)): 
    Cpara[i]=(Epara(i + 1) - Epara(i))/(T(i + 1) - T(i)) 

#plot 
x = k*T/epsilon 
plt.plot(x(np.arange(1,6000)),Cpara(np.arange(1,6000))/k,'r') 
plt.axis([0, 7, 0, 1.5]) 
plt.ylabel('C_v/k') 
plt.xlabel('kT/eps') 
plt.show() 

這應該是計算（近似）這個問題，因爲解析表達式的方法最簡單的方法更多地參與。我是Python的新手，所以任何建議或更正，讚賞。

Answer 1

我同意@rayryeng這個問題是離題。但是，由於我對matlab，python，和理論物理感興趣，所以我花時間仔細查看了代碼。

它有多個語法問題，還有多個語義問題。在Python中總是應該使用[]來訪問數組，通常你會嘗試使用()。數組的自然索引從0開始，與matlab不同。

這裏是原來的代碼語法和語義正確版本：

import numpy as np 
import matplotlib.pyplot as plt 
#import math #use np.* if you have it already imported 

epsilon=0.0038*1.60217662*10**-19 
k = 1.38*10**-23 
T = np.arange(1,2000,0.1) 
beta = 1.0/(k*T) #changed to 1.0 for safe measure; redundant 

#partitionfunction 
svec=np.arange(1,31,2) 
p=np.zeros(max(svec)) #added pre-allocation 
Zodd=np.zeros(len(T)) #added pre-allocation 
for i in np.arange(len(T)): #changed to index Zodd from 0 
    for s in svec: #changed to avoid magic numbers 
     p[s-1] = 3*(2*s+1)*np.exp(-(s**2+s)*epsilon*beta[i]) #changed to index p from 0; changed beta(i) to beta[i]; changed to np.exp 
    Zodd[i] = sum(p) 

#energy 
ln_Zodd = np.log(Zodd) #changed to np.log 

Epara=np.zeros(len(T)-2) #added pre-allocation 
for i in np.arange(len(T) - 2): #changed to index Epara from 0 
    Epara[i]=- (ln_Zodd[i + 1] - ln_Zodd[i])/(beta[i + 1] - beta[i]) #changed bunch of() to [] 

#heat capacity 
Cpara=np.zeros(len(T)-3) #added pre-allocation 
for i in np.arange(len(T) - 3): #changed to index Cpara from 0 
    Cpara[i]=(Epara[i + 1] - Epara[i])/(T[i + 1] - T[i]) 

#plot 
x = k*T/epsilon 
plt.plot(x[:6000],Cpara[:6000]/k,'r') #fixed and simplified array indices 
plt.axis([0, 7, 0, 1.5]) 
plt.ylabel('C_v/k') 
plt.xlabel('kT/eps') 
plt.show() 

花時間通過我提出的意見看，他們在那裏指導你。如果不清楚，請詢問澄清:)

但是，這段代碼遠沒有效率。特別是你的雙循環需要很長時間來運行（這可能解釋你爲什麼認爲它掛起）。所以我也做了非常基於numpy。

這裏的結果：

import numpy as np 
import scipy.constants as consts 
import matplotlib.pyplot as plt 

epsilon=0.0038*consts.eV #changed eV 
k = consts.k #changed 
T = np.arange(1,2000,0.1) 
beta = 1.0/(k*T) #changed to 1.0 for safe measure; redundant 

#partitionfunction 
s=np.arange(1,31,2)[:,None] 
Zodd = (3*(2*s+1)*np.exp(-(s**2+s)*epsilon*beta)).sum(axis=0) 

#energy 
ln_Zodd = np.log(Zodd) #changed to np.log 
#Epara = - (ln_Zodd[1:]-ln_Zodd[:-1])/(beta[1:]-beta[:-1]) #manual version 
Epara = - np.diff(ln_Zodd)/np.diff(beta) 

#heat capacity 
Cpara=np.diff(Epara)/np.diff(T)[:-1] 

#plot 
x = k*T/epsilon 
plt.plot(x[:len(Cpara)],Cpara/k,'r') #fixed and simplified array indices 
plt.axis([0, 7, 0, 1.5]) 
plt.ylabel('C_v/k') 
plt.xlabel('kT/eps') 
plt.show() 

請再次檢查所做的更改。我利用scipy.constants模塊將物理常數導入到高精度。我還使用了數組廣播，這使我可以將你的雙循環變成一個矩陣的一個總和（就像你應該在matlab中完成它一樣;你的原始matlab代碼也遠沒有效率）。

這裏的共同作用的結果：

你可以看到，它似乎是正確的：在高溫下你得到的獨龍族 - 佩蒂特的行爲，並在T->0我們得到了零極限按照熱力學第三定律。熱容量呈指數級下降，但這應該是有意義的，因爲您的能隙有限。