在试图积分的同时溢出指数函数
在试图积分的同时溢出指数函数
提问于 2021-02-19 12:05:53
我想数值积分一个离散的数据集(给定的熊猫序列) -here橙色，它与给定的解析指数函数(费米-狄拉克分布的导数) -here蓝色相乘。然而，当指数变大时(例如对于小T)，当导数fermi_dT(E, mu, T)爆炸时，我就失败了。我找不到用合适的方式重写fermi_dT(E, mu, T)的方法。
下面是一个很小的例子(不是熊猫系列)，我用高斯模拟了数据集。
如果T<30.我会弄到溢出的。有没有人想出一种聪明的方法来四处走动？
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
scale_plot = 1e6
kB = 8.618292134831462e-5 #in eV
Ef = 2.0
def gaussian(E, amp, E0, sig):
return amp * np.exp(-(E-E0)**2 / sig)
def fermi_dT(E, mu, T):
return ((np.exp((E - mu) / (kB * T))*(E-mu)) / ((1 + np.exp((E - mu) / (kB * T)))**2*kB*T**2))
T = 100.0
energies = np.arange(1.,3.,0.001)
plt.plot(energies, (energies-Ef)*fermi_dT(energies, Ef, T))
plt.plot(energies, gaussian(energies, 1e-5, 1.8, .01))
plt.plot(energies, gaussian(energies, 1e-5, 1.8, .01)*(energies-Ef)*fermi_dT(energies, Ef, T)*scale_plot)
plt.show()
cum = integrate.cumtrapz(gaussian(energies, 1e-5, 1.8, .01)*(energies-Ef)*fermi_dT(energies, Ef, T), energies)
print(cum[-1])
回答 1
Stack Overflow用户
回答已采纳
发布于 2021-02-19 13:04:12
这类数值问题在处理指数导数时是很常见的。诀窍是先计算日志，然后再应用指数：
log(a*exp(b) / (1 + c*exp(d)) ** k) = log(a) + b - k * log(1 + exp(log(c) + d)))
现在，您需要找到一种精确计算log(1 + exp(x))的方法。根据这个post的说法，幸运的是，人们以前也做过这样的事情。所以也许您可以使用fermi_dT重写log1p。
import numpy as np
def softplus(x, limit=30):
val = np.empty_like(x)
val[x>=limit] = x[x>=limit]
val[x<limit] = np.log1p(np.exp(x[x<limit]))
return val
def fermi_dT(E, mu, T):
a = (E - mu) / (kB * T ** 2)
b = d = (E - mu) / (kB * T)
k = 2
val = np.empty_like(E)
val[E-mu>=0] = np.exp(np.log(a[E-mu>=0]) + b[E-mu>=0] - k * softplus(d[E-mu>=0]))
val[E-mu<0] = -np.exp(np.log(-a[E-mu<0]) + b[E-mu<0] - k * softplus(d[E-mu<0]))
return val页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/66284316
复制
相似问题
java.lang.IllegalStateException:未知标签！pos=40 poolCount = 47
JBoss EAP7.0 java.lang.IllegalStateException:未知标签！pos=3 poolCount = 20警告
无法索引类模块-info.class
jbas015852无法索引类java.lang.IllegalStateException:未知标记
Weblogic 12c无法解析类文件hikaricp-3.4.5模块-info.class
添加站长 进交流群
领取专属 10元无门槛券
AI混元助手 在线答疑
关注 腾讯云开发者公众号
洞察 腾讯核心技术
剖析业界实践案例
腾讯云开发者
