尝试积分时指数函数溢出

Question

我想对离散数据集（给定的广告熊猫系列）进行数值积分 - 此处为橙色 - 乘以给定的分析指数函数（费米-狄拉克分布的导数）-此处为蓝色-。但是，当指数变大（例如对于小 T）以及导数 fermi_dT(E, mu, T)explodes 时，我会失败。我找不到以适当的方式重写fermi_dT(E, mu, T)来完成它的方法。

下面是一个最小的例子（不是熊猫系列），我用高斯模拟了数据集。

如果 T

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