Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question

1 Answer 1

up vote 4 down vote accepted

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.

share|improve this answer
    
Nice answer, Michael, +1 –  rasher Mar 15 at 0:41
    
@rasher Thanks. I added an alternate method, which is probably better, and in any case, it confirms the first one. –  Michael E2 Mar 15 at 2:16

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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