9
$\begingroup$

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}]
$\endgroup$

3 Answers 3

12
$\begingroup$

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 *)
$\endgroup$
1
  • $\begingroup$ (+1) beat me to it. Nicely done!!! $\endgroup$
    – bmf
    Mar 19 at 20:16
10
$\begingroup$

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

opfunction

If you execute

modf = Rationalize[FullSimplify[func1], 0]

you get

modf

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]

plotf

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.

$\endgroup$
6
$\begingroup$

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  *)
$\endgroup$
4
  • $\begingroup$ The critical points approach...such a beauty $\endgroup$
    – bmf
    Mar 20 at 5:05
  • $\begingroup$ @bmf Thanks. :) $\endgroup$
    – Michael E2
    Mar 20 at 5:43
  • 3
    $\begingroup$ Here's another good approach. Since the function is analytic on the interval of integration (Abs[w] is a red herring), use a high-order Clenshaw-Curtis or Gauss-Kronrod rule, which are rapidly converging on analytic functions: NIntegrate[func1, {w, 0, 1.377249681416402*^8}, Method -> {"ClenshawCurtisRule", "Points" -> 65}]. $\endgroup$
    – Michael E2
    Mar 20 at 5:49
  • $\begingroup$ That's f... IMPRESSIVE. Thanks for this. I will keep it mind for future purposes :-) $\endgroup$
    – bmf
    Mar 20 at 6:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.