# Yet another slow plot involving numerical integration

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?

• You put high AccuracyGoal, but your numerical constants are just in machine precision. Jul 24 at 20:49
• Should I remove the accuracy goal? Jul 24 at 20:50
• 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 Jul 24 at 20:56
• @yarchik I thought that it was better not to use := when it comes to plotting 'resource-demanding' functions. Jul 24 at 21:41
• NIntegrate uses AccuracyGoal -> Infinity by default. It means that PrecisionGoal determines convergence. Jul 25 at 2:16

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.)

• You have an extra factor of 3 on the second term in your integral. Jul 26 at 15:40
• @bRost03 Good catch! Thanks. Jul 26 at 16:16