2017-02-23 166 views
1

我試圖運行在Matlab以下Python代碼:運行Python代碼在Matlab

from numba import autojit 
import numpy as np 
import matplotlib.pyplot as plt 
import time 
from scipy.integrate import quad 
from scipy.special import jv 

def make_p(p_max):  
    p_list = [] 
    for p in range(p_max+1): 
     temp = [] 
     for u in range(p+1): 
      for r in range(p+1): 
       for i in range(p+1): 
        if u+r+i == p: 
         temp.append([u,r,i]) 
     p_list.append(temp) 
    return p_list 

@autojit 
def p_find(e,Ae,Aho,k,q,p_max,p_list): 
    bessels_sqr_0 = (jv(0,Ae*k*q/2)*jv(0,Aho*k*q/4)*jv(0,Aho*k*q/4))**2  
    if bessels_sqr_0 > 0.99: 
     return 0 
    else: 
     for i in range(1,len(p_list)): 
      bessels_sqr = 0 
      for triplets in p_list[i]: 
       n,j,l = triplets[0], triplets[1], triplets[2] 
       bessels_sqr += (jv(n,Ae*k*q/2)*jv(j,Aho*k*q/4)*jv(l,Aho*k*q/4))**2 
      if bessels_sqr_0 + bessels_sqr > 0.99: 
       return i 
      return p_max 

@autojit 
def integrand1(x,Dm,O_sqr,G_sqr,Dl,h,Ae,Aho,o_rf,o_ho,k,q): 
    sum_func = 0 
    for u in h: 
     n,j,l = u 
     sum_func += (O_sqr*(jv(n,Ae*k*q/2)*jv(j,Aho*k*q/4)*jv(l,Aho*k*q/4))**2)/(G_sqr+2*O_sqr*((jv(n,Ae*k*q/2)*jv(j,Aho*k*q/4)*jv(l,Aho*k*q/4))**2)+4*(Dl-o_rf*(n-l+j)-o_ho*(l+j)+Dm*x)**2) 
    return (x/(np.pi*np.sqrt(1-x**2)))*sum_func 

@autojit 
def dEdt1(O_sqr,G_sqr,Dl,Dm,hb,h,Ae,Aho,o_rf,o_ho,k,q): 
    return -hb*np.sqrt(G_sqr)*Dm*quad(integrand1,-1,1,args=(Dm,O_sqr,G_sqr,Dl,h,Ae,Aho,o_rf,o_ho,k,q), epsabs = 1e-35)[0] 

@autojit 
def integrand2(x,Dm,O_sqr,G_sqr,Dl): 
    return (x/(np.pi*np.sqrt(1-x**2)))*(O_sqr/(G_sqr+2*O_sqr+4*(Dl+Dm*x)**2)) 

@autojit 
def dEdt2(O_sqr,G_sqr,Dl,Dm,hb): 
    return -hb*np.sqrt(G_sqr)*Dm*quad(integrand2,-1,1,args=(Dm,O_sqr,G_sqr,Dl), epsabs = 1e-35)[0] 

def main(): 
    kb = 1.38064852e-23 
    hb = 1.054571800e-34 
    s = 0.5 
    G_sqr = (2*np.pi*21.5e6)**2 
    k = 2*np.pi/(422e-9) 
    O_sqr = s*G_sqr/2 
    dt = 1e-5 
    m = 87.9*1.660539040e-27 
    max_t = 2e-3 
    tt = np.arange(0,max_t,dt) 
    Dl = -2*np.pi*10e6 
    T0 = 1 

    E_array = np.zeros(len(tt)) 
    E_array[0] = kb*T0 

    p_max = 5 
    o_rf = 2*np.pi*26.51e6 
    a = 0 
    q = 0.1 
    o_ho = o_rf/(2*np.sqrt(a+q**2/2)) 
    Ae = 1e-9 

    p_list = make_p(p_max)   

    for i in range(len(E_array)-1): 
     if i%int(len(E_array)/10) == 0: 
      print(i/(len(E_array)-1)/10)   
     Aho = np.sqrt(2*E_array[i]/(m*o_ho**2)) 
     p_min = p_find(E_array[i],Ae,Aho,k,q,p_max,p_list) 
     if p_min > 0: 
      print(p_min) 
      for h in p_list[:p_min+1]: 
       if i < len(E_array):    
        k1 = dEdt1(O_sqr,G_sqr,Dl,k*np.sqrt(2*E_array[i]/m),hb,h,Ae,Aho,o_rf,o_ho,k,q) 
        k2 = dEdt1(O_sqr,G_sqr,Dl,k*np.sqrt(2*(E_array[i]+(dt/2)*k1)/m),hb,h,Ae,Aho,o_rf,o_ho,k,q) 
        k3 = dEdt1(O_sqr,G_sqr,Dl,k*np.sqrt(2*(E_array[i]+(dt/2)*k2)/m),hb,h,Ae,Aho,o_rf,o_ho,k,q) 
        k4 = dEdt1(O_sqr,G_sqr,Dl,k*np.sqrt(2*(E_array[i]+dt*k3)/m),hb,h,Ae,Aho,o_rf,o_ho,k,q) 
        E_array[i+1] = E_array[i]+dt/6*(k1+2*k2+2*k3+k4)     
     else: 
      if i < len(E_array):    
       k1 = dEdt2(O_sqr,G_sqr,Dl,k*np.sqrt(2*E_array[i]/m),hb) 
       k2 = dEdt2(O_sqr,G_sqr,Dl,k*np.sqrt(2*(E_array[i]+(dt/2)*k1)/m),hb) 
       k3 = dEdt2(O_sqr,G_sqr,Dl,k*np.sqrt(2*(E_array[i]+(dt/2)*k2)/m),hb) 
       k4 = dEdt2(O_sqr,G_sqr,Dl,k*np.sqrt(2*(E_array[i]+dt*k3)/m),hb) 
       E_array[i+1] = E_array[i]+dt/6*(k1+2*k2+2*k3+k4) 
    return E_array, tt, kb 


start = time.time() 
E_array, tt, kb = main() 
end = time.time() 
print('runtime: {:.3f}'.format(end-start)) 
plt.figure() 
plt.title('Temperature vs time') 
plt.plot(tt*1e3,E_array/kb, linestyle = '', marker = 'o', color='black') 
plt.show() 

的問題是,當我嘗試打電話py.matcode.main(),我收到以下錯誤信息:Python的錯誤:

main() missing 14 required positional arguments: 'E_array', 'k', 'm', 'O_sqr', 'G_sqr', 'Dl', 'hb', 'kb', 'TD', 'method', 'p_list', 'o_ho', 'Ae', and 'p_max'

即使Matlab顯示Python運行我的代碼時出現問題,代碼在Python中運行良好。

而且,I am能夠呼叫py.matcode.make_p(py.int(5)),它給出了正確的結果。

任何人都可以請幫助試圖找到我做錯了什麼?

+1

這很有趣,因爲'的main()'似乎並沒有要求任何參數......也許嘗試['pyargs'(https://開頭www.mathworks.com/help/matlab/ref/pyargs.html)獲取_required位置參數_。哪個MATLAB版本是這個? –

+0

@ Dev-iL,由於main()不接受任何參數,我將使用pyargs嗎?而且,如前所述,python可以很好地與這段代碼搭配使用。 Matlab版本2015a –

+0

只需查看缺失參數列表並指定所有參數......還有什麼? –

回答

0

更改代碼的底部:

if __name__ == "__main__": 
    start = time.time() 
    E_array, tt, kb = main() 
    end = time.time() 
    print('runtime: {:.3f}'.format(end-start)) 
    plt.figure() 
    plt.title('Temperature vs time') 
    plt.plot(tt*1e3,E_array/kb, linestyle = '', marker = 'o', color='black') 
    plt.show() 
+0

我很抱歉,但沒有幫助。它是否解決了您的Matlab版本問題? –