Mathematica Asked on August 6, 2021
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?
Update: Typo in the original code made it work more easily.
"LocalAdaptive"
-- it's usually slower except when it isn't (and sometimes it isn't).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
.WorkingPrecision
.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.)
Correct answer by Michael E2 on August 6, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP