NIntegrate does not converge for this integrand but Simpson's Rule does

Question

I am trying to numerically integrate a function that flips sign at a certain location. In the vicinity of this location the function is approximately asymmetric so there are large negative and positive contributions to the integral which presumably cancel. NIntegrate won't converge, but I have written a separate simple Python code using Simpson's Rule which converges fine to the analytic result (the integral can be done analytically). Is there a technique I can use in Mathematica to converge, and/or why does my simple code perform better than NIntegrate? The function is:

func1 = (0.00034677878016334127*(4.028933642832898*^12 - 2.0072104720189415*^6*w)*w^2)/
   (E^(9.489583333333334*^-16*w^2)*(4.028933642832898*^12 + w^2 - 4.014420944037883*^6*Abs[w])^(3/2))

and this is to be integrated as

NIntegrate[func1, {w, 0, 1.377249681416402*^8}]

Answer 1

A numerical integration routine like NIntegrate must sample the function. If the function changes a lot over a very small region, it may happen that the sample process overlooks this irregularity. You can help by specifying a method that samples the function not uniformly, but calculates more samples points where the function changes fast. This does e.g. the method "LocalAdaptive". With this your integral is calculated without problems

func1 = (0.00034677878016334127*(4.028933642832898*^12 - 
      2.0072104720189415*^6*w)*
    w^2)/(E^(9.489583333333334*^-16*w^2)*(4.028933642832898*^12 + 
       w^2 - 4.014420944037883*^6*Abs[w])^(3/2))

NIntegrate[func1, {w, 0, 1.377249681416402*^8}, 
 Method -> "LocalAdaptive"]

(* -5.25938*10^10 *)

Answer 2

You have

func1 = (0.00034677878016334127*(4.028933642832898*^12 - 
      2.0072104720189415*^6*w)*
    w^2)/(E^(9.489583333333334*^-16*w^2)*(4.028933642832898*^12 + 
       w^2 - 4.014420944037883*^6*Abs[w])^(3/2))

If you execute

modf = Rationalize[FullSimplify[func1], 0]

you get

This actually avoids precision problems. The same can be done for the endpoint of integration

1.377249681416402*^8 // Rationalize[#, 0] &

953883129349/6926

You can plot the function to have a look

Plot[modf, {w, 0, 953883129349/6926}, PlotRange -> All]

Why did I point this out?

You can compare the outcome of

NIntegrate[func1, {w, 0, 10^6}, PrecisionGoal -> 100, 
 AccuracyGoal -> 100, WorkingPrecision -> 1000]
NIntegrate[modf, {w, 0, 10^6}, PrecisionGoal -> 100, 
 AccuracyGoal -> 100, WorkingPrecision -> 1000]

both of which give the same number, but the former comes errors, while the latter does not; I know that 10^6 is not the initial endpoint of integration.

Finally, you can use one of the built-in Methods to get the answer:

NIntegrate[modf, {w, 0, 953883129349/6926}, PrecisionGoal -> 100, 
 AccuracyGoal -> 100, WorkingPrecision -> 1000, 
 Method -> "LocalAdaptive"]

As an extra bonus, you can use Alex Trounev's answer to have a function that performs SimpsonIntegral

SetAttributes[SimpsonIntegral, HoldAll]
SimpsonIntegral[f_, x_, xmin_, xmax_, 
   steps_] := (xmax - xmin)/(3 steps) (Sum[
      f /. {x -> xmin + (xmax - xmin)/steps (2*y - 2)}, {y, 1, 
       steps/2}] + 
     4*Sum[f /. {x -> xmin + (xmax - xmin)/steps (2*y - 1)}, {y, 1, 
        steps/2}] + 
     Sum[f /. {x -> xmin + (xmax - xmin)/steps (2*y)}, {y, 1, 
       steps/2}]);

SimpsonIntegral[modf, w, 0, 953883129349/6926, 5000] // N

gives back

-9.68443*10^11

and I think as you increase the steps the result should be more trustworhty - I did not check thoroughly.

Answer 3

Couple more ways:

NIntegrate[func1, {w, 0, 1.377249681416402*^8}, MinRecursion -> 6]
(*  -5.25938*10^10  *)

The critical point (a minimum) of the denominator is where the standard global-adaptive runs into trouble:

cp = NSolveValues[
  D[(E^(9.489583333333334*^-16*w^2)*(4.028933642832898*^12 + w^2 - 
          4.014420944037883*^6*w)^(3/2)), w] == 0 && 
   0 < w < 1.377249681416402*^8, w]
(*  {2.00721*10^6}  *)

NIntegrate[func1, {w, 0, Sequence @@ cp, 1.377249681416402*^8}]
(*  -5.25938*10^10  *)

Alternatively, use the critical points of the function itself:

cp2 = NSolveValues[
  D[func1 /. Abs[w] -> w, w] == 0 && 0 < w < 1.377249681416402*^8, w]
(*  {2.00278*10^6, 2.0117*10^6}  *)

NIntegrate[func1, {w, 0, Sequence @@ cp2, 1.377249681416402*^8}]
(*  -5.25938*10^10  *)