I'm trying to integrate the Airy function (of the first kind, that I'll call Ai from now) squared, between 0 and infinity.
While scipy.special.itairy() works flawlessly for integrating Ai between 0 and infinity, I haven't found a routine for Ai². I'm thus trying to calculate the integral using scipy.integrate.quad().
However, for large values of the upper bound, the value of the integral starts suddenly decreasing (from ~66.99e-3 to values close to zero), and keeps decreasing as the upper bound increases, which is of course impossible for a positive function.
I assume this is due to quad() reaching its absolute error tolerance (set by the parameter espabs, default value 1.49e-8) and simply "crashing" as a result. This is evidenced by the upper subplot on this figure.
Unlike epsabs, changing epsrel or limit doesn't seem to move the threshold after which the integral starts decreasing.
On all 3 subplots, the green curves are obtained with the default value of epsabs. The bottom left subplot shows that the integral value decreases exponentially with the upper bound once the threshold is crossed, and the bottom right subplot shows that the threshold is crossed when the upper bound reaches x ~ 5226, using the default value of espabs.
How can I overcome this problem?
The remaining part of the integral (between the threshold and +infinity) is, most likely, very close to zero and could potentially be discarded, but this feels like a jerry-rigged solution that could cause problems or result in very case-specific code, as I need to calculate the integral of Ai²(ax + b) further down the line, for different values of (a, b).
Here is the code to generate the figure :
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")
If the goal is simply to get the value of the integral, then, as per the documentation, you should use np.inf for the upper integration limit. The underlying code will use different techniques depending on the domain of the integral, and using very large numbers doesn't "activate" the techniques used for semi-infinite or infinite domains.
Here is a short example code to get the integral value:
import numpy as np
from scipy.special import airy
from scipy.integrate import quad
res = quad(lambda x: airy(x)[0]**2, 0, np.inf)
print(res) # (0.066987483779594, 1.4764521990609899e-08)
This result matches WolframAlpha's result of 0.0669875.