scipy quadrature计算时Float overflow问题求解

你需要计算的积分可表示为：
$$
\frac{\gamma}{2\pi G} \int_0^\infty \frac{\cos(sx)}{\frac{s(1+2s^2+\cosh2s)}{\sinh2s-2s} + \frac{\upsilon}{2Gh}s^2} ds
$$
其中 $x=(r-R)/h$。当前用scipy.integrate.quad实现时，当$s>300$左右，$\sinh(2s)$和$\cosh(2s)$会因指数增长触发浮点数溢出，以下是可行的解决思路和实现方案：

核心思路：渐近近似避免大指数计算 #

当$s$足够大时，$\sinh(2s) \approx \cosh(2s) \approx \frac{e^{{2s}}{2}$，因此分式项$\frac{1+2s}2+\cosh2s}{\sinh2s-2s}$会趋近于$\frac{\cosh2s}{\sinh2s} = \coth2s$，而$\coth2s$在$s$很大时几乎等于1（因为$e^{-2s}$趋近于0）。此时分母可近似为：
$$
s + \frac{\upsilon}{2Gh}s^2 = s\left(1 + \frac{\upsilon}{2Gh}s\right)
$$
基于这个近似，我们可以将积分拆分为低s区间（用原表达式计算）和高s区间（用近似表达式计算），既保证精度又避免溢出。

修改后的代码实现 #

from scipy.integrate import quad
import numpy as np

h = 0.0005
G = 1
R = 0.002
upsilon = 0.06
gamma = 0.072
B = upsilon / (2 * G)
k = B / h
S = 300  # 设定阈值，大于该值使用渐近近似

def style_s(r, gamma, upsilon, G):
    x_ = (r - R) / h
    
    # 低s区间[0, S]：使用原表达式计算
    def integrand_low(s_):
        numerator = 1 + 2 * s_**2 + np.cosh(2 * s_)
        denominator = np.sinh(2 * s_) - 2 * s_
        term1 = (numerator / denominator) * s_
        term2 = k * s_**2
        return np.cos(s_ * x_) / (term1 + term2)
    
    # 高s区间[S, ∞)：使用渐近近似后的表达式
    def integrand_high(s_):
        denom = s_ * (1 + k * s_)
        return np.cos(s_ * x_) / denom
    
    # 计算两部分积分并求和
    int_low, err_low = quad(integrand_low, 0, S)
    int_high, err_high = quad(integrand_high, S, np.inf)
    
    total_integral = int_low + int_high
    return gamma / (2 * np.pi * G) * total_integral

补充说明 #

1. 阈值选择：$S=300$是基于$\sinh(600)$已经远超浮点数范围设定的，你可以根据实际计算精度需求调整这个值，比如测试$s=200$时原分式与1的差距，只要误差足够小即可。
2. 精度验证：当$s>300$时，$\coth(2s) \approx 1 + 2e^{{-4s}$，这个修正项的数值极小（$e}{-1200}$几乎为0），因此近似带来的误差可以忽略不计。
3. 替代实现方式：也可以在单个被积函数中加入条件判断，当$s$超过阈值时自动切换到近似表达式，无需拆分积分区间，效果一致。

内容的提问来源于stack exchange，提问作者Raphael