Yet another slow plot involving numerical integration

Question

Well, after reading a lot about how plotting expressions involving NIntegrate can take a lot of time, and how to overcome this issue with DSolve, I still have problems when plotting this function:
myfunction[delta_] = 
 5*(2 + 3*NIntegrate[(E^(-(1/2) (-((3 Sqrt[5/2] delta)/(-1 + delta)) +
                t)^2) (1 + 
           E^(-((3 Sqrt[10] delta t)/(-1 + delta)))) GammaRegularized[
          9/2, 0, (0.222222 (4.74342 + (-1. + 1. delta) t)^2)/(-1. + 
              delta)^2])/Sqrt[2 [Pi]], {t, 
       0, (3*Sqrt[10])/(2*(1 - delta))}, AccuracyGoal -> Infinity, 
          Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
      MaxRecursion -> 100] - 
      3*NIntegrate[(E^(-(1/2) (-((3 Sqrt[5/2] delta)/(-1 + delta)) + 
               t)^2) (1 + 
           E^(-((3 Sqrt[10] delta t)/(-1 + delta)))) GammaRegularized[
          9/2, 0, 1/18 ((3 Sqrt[5/2])/(-1 + delta) + t)^2])/
       Sqrt[2 [Pi]], 
          {t, 0, (3*Sqrt[10])/(2*(1 - delta))}, 
      AccuracyGoal -> Infinity, 
          Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
      MaxRecursion -> 100])

Plot[myfunction[delta], {delta, 0, 1}]

It won't finish after serval minutes.
I used AccuracyGoal -> Infinity and so because, when calculating and plotting the result of the integration separately, I get better results.
So, is there any way to speed up this plot (and the calculation itself, as a matter of fact)? What am I missing?

Michael E2 · Accepted Answer

Update: Typo in the original code made it work more easily.

Compute one integral, not two.
Don't use "LocalAdaptive" -- it's usually slower except when it isn't (and sometimes it isn't).
Use MaxRecursion -> 0 in Plot to experiment with possible solutions. (It showed rather quickly there was a discontinuity around 0.8. Suddenly, the plot jumped to zero. Numerics problem? Underflow? Increase working precision?)
EvaluationMonitor confirmed that when MaxRecursion is raised a little, Plot gets bogged down around 0.8.
I guessed where the floating-point numbers came from and substituted exact formulas, so that I could raise WorkingPrecision.
(Update.) The actual problem exceeds the limits of bignums as delta approaches 1 (around delta == 1 = 1*^-3).  It appears this causes the problems with the integration. If we reflect the interval and replace t by Exp[t], the effective support of the integrand becomes a larger proportion of the interval of integration -- and whether or not those are the reasons, the integration is done more easily by NIntegrate until around delta == 1 = 1*^-7.  Not much gain I guess, unless you're interested in delta -> 0.  The plot below seems to work, but that's probably because it does not get so close to the end point delta == 1.

Code:
integrand[t_, delta_] := ((E^(-(1/2) (-((3 Sqrt[5/2] delta)/(-1 + delta)) + 
             t)^2) (1 + 
         E^(-((3 Sqrt[10] delta t)/(-1 + delta)))) GammaRegularized[9/
        2, 0, (2 (3 Sqrt[5/2] + (-1 + 1 delta) t)^2)/(
        9 (-1 + delta)^2)])/Sqrt[2 π]) -
   ((E^(-(1/2) (-((3 Sqrt[5/2] delta)/(-1 + delta)) + t)^2) (1 + 
         E^(-((3 Sqrt[10] delta t)/(-1 + delta)))) GammaRegularized[
        9/2, 0, 1/18 ((3 Sqrt[5/2])/(-1 + delta) + t)^2])/
     Sqrt[2 π]);

myfunction[delta_] := 
  5*(2 + 3 (NIntegrate[
       Exp[t] integrand[(3*Sqrt[10])/(2*(1 - delta)) - Exp[t], delta],
       {t, -Infinity, (3*Sqrt[10])/(2*(1 - delta)) // Log},
       WorkingPrecision -> 32, PrecisionGoal -> 6]));

Plot[myfunction[delta], {delta, 0, 1}(*, MaxRecursion->0*), 
  WorkingPrecision -> 32, 
  PlotRange -> {myfunction[0] - 0.5, All}] // AbsoluteTiming

(Plot gives a precision warning because, even though WP -> 32 is specified, it tests the function by plugging in a machine-precision float.)

Yet another slow plot involving numerical integration

One Answer

Add your own answers!

Ask a Question