如何结合scipy interp1d与mpmath quadosc实现高振荡积分

针对你遇到的mpmath.quadosc不支持数组输出的问题，以及高振荡积分失效的情况，提供以下几种可行方案：

方案1：逐个计算每个q对应的积分（适配mpmath.quadosc） #

quadosc不支持返回数组类型的结果，因此可以循环遍历每个q值，单独计算积分。这样既利用了quadosc的振荡积分优化，又避开了数组输出的限制。

修改后的代码示例：

import numpy as np
from mpmath import besselj, sqrt, pi, besseljzero, inf, quadosc
from scipy.interpolate import interp1d

n = 1
q = np.geomspace(1e-7, 500, 1000)
# 构造模拟密度函数
x = np.geomspace(1e-7, 10, 1000)
y = np.exp(-(x-5)**2)
density = interp1d(x, y, kind='cubic', fill_value=0, bounds_error=False)

def spherical_jn_density_single(x, q_val, n=n):
    arg = q_val * x
    return besselj(n + 1/2, arg) * sqrt(pi / (2 * arg)) * density(x)

# 逐个计算每个q对应的积分
vals_density = []
for q_val in q:
    val = quadosc(
        lambda x: spherical_jn_density_single(x, q_val),
        [0, inf],
        zeros=lambda m: besseljzero(n + 1/2, m)
    )
    vals_density.append(float(val))
vals_density = np.array(vals_density)

方案2：使用SciPy的向量值振荡积分（更高效） #

SciPy的scipy.integrate.quad_vec支持向量值被积函数，且可以通过指定振荡函数的零点来优化积分精度。结合scipy.special.jn_zeros获取球贝塞尔函数的零点，替代mpmath的实现：

import numpy as np
from scipy.integrate import quad_vec
from scipy.interpolate import interp1d
from scipy.special import jn, jn_zeros

n = 1
q = np.geomspace(1e-7, 500, 1000)
x = np.geomspace(1e-7, 10, 1000)
y = np.exp(-(x-5)**2)
density = interp1d(x, y, kind='cubic', fill_value=0, bounds_error=False)

# 用SciPy实现球贝塞尔函数（避免mpmath的类型问题）
def spherical_jn(x, n=n):
    return np.sqrt(np.pi/(2*x)) * jn(n + 1/2, x)

def integrand(x):
    arg = q[:, None] * x
    return spherical_jn(arg) * density(x)

# 获取球贝塞尔函数的前N个零点（用于振荡积分优化）
num_zeros = 200
zeros = jn_zeros(n + 1/2, num_zeros)
# 由于积分区间是[0, inf]，截断到密度函数可忽略的区间（比如x=10，因为你的密度在x>10时几乎为0）
truncated_zeros = zeros[zeros <= 10 * q.max()]  # 对应x=zeros/q <=10

# 执行向量积分，指定振荡零点
result, error = quad_vec(
    integrand,
    0, 10,  # 截断区间，因为密度在x>10时已衰减到0
    points=truncated_zeros / q.max()  # 转换为x区间的零点
)

方案3：利用密度衰减特性截断积分区间+自适应振荡积分 #

由于密度函数在超过阈值后快速衰减，可以直接截断积分区间（比如x=10），然后使用SciPy的quad手动处理振荡特性，对每个q单独计算：

import numpy as np
from scipy.integrate import quad
from scipy.interpolate import interp1d
from scipy.special import jn

n = 1
q = np.geomspace(1e-7, 500, 1000)
x = np.geomspace(1e-7, 10, 1000)
y = np.exp(-(x-5)**2)
density = interp1d(x, y, kind='cubic', fill_value=0, bounds_error=False)

def spherical_jn(x, q_val, n=n):
    arg = q_val * x
    return np.sqrt(np.pi/(2*arg)) * jn(n + 1/2, arg)

vals_density = []
for q_val in q:
    # 截断积分区间到x=10，因为密度在此处已可忽略
    val, err = quad(
        lambda x: spherical_jn(x, q_val) * density(x),
        0, 10,
        epsabs=1e-8, epsrel=1e-8
    )
    vals_density.append(val)
vals_density = np.array(vals_density)

关键说明 #

• 方案1兼容mpmath的quadosc，但循环计算效率较低，适合q数量不多的场景。
• 方案2使用SciPy的向量积分，效率更高，适合大规模q网格，且利用零点优化振荡积分精度。
• 方案3最简单，利用密度衰减截断区间，适合密度衰减极快的场景，精度依赖截断区间的选择。

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