尝试积分时指数函数溢出

Question

我想数值积分一个离散数据集（给定的ad系列）-这里是橙色-它与给定的分析指数函数（费米-狄拉克分布的导数）-这里是蓝色-相乘。然而，当指数变大（例如，对于小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])

Answer 1

这种数值问题在处理指数导数时非常常见。诀窍是首先计算日志，然后应用指数：

log(a*exp(b) / (1 + c*exp(d)) ** k) = log(a) + b - k * log(1 + exp(log(c) + d)))

现在，您需要找到一种精确计算log(1 + exp(x))的方法。幸运的是，根据这个post，人们以前做过。因此，也许您可以使用log1p重写fermi_dT：

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