0
$\begingroup$

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?

$\endgroup$
5
  • $\begingroup$ 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 . $\endgroup$
    – user64494
    Jul 15, 2021 at 17:47
  • $\begingroup$ {{... , {-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}, ....} $\endgroup$
    – user64494
    Jul 15, 2021 at 17:48
  • $\begingroup$ thx but that's even worse ;) $\endgroup$
    – Luke
    Jul 15, 2021 at 19:39
  • $\begingroup$ Luke (@does not work.) : Sorry, don't understand you. Can you present the results to compare? TIA and regard, $\endgroup$
    – user64494
    Jul 15, 2021 at 20:01
  • $\begingroup$ I guess the issue is that I evaluated from -20000 to 1000, you somehow lost a 0 $\endgroup$
    – Luke
    Jul 15, 2021 at 21:07

1 Answer 1

3
$\begingroup$
  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  *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.