Python贝塞尔函数数值积分异常：结果随采样点数量放大

核心问题：变量复用导致的误解 #

你的代码中numpoints同时用于生成z_values的采样点数量和积分的采样点数量。当你增大numpoints时：

• 单个z的积分结果会趋近于真实值（而非成比例放大）
• z_values的点数同步增多，np.sum(np.square(outputspinwave))的求和项数成比例增加，最终导致总和看起来随numpoints放大

准确积分的实现方案 #

要兼顾精度和速度，可从以下几个方向优化：

1. 分离变量控制 #

将z的采样点数与积分的采样点数分开定义，避免混淆：

z_numpoints = 1000
integral_numpoints = 200  # 可单独调整积分精度

z_values = np.linspace(1e-10, 1-1e-10, z_numpoints)    
outputspinwave = spinwave_calculation(z_values, E_c_val, E_p_val, integral_numpoints)

2. 升级积分方法 #

np.trapz（梯形法）精度有限，针对不同场景可选择更高效的方法：

辛普森法（Scipy Simpson）

适合等间距采样，精度远高于梯形法，速度接近：

from scipy.integrate import simpson

def spinwave_calculation(z_values, Ec, Ep, numpoints = 200):
    wdummy_values = np.linspace(1e-10, 1-1e-10, numpoints)
    readin_values = readinKernel(wdummy_values, z_values, Ec, Ep)
    # 辛普森法对偶数采样点数支持更优，也可自动处理奇数点数
    readin_integrals = simpson(readin_values, wdummy_values, axis=1)
    return readin_integrals

自适应向量积分（Quad Vec）

如果被积函数存在局部剧烈变化，自适应积分精度最高，quad_vec支持批量处理z值，比循环调用quad效率更高：

from scipy.integrate import quad_vec

def integrand(w, z, Ec, Ep, kval=1, ic=1):
    return Ec * kval * special.jv(0, 2*Ec*kval*np.sqrt(z*(1 - w))) * np.sqrt(ic)*Ep

def spinwave_calculation(z_values, Ec, Ep):
    # 对批量z值自适应积分，返回积分结果和误差
    integrals, _ = quad_vec(integrand, 1e-10, 1-1e-10, args=(z_values, Ec, Ep))
    return integrals

Romberg积分

适合等间距采样，精度极高，但需要采样点数为2^n + 1（如33、65、129等）：

from scipy.integrate import romb

def spinwave_calculation(z_values, Ec, Ep, numpoints = 129):  # 2^7 +1=129
    wdummy_values = np.linspace(1e-10, 1-1e-10, numpoints)
    readin_values = readinKernel(wdummy_values, z_values, Ec, Ep)
    step = (1 - 2e-10)/(numpoints -1)  # Romberg需要指定采样步长
    readin_integrals = romb(readin_values, dx=step, axis=1)
    return readin_integrals

3. 变量替换优化 #

针对被积函数中的sqrt(1 - w)，做变量替换t = sqrt(1 - w)（即w=1-t²，dw=-2t dt），可消除边界数值问题，同时让被积函数更平滑：

def integrand(t, z, Ec, Ep, kval=1, ic=1):
    # 变量替换后，积分区间从t=0到1，dw=2t dt
    return Ec * kval * special.jv(0, 2*Ec*kval*np.sqrt(z)*t) * np.sqrt(ic)*Ep * 2*t

def spinwave_calculation(z_values, Ec, Ep):
    integrals, _ = quad_vec(integrand, 0, 1, args=(z_values, Ec, Ep))
    return integrals

4. 精度与速度的平衡 #

• 平滑函数场景：辛普森法+200-500个采样点，可在速度和精度间取得最优平衡
• 复杂函数场景：自适应积分quad_vec精度最高，速度略慢，适合高精度需求场景
• 可通过少量z点的高精度积分结果（如quad）作为基准，验证所选方法的精度

最终修改示例（辛普森法+分离变量） #

import numpy as np
from scipy import special
from scipy.integrate import simpson

def readinKernel(wdummy, z, Ec, Ep, kval=1, ic = 1):
    return (Ec * kval * special.jv(0, 2* Ec * kval * np.sqrt(np.outer(z, (1 - wdummy))))
            * np.sqrt(ic)*Ep )

def spinwave_calculation(z_values, Ec, Ep, integral_numpoints = 200):
    wdummy_values = np.linspace(1e-10, 1-1e-10, integral_numpoints)
    readin_values = readinKernel(wdummy_values, z_values, Ec, Ep)
    readin_integrals = simpson(readin_values, wdummy_values, axis=1)
    return readin_integrals

# 分离控制z点数和积分点数
z_numpoints = 1000
integral_numpoints = 200

z_values =  np.linspace(1e-10, 1-1e-10, z_numpoints)    
E_c_val = 1
E_p_val = 1

outputspinwave = spinwave_calculation(z_values, E_c_val, E_p_val, integral_numpoints)
output = np.sum(np.square(outputspinwave))
print(output)

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