# Highly oscillatory Integral over Airy function

I have to (numerically) integrate this function over all p from to get a function g(x):

f = 0.102354 Exp[(11876.9 p^2 -
154.386 p x + (0.000556163 + 0.501713 x) x)]
AiryAi[0.0275599 + 89099. p^2 + 0.00335015 x - 1158.19 p x +
3.7638 x^2]


I tried this:

tabx = Table[{x1,
NIntegrate[f /. x -> x1, {p, -Infinity, Infinity},
MinRecursion -> 10, MaxRecursion -> 20]}, {x1, -20000, 1000,
100}];


However, it takes a long time, I get a lot of error messages and the solution doesn't look like it should, especially for large negative values of x. Is there a clever way to play with methods (or any other way), to make this more accurate and faster?

• Try Table[{x1, NIntegrate[f /. x -> x1, {p, -Infinity, Infinity}, Method -> "ExtrapolatingOscillatory", AccuracyGoal -> 9, PrecisionGoal -> 9, MinRecursion -> 10, MaxRecursion -> 20]}, {x1, -2000, 1000, 100}].This results in . Jul 15, 2021 at 17:47
• {{... , {-1600, -2.34539*10^-46}, {-1500, 7.44701*10^-163}, {-1400, 0.}, {-1300, 5.80957*10^-222}, {-1200, -2.76047*10^-148}, {-1100, \ -2.64119*10^-70}, {-1000, -1.98239*10^-104}, {-900, 3.56913*10^-77}, {-800, 0.000755784}, {-700, 1.38624*10^-20}, {-600, 0.000679051}, {-500, 0.000800901}, {-400, 0.000715841}, {-300, 0.000555767}, {-200, 0.000407082}, {-100, 0.000296536}, {0, 0.000222761}, {100, 0.000177177}, {200, 0.000152183}, {300, 0.000143298}, {400, 0.00014962}, {500, 0.000174703}, {600, 0.000229498}, {700, 6.98312*10^-27}, {800, 0.000571349}, ....} Jul 15, 2021 at 17:48
• thx but that's even worse ;)
– Luke
Jul 15, 2021 at 19:39
• Luke (@does not work.) : Sorry, don't understand you. Can you present the results to compare? TIA and regard, Jul 15, 2021 at 20:01
• I guess the issue is that I evaluated from -20000 to 1000, you somehow lost a 0
– Luke
Jul 15, 2021 at 21:07

1. You should try plotting your functions before assuming they're oscillatory. Some are not. Some have numerics issues.
2. You're asking a lot if you want people to help you with 211 slow integrals. A MWE should be minimal.

Here's an idea to help with a few of the 211 integrals: Figure out where the support of the integral lies, and make sure it's sampled.

rat = Rationalize[Rationalize@#, 0] &;
integrand = Rationalize[rat@f, 0];
min = p /.
First@
Solve[
D[0.0275599 + 89099. p^2 + 0.00335015 x - 1158.19 p x +
3.7638 x^2 // rat, p] == 0, p];
dp = 1/D[
0.0275599 + 89099. p^2 + 0.00335015 x - 1158.19 p x +
3.7638 x^2 // rat, p, p];
roots = p /.
NSolve[
SetPrecision[
0.0275599 + 89099. p^2 + 0.00335015 x - 1158.19 p x +
3.7638 x^2 == 500, 32], p, WorkingPrecision -> 32];
g[x0_ /; -8.180990019212768 < x0 < 1479.8684784944994] :=
Block[{x = x0},
NIntegrate[
integrand, {p, -Infinity, min - 10 dp, min, min + 10 dp,
Infinity}, Method -> "GaussKronrodRule", WorkingPrecision -> 16]
];
g[x0_] := Block[{x = x0},
With[{central =
NIntegrate[integrand, {p, roots[[1]], Mean@roots, roots[[2]]},
Method -> "LevinRule", WorkingPrecision -> 32,
PrecisionGoal -> 6]},
central +
NIntegrate[integrand, {p, -Infinity, roots[[1]]},
Method -> "GaussKronrodRule", WorkingPrecision -> 32,
PrecisionGoal -> 6, AccuracyGoal -> 16 - Min[Log10@central, 0]] +
NIntegrate[integrand, {p, roots[[2]], Infinity},
Method -> "GaussKronrodRule", WorkingPrecision -> 16,
AccuracyGoal -> 16 - Min[Log10@central, 0]]
]];


Examples:

g[-20000]
(*  1.0286007852287060785130193502038*10^572  *)

g[-2000]
(*  37.208582874780115313397802617944  *)

g[200]
(*  0.0001521828356334408  *)
`