2
$\begingroup$

I am doing:

  NIntegrate[Sin[Exp[(x^4)]], {x, 2, Infinity}, PrecisionGoal -> 12]

It prints out a host of warnings, but also shows the results as:

$$-0.0113737$$

This is a highly oscillatory integral and weirdness happens to the numerical routines as they are doing their thing. For example:

  • integration with LevinRule failed.
  • converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small
  • "NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {4.531228269745986`}"
  • ...

You can see a plot as:

 Plot[Sin[Exp[(x^4)]], {x, 2, 5}]

As an aside, the same input (exactly) using WA, produces:

$$-0.00644455$$

Are those warnings/errors telling me not to trust the result as it is suspect?

$\endgroup$

1 Answer 1

6
$\begingroup$

The Levin rule fails because the argument Exp[x^4] cause an overflow at one of the sample points:

Exp[(x^4)] /. x -> 5.734404367911653`*^7

General::ovfl: Overflow occurred in computation. >>

Overflow[]

First solution

So try a substitution:

Sin@Exp[(x^4)] Dt[x] /. x -> (Log[u])^(1/4)
(*
   (Dt[u] Sin[u])/(4 u Log[u]^(3/4))
*)

Then:

NIntegrate[Sin[u]/(4 u Log[u]^(3/4)), {u, E^16, Infinity}, 
 Method -> "LevinRule", MaxRecursion -> 30, PrecisionGoal -> 12, 
 AccuracyGoal -> 12, WorkingPrecision -> 20]
(*
  -3.0179581945309077605*10^-9
*)

Increasing the AccuracyGoal, PrecisionGoal, and WorkingPrecision shows the first six digits are stable.


Another approach

Looking that the integrand, Sin[u]/(4 u Log[u]^(3/4)), one can see the "amplitude" 1/(4 u Log[u]^(3/4)) is monotonic decreasing. We can treat the integral as an alternating series of areas above and below the axis and bound the terms by a sine function of amplitude 1/(4 u0 Log[u0]^(3/4)), where u0 is the left endpoint when the graph is above the axis and the right endpoint when it is below. The integral of one cycle is less than the difference in the areas. Further, areas adjacent at the point where the graph crosses from below to above are congruent and cancel (see the green shaded regions below). Thus we have a telescoping series and the integral from the yellow region (at any u0) to infinity is bounded by the area of the yellow region, which is 2/(4 u0 Log[u0]^(3/4)).

Mathematica graphics

So if we integrate from x = 2 to x = (Log[10^50 Pi])^(1/4), that is u = 10^50 π, the difference with integral out to infinity is less than 10^-52:

2./(4 u Log[u]^(3/4)) /. u -> 10^50 π
(*
   {4.49478*10^-53}
*)

So we have

NIntegrate[Sin[Exp[(x^4)]], {x, 2, (Log[10^50 Pi])^(1/4)}, 
  PrecisionGoal -> 20, MaxRecursion -> 30, WorkingPrecision -> 30]
(*
   -3.01795244987123683885173256511*10^-9
*)

In this case increasing the precision settings shows that all 30 digits are stable.

$\endgroup$
2
  • $\begingroup$ Nice answer, Michael, +1 $\endgroup$
    – ciao
    Mar 15, 2014 at 0:41
  • 1
    $\begingroup$ @rasher Thanks. I added an alternate method, which is probably better, and in any case, it confirms the first one. $\endgroup$
    – Michael E2
    Mar 15, 2014 at 2:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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