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

Question

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

Can anybody please help? Any suggestions appreciated.

Answer 1

The integral diverges, please see Michael's answer for more detail. 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.

Answer 2

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[] (* Asymptotics`ClassicLimit 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],
  _Asymptotics`ClassicLimit,
  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

Original answer

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