当我使用mpmath(正函数)增大积分范围时,积分值减小

2024-05-08 14:11:33 发布

您现在位置:Python中文网/ 问答频道 /正文

我正在使用python上的mpmath库计算此积分: enter image description here

enter image description here

问题是,当我增大积分范围时,积分值会慢慢减小到0。。当被积函数为正时!!我给你一些价值观:

积分范围:[-1.17,-1.11],积分值:1.28479400889196e-8

积分范围:[-1.2,0],积分值:2.87226592130464e-29

积分范围:[-10,10],积分值:0.0

除了只选择我的被积函数不是0的范围外,发生了什么以及如何解决这个问题

这是我的代码(费米-狄拉克常数和导电率已经为PMO进行了调整,所以它不完全是我所展示的积分)

import numpy as np
from math import pi, sqrt, exp, log
from scipy.integrate import quad
import scipy.constants as phys
import matplotlib.pyplot as plt
from decimal import Decimal, getcontext
getcontext().prec = 1000

import mpmath as mp


k = (phys.Boltzmann)              # Boltzmann constant
                                  #Parenthesis left due to previous decimal conversion
m = (9.1094e-31)     #electron mass

h = (phys.Planck)     #Planck constant
h_b = h/(2 * pi)        #reduced

q = (1.602176634e-19)  
e =   (1)      
dE = 3e-3 * e

pi = (pi)


m_e = m * (1.06)   #electron effective mass
m_h = m * 0.59    #hole effective mass
A_2D = m_h / (pi*h_b**2) 


###### INPUTS
T = (4.2)
V_d = (0.02)
chi = (4.05)*e
E_g = (1.16)*e


g = 4                    #degeneracy
μ = (1e3)              #mobility, cm**2/Vs
A = A_2D

W = 1.00E-05
L = 1.00E-07

#########
E_c = 0
E_v = E_c - E_g

Phi_f = (0.55)
E_f = ( E_c - E_g/2 - e*Phi_f)
Ef_d = ( E_f - e * V_d)



#######

def fd(E,E_f):             #Fermi-Dirac distribution for holes
    res = 1/(1 + mp.exp(((E - E_f)*q/(k*T))))
    return 1-res

def sigma(E):
    if E_v-E < 2 :
        res = q * A * μ *1e-4 *  dE*q *  log( 1 + mp.exp((E_v-E)/ dE)) 
    else:
        res = - q * A * μ *1e-4 *q *(E-E_v)
    return res

def integrand(E):
    return sigma(E)*( fd(E,Ef_d) -fd(E,E_f))

X = np.linspace(E_v-5.05,5,1000)

FD = np.vectorize(fd)

Sigma = np.vectorize(sigma)



Integrand = np.vectorize(integrand)
#mp.quad
Int = mp.quad(integrand, [-10,10])
I_c = W/(q*L)*Int

fig, ax1 = plt.subplots()
ax1.plot(X,FD(X,Ef_d) - FD(X,E_f), color = 'g', label='fd stat');

ax2 = ax1.twinx()

ax2.semilogy() ;
ax2.plot(X,Sigma(X), label='sigma');
ax2.plot(X,Integrand(X), label='integrand');
ax1.legend()

ax2.legend(loc='lower right')
print('Int = ',Int)

谢谢