2
$\begingroup$

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?

$\endgroup$
9
  • 1
    $\begingroup$ You put high AccuracyGoal, but your numerical constants are just in machine precision. $\endgroup$
    – yarchik
    Jul 24 at 20:49
  • $\begingroup$ Should I remove the accuracy goal? $\endgroup$
    – Vicent
    Jul 24 at 20:50
  • 1
    $\begingroup$ I would say so, but it does not matter much. In any case, I can do plot in 6 seconds in the range (0,0.95) and 7 seconds in (0,1). There is a problem at 1. I have MA11 $\endgroup$
    – yarchik
    Jul 24 at 20:56
  • $\begingroup$ @yarchik I thought that it was better not to use := when it comes to plotting 'resource-demanding' functions. $\endgroup$
    – Vicent
    Jul 24 at 21:41
  • 1
    $\begingroup$ NIntegrate uses AccuracyGoal -> Infinity by default. It means that PrecisionGoal determines convergence. $\endgroup$
    – Michael E2
    Jul 25 at 2:16
3
$\begingroup$

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


  1. Compute one integral, not two.
  2. Don't use "LocalAdaptive" -- it's usually slower except when it isn't (and sometimes it isn't).
  3. 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?)
  4. EvaluationMonitor confirmed that when MaxRecursion is raised a little, Plot gets bogged down around 0.8.
  5. I guessed where the floating-point numbers came from and substituted exact formulas, so that I could raise WorkingPrecision.
  6. (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.)

$\endgroup$
2
  • $\begingroup$ You have an extra factor of 3 on the second term in your integral. $\endgroup$
    – bRost03
    Jul 26 at 15:40
  • $\begingroup$ @bRost03 Good catch! Thanks. $\endgroup$
    – Michael E2
    Jul 26 at 16:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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