解决scipy.integrate.quad()积分Airy函数平方的数值异常问题

我需要计算第一类Airy函数（下文记为Ai）的平方在0到正无穷区间的积分。scipy.special.itairy()可完美计算Ai的积分，但暂无Ai²的现成积分函数，因此尝试用scipy.integrate.quad()计算。然而当积分上限增大到一定值时，积分值会从约66.99e-3突然下降至接近0，且随上限继续增大持续降低——这对于正函数的积分而言显然不合理。

推测该问题源于quad()达到了绝对误差容限参数epsabs（默认值1.49e-8），改变epsrel或limit参数无法改变积分值开始下降的阈值。默认epsabs下，阈值出现在x≈5226，此后积分值呈指数下降。

虽然阈值到正无穷的剩余积分值接近0，可直接截断，但这属于临时方案，后续还需计算不同(a,b)下Ai²(ax+b)的积分，需要通用解决方法。

生成相关图像的代码 #

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
from scipy.special import airy

fig = plt.figure()

ax2 = plt.subplot(223)
ax3 = plt.subplot(224)
ax1 = plt.subplot(211)

# subplot 1
xmax1 = np.linspace(1000, 10000, 100)
for error,clr in zip([10e-5, 10e-7, 1.49e-8, 10e-11, 10e-13, 10e-16],
                     ['red','darkorange','forestgreen','cornflowerblue','blue']
                     ):
    airy_square_int1 = [integrate.quad(lambda x: (airy(x)[0])**2, 
                                       0, 
                                       end, 
                                       epsabs = error,
                                       # epsrel = 10**(-15),
                                       limit = 20
                                       )[0] 
                        for end in xmax1
                        ]

    ax1.plot(xmax1, 
             airy_square_int1, 
             label ='epsabs = {}'.format(error),
             color = clr
             )

ax1.set_xlabel("x")
ax1.set_ylabel("integral of Ai² from 0 to x")
plt.legend()

# subplot 2
xmax2 = np.linspace(6000, 20000, 100)
airy_square_int2 = [integrate.quad(lambda x: (airy(x)[0])**2,
                                   0, 
                                   end
                                   )[0] 
                    for end in xmax2
                    ]

ax2.plot(xmax2, airy_square_int2, color = 'forestgreen')
ax2.set_xlabel("x")
ax2.set_ylabel("integral of Ai² from 0 to x")

# subplot 3
xmax3 = np.linspace(5226, 5227, 100)
airy_square_int3 = [integrate.quad(lambda x: (airy(x)[0])**2, 
                                   0, 
                                   end
                                   )[0] 
                    for end in xmax3
                    ]

ax3.plot(xmax3 - 5226, airy_square_int3, '+', color = 'forestgreen')
ax3.set_xlabel("x - 5226")

plt.tight_layout(pad = 0.5)
plt.subplots_adjust(hspace = 0.5)
plt.savefig("integral_AiSquared_behavior.png")

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