有人知道MATLAB和Python吗？ （代码转换MATLAB> Python）

Question

我正在尝试在Python中重新编写这个MATLAB程序。我还没有成功获得相同的Python输出。但是我的尝试是在MATLAB代码下面给出的。该代码不需要任何额外的文件/信息运行。所以这应该在你的MATLAB上运行。而且，如果一切按计划进行，巨蟒还...

摘要的代码做什么： 执行在参数不可或缺回吐eV和t，对变量eV数组。考虑替代E也有更复杂的事情。但是，那些熟悉这两个守则的人应该能够跟随。

请随时询问您是否有任何问题，并非常感谢您提供任何帮助/提示/解决方案。

MATLAB代码：

的main.m

clear all %Remove items from MATLAB workspace and reset MuPAD engine 
clc %Clear command window 
clf %Clear figure window 

global d1 d2 T %Declare global variables 

T = 0.02; %Temperature value (K) 
d1 = 1; %Energy gap in electrode 1. 
d2 = 0.5; %Energy gap in electrode 2. 

small = 1e-9; 
eV_values = linspace(0, 2.5, 2e3); %Row vector of 2e3 points linearly spaced between 0 and 0.25. These are the voltage values. 

current = zeros(size(eV_values)); %Zeros creates array all of zeros. size gives size of dataset array. 

tic %Start clock to measure performance 
for x = 1:numel(eV_values) %numel gives number of elements in array ev_values). 
    eV = eV_values(x); 

    clc %Clear command window 
    disp(x) %Display array 

    current(x) = quad(@(t)integrand(t, eV), -1 + small, 1 - small); 
end 
toc %End clock 

clf %Clear figure window 
figure(1) %Create graphics object 

hold on %Retain current graph when adding new graphs/Delay evaluation. 
box on %Display the boundary of the current axes. 
plot(eV_values, real(current), 'b') 
plot(eV_values, imag(current), 'r') 
title('S-S') 
xlabel('eV/\Delta') 
ylabel('I(eV)') 

Integrand.m

function x = integrand(t, eV) %Declare function name and inputs 

global d1 d2 T %Declare global variable 

E = t./(1 - t.^2); %Variable substitution 

x = abs(E)./sqrt(E.^2 - d1^2).*abs(E + eV)./sqrt((E + eV).^2 - d2^2).*... 
    (1./(1 + exp(E./T)) - 1./(1 + exp((E + eV)./T))); 

x = x.*heaviside(E.^2 - d1^2).*heaviside((E + eV).^2 - d2^2); 
x = x.*(1 + t.^2)./(1 - t.^2).^2; 

%heaviside step function 

Python代码：

from numpy import * 
import pylab as pl 
import array 
from scipy import integrate 

T = 0.02 # Global variable - Temperature (K) 
d1 = 1  # Global variable - Energy gap in electrode 1. 
d2 = 0.5 # Global variable - Energy gap in electrode 2. 

small = 1e-9 
eV_values = linspace(0.0, 2.5, num=10) 

def heaviside(x): 
    # Return 0 for x<0, 1 for x>0, 0.5 for x=0. # 
    if x == 0: 
     return 0.5 
    return 0 if x < 0 else 1 

def integrand(t, eV): 
    #print(" t: %s, eV: %s" % (t, eV)) 
    E = t/(1 - t*t) # E substitution. 
    x1 = (abs(E)/sqrt(E*E - d1*d1)) * (abs(E + eV)/(sqrt((E + eV)**2 - d2*d2))) * (1/(1 + exp(E/T)) - 1/(1 + exp((E + eV)/T))) 
    x2 = x1*(heaviside(E*E - d1*d1)*heaviside((E + eV)**2 - d2*d2)) 
    x = x2*((1 + t*t)/(1 - t*t)**2) 
    return x 

current = [] 
for eV in eV_values: 
    integral, err = integrate.quad(integrand, (-1 + small), (1 - small), args=(eV,)) 

# print(eV, integral) 
    print(eV, integral, err) 
    current.append(integral) 

#print('current values') 
print(current) 

#pl.plot(eV_values,current,'b') 
#pl.plot(eV_values,imag(current),'r') 
#pl.title('S-S') 
#pl.xlabel(r'eV/$\Delta$') 
#pl.ylabel('I(eV)') 
#pl.show() 

值得注意的问题：

• MATLAB代码考虑quad中的虚/实电流值。 Python代码目前没有，但应该获得相同的输出。
• 当前Python代码输出：为integral和err提供nan值。再次，这可能是由于程序没有考虑积分中的虚构值和实值。

In [3]: run IV.py IV.py:22: RuntimeWarning: invalid value encountered in sqrt x1 = (abs(E)/sqrt(E*E - d1*d1)) * (abs(E + eV)/(sqrt((E + eV)**2 - d2*d2))) * (1/(1 + exp(E/T)) - 1/(1 + exp((E + eV)/T))) IV.py:22: RuntimeWarning: overflow encountered in exp x1 = (abs(E)/sqrt(E*E - d1*d1)) * (abs(E + eV)/(sqrt((E + eV)**2 - d2*d2))) * (1/(1 + exp(E/T)) - 1/(1 + exp((E + eV)/T))) (0.0, nan, nan) (0.27777777777777779, nan, nan) (0.55555555555555558, nan, nan) (0.83333333333333337, nan, nan) (1.1111111111111112, nan, nan) (1.3888888888888888, nan, nan) (1.6666666666666667, nan, nan) (1.9444444444444446, nan, nan) (2.2222222222222223, nan, nan) (2.5, nan, nan) [nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]

Answer 1

numpy.sqrt不能处理负值。改为使用numpy.emath.sqrt。

并与浮点数计算，因为不准确的，这是更好地返回x的实部：

return x.real 

此外，numpy的标量计算是缓慢的。最好使用Python标准模块math & cmath来执行integrand()中的计算。