尝试积分时指数函数溢出

Question

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

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

如果 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

尝试积分时指数函数溢出

overflow in exponential function while trying to integrate

python

numpy

integer-overflow

scipy