# Code that produces plot in V5 doesn't work in later versions

I have problem in plotting Integral function.

I can compute/plot the graph of this integration below in Mathematica 5.0, but it is not possible to plot it in higher Mathematica versions.

My code is:

w = 100
b = 438095.2 (2.5*10^6 + 547826*y^2)
d = Sqrt[-8.8*10^21 *y^2 + 1.92*10^11 *(2.52*10^6 + 547826*y^2)^2]
p = (-b - d)/2
p1 = (-b + d)/2
T = 10^18 *Sqrt[π/2]*NIntegrate[BesselJ[0, y]/Sqrt[(y^2 + w^2)^3]*1/
d*(1 - Erf[(p*10^-15 - 1)/Sqrt[2]])*
Exp[(p*10^-15)^2/2 - p*10^-15*(1 - τ)]*y*w, {y, 0, ∞},
Method -> Oscillatory, MaxRecursion -> 12]
Plot[T, {τ, 0, 6*10^3}]

The plot is generated by Mathematica 5.0:

Error in later versions is NIntegrate::nconv

• Post your code rather than a picture so that the code can be copied and pasted into a Mathematica notebook. – Bob Hanlon Jul 3 at 5:35
• why is $p_1$ defined but not used anywhere? Add you also do not a value for $\tau$ defined anywhere but using it inside a numerical integration. How did this work in V 5? And there is no need to use subscripts, makes the code hard to read. Simply p2write instead Subscript[p, 2] – Nasser Jul 3 at 8:28
• p1 is defined for another integral, I didn't mention it. Yes it works in V5, because at last, in Plot part all the variables are integrated. Can you suggest me another way of writing for this? I am a new user of Mathematica don't pay attention at subscripts :) – Ismatov Tolib Jul 3 at 9:12
• Yes it works in V5 how could numerical integral work when there is undefined symbol (in your case $\tau$) inside it, other than the integration variable? This makes little sense to me as numerical integration needs numerical values for all parameters. Could you show screen shot, from clean kernel showing it works? This will help clarify things. – Nasser Jul 3 at 9:46
• @Nasser you can realised it like this Plot[10^18 Sqrt[[Pi]/2] NIntegrate[ BesselJ[0, y]/Sqrt[(y^2 + w^2)^3]*1/ d*(1 - Erf[(p*10^-15 - 1)/Sqrt[2]])* Exp[(p*10^-15)^2/2 - p*10^-15*(1 - [Tau])]*y*w, {y, 0, [Infinity]}], {[Tau], 0, 6*10^3}] – Ismatov Tolib Jul 3 at 9:59

## The integral diverges, please see Michael's answer for more detal. This is just an answer reproducing result of v5.

Unlike higher versions, v5 is using very few points for the generation of graphic:

plot = Plot[T, {τ, 0, 6*10^3}]; // AbsoluteTiming
(* {32.1406250 Second, Null} *)

plot[[1, 1, 1, 1]]
% // Length

To obtain the result in a reasonable time in higher versions, just limit the points for plotting:

Quiet@Plot[T, {τ, 0, 6*10^3}, MaxRecursion -> 0] // AbsoluteTiming

Still a bit slower than v5, but acceptable in my view.

## Update: A More General Solution

The solution above works at least in v12.0.1, v12.1.1 and v11.3, but in v8.0.4 and v9.0.1 samples like

T /. τ -> 1 // AbsoluteTiming

returns unevaluated (to be precise, NIntegrate[…] therein returns unevaluated) after nconv warning generates. If one still needs the result in v5, a possible approach is to implement the "ExtrapolatingOscillatory" method ourselves as I've done here:

Clear[int, separateint]
zero[i_] = Piecewise[{{BesselJZero[0, i], i > 0}}];
separateint[f_, t_, i_?NumericQ, prec_] :=
NIntegrate[BesselJ[0, y] f[y, t], {y, zero@i, zero[i + 1]}, WorkingPrecision -> prec,
MaxRecursion -> 40];
int[f_, t_?NumericQ, prec_ : MachinePrecision] :=
NSum[separateint[f, t, i, prec], {i, 0, Infinity}, Method -> "AlternatingSigns",
WorkingPrecision -> prec];

The usage of int is as follows:

Clear@func;
func[y_, τ_] =
1/Sqrt[(y^2 + w^2)^3]*1/d*(1 - Erf[(p*10^-15 - 1)/Sqrt[2]])*
Exp[(p*10^-15)^2/2 - p*10^-15*(1 - τ)]*y*w;

T = 10^18*Sqrt[π/2]*int[func, τ];

T /. τ -> 1
(* 84.0182 *)

Plot[T, {τ, 0, 6*10^3}, MaxRecursion -> 0] // AbsoluteTiming

The resulting graphic is the same as shown above, but the solution also works in v8.0.4 and v9.0.1.

This solution also works for the Tc and Td in your new question, on which NIntegrate of v12 again returns unevaluated.

• XZCZD thanks for advice. Can you give me your contact info if I send you another case of the problem. I'm really confused. – Ismatov Tolib Jul 4 at 10:56
• @IsmatovTolib I'm sorry, but it's my stupid insistence that I only discuss Mathematica in public. Just ask your question in this site, if I find it interesting and know the solution, I'll answer. (Do remember to make the question on-topic. ) – xzczd Jul 4 at 11:38
• I asked another question about my whole problem. The same problembut a little bit complex function. Can you see please. Thanks for your effort. – Ismatov Tolib Jul 4 at 19:30
• XZCZD by using which math version you produced the graph above? – Ismatov Tolib Jul 6 at 19:01
• I should add the note that yes, the old Oscillatory method in pre-version 6, and the equivalent "ExtrapolatingOscillatory" method in newer versions, are able to evaluate divergent oscillatory integrals in the Abel sense, so depending on what you're doing, you have to be careful. Here's one of my favorite examples: NIntegrate[Log[x] Sin[x], {x, 0, ∞}, Method -> "ExtrapolatingOscillatory", WorkingPrecision -> 25] and compare with Limit[LaplaceTransform[Log[x] Sin[x], x, t], t -> 0]. – J. M.'s technical difficulties Jul 8 at 12:40

## This is not a bug. It's an improvement. The integral is divergent.

The V5 Oscillatory method is defunct. NIntegrate chooses the "ExtrapolatingOscillatory" method (which is the method it chooses for this integral if Method -> Automatic). This method checks convergence, and the amplitude of the oscillations goes to infinity. Therefore you cannot apply this method. Whether or not one can make sense of the integral or series used in the extrapolating oscillatory method by applying methods for divergent series, I did not look further into. The nice shape of the graph suggests something like Euler's approach to $$\sum (-1)^n n!x^{n+1}$$ might be possible.

Analytic divergence. Let's define the integrand two ways, the original machine-precision one as integrand0 and the rationalized, exact one as integrand.

integrand = Rationalize[
Rationalize[
integrand0 =
BesselJ[0, y]/Sqrt[(y^2 + w^2)^3]*1/
d*(Erfc[(p*10^-15 - 1)/Sqrt[2]])*
Exp[(p*10^-15)^2/2 - p*10^-15*(1 - τ)]*y*w
],
0];

The aymptotic expansion of BesselJ[0, y] has a leading term proportional to 1/Sqrt[y]:

Normal@Series[BesselJ[0, y], {y, Infinity, 0}]
(*  (Sqrt[2/π] Cos[π/4 - y])/Sqrt[y]  *)

To check convergence, the "ExtrapolatingOscillatory" strategy replaces BesselJ[0, y] by 1/Sqrt[y] and takes the limit at infinity. This results in ComplexInfinity, even if τ does not have a numeric value:

Limit[integrand0 /. _BesselJ -> 1/Sqrt[y], y -> Infinity]
(*  ComplexInfinity  *)
Limit[integrand /. _BesselJ -> 1/Sqrt[y], y -> Infinity,
Assumptions -> τ > 0]
(*  Infinity  *)

(The exact limit is difficult for Limit to evaluate without some helpful assumption about τ.) We can see the actual limit use in the following trace:

ClearSystemCache[] (* AsymptoticsClassicLimit caches some results *)
Trace[
10^18*Sqrt[π/2]*
NIntegrate[
BesselJ[0, y]/Sqrt[(y^2 + w^2)^3]*1/
d*(Erfc[(p*10^-15 - 1)/Sqrt[2]])*
Exp[(p*10^-15)^2/2 - p*10^-15*(1 - τ)]*y*w,
{y, 0, ∞},
Method -> Oscillatory,(* defunct method *)
MaxRecursion -> 12],
_AsymptoticsClassicLimit,
TraceInternal -> True, TraceForward -> True] /.
x_Times /; Simplify[x - integrand0 /. _BesselJ -> 1/Sqrt[y]] == 0 ->
"integrand0"

Numerical divgence 1. We've seen that things in Mathematica can go wrong, so we should look for further verification.

Here is an evaluation of the integrand at increasing powers of 2.

Block[{τ = 1000},
integrand /. y -> N[2^Range@14, 16] // RealExponent
]
(*
{-17.3277, -18.4065, -23.9926, -44.2635, -124.654, -444.958,
-1723.93, -6791.65, -26461., -95463.2, -216725., 1.77412*10^6,
4.93517*10^7, 8.7349*10^8}
*)

Note these are the exponents (base $$10$$). Initially, the integrand seems to converge very rapidly to zero (down to $$10^{-216725}$$). This is why NIntegrate misses the divergence in Plot: Plot applies N to NIntegrate when it first returns unevaluated, and N[NIntegate[..]] somehow forces the fall-back method of the "ExtrapolatingOscillatory" strategy, which by default is "GlobalAdaptive" with "GaussKronrodRule". The "GaussKronrodRule" fails to detect divergence for smaller values of τ because it samples only in the apparently-convergence part of the domain.

However, it can be seen that the amplitude takes off, reaching $$10^{10^8}$$ at $$y = 16384$$.

Numerical divergence 2. As another confirmation, we can modify @xzczd's implementation of the extrapolating oscillatory method. We add NSumTerms -> 100:

int[f_, t_?NumericQ, prec_ : MachinePrecision] :=
NSum[separateint[f, t, i, prec], {i, 0, Infinity}, Method -> "AlternatingSigns",
NSumTerms -> 100
WorkingPrecision -> prec];

Then we get a much different result:

T /. τ -> 1
(*  2.09192*10^153  *)

My original answer is only partly right, but using a proper definition turned out to be a minor side issue. Here's an easier way to get plot of the "value" of the integral. Erelyi's approach to Euler's sum for $$\sum (-1)^n n!x^{n+1}$$ in Asymptotic Expansions is to observe that for small $$x$$, "the terms of the series at first decrease quite rapidly, and an approximate numerical value of [the sum] may be computed." Likewise we can truncate the interval of integration near the minimum amplitude and cut off the divergent part. No analysis has been made of whether this value can be considered equivalent to the integral, whether in the sense of Euler's paper or in some other way.

ClearAll[T];
T[τ0_?NumericQ] :=
Block[{τ = SetPrecision[τ0, Infinity]},
10^18*Sqrt[π/2]*NIntegrate[integrand,
Evaluate@{y, 0,
Max[15, First@
FindArgMin[
integrand /. _BesselJ -> 1/Sqrt[y], {y, 50 + τ, 0,
10000}, WorkingPrecision -> 25]
]}, MaxRecursion -> 12]
];

Plot[T[τ1], {τ1, 0, 6*10^3},
MaxRecursion -> 1] // AbsoluteTiming

The problem (for V12.1.1) is that T is not properly defined. Making it a function of τ gets rid of all error/warning messages:

ClearAll[T];
T[τ_?NumericQ] :=
10^18*Sqrt[π/2]*
NIntegrate[
BesselJ[0, y]/Sqrt[(y^2 + w^2)^3]*1/
d*(Erfc[(p*10^-15 - 1)/Sqrt[2]])*
Exp[(p*10^-15)^2/2 - p*10^-15*(1 - τ)]*y*w,
{y, 0, ∞}, MaxRecursion -> 12];
Plot[T[τ], {τ, 0, 6*10^3},
MaxRecursion -> 1] // AbsoluteTiming
• The way to define T isn't too big a problem here. I've tested both, the timing is similar. (One needs to test with a fresh kernel, the timing reduces significantly after the first run. ) I think the real problem is a bug of NIntegrate, see the update to my answer. BTW the warning nconv is still there using this method. Tested on v12.1.1, Win10. – xzczd Jul 7 at 3:23
• @xzczd What's the timing got to do with it? I thought the OP's problem was producing the graph. I get no error messages. (You can ClearSystemCache[] instead of restarting the kernel.) – Michael E2 Jul 7 at 3:27
• I just mean τ_?NumericQ actually helps little, the original T also works in v12.1.1, and the performance is almost the same. – xzczd Jul 7 at 3:33
• @xzczd I made a mistake, but I think we both did. I don't think it's a bug. – Michael E2 Jul 7 at 17:50
• Oh… I did suspect the integral diverges, but the moment I see the result of self-implementation of ExtrapolatingOscillatory I forgot everything… I've edited my answer accordingly, thx for pointing out, but…"This is why NIntegrate misses the divergence in Plot: Plot applies N to NIntegrate when it first returns unevaluated, and N[NIntegate[..]] somehow forces the fall-back method " This isn't quite right, because NIntegrate[integrand0 /. τ -> 1, {y, 0, ∞}] evaluates. Or you actually mean NIntegrate[integrand, …? – xzczd Jul 8 at 2:40