我試圖運行在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))
,它給出了正確的結果。
任何人都可以請幫助試圖找到我做錯了什麼?
這很有趣,因爲'的main()'似乎並沒有要求任何參數......也許嘗試['pyargs'(https://開頭www.mathworks.com/help/matlab/ref/pyargs.html)獲取_required位置參數_。哪個MATLAB版本是這個? –
@ Dev-iL,由於main()不接受任何參數,我將使用pyargs嗎?而且,如前所述,python可以很好地與這段代碼搭配使用。 Matlab版本2015a –
只需查看缺失參數列表並指定所有參數......還有什麼? –