Highly oscillating integral

Question

I am trying to integrate a highly oscillatory integral. the integrand is a function of 4 variables:

f1, f2, x2, y

the integrand is

int= (768*E^(8*I*y)*f1*(1 - f1 - f2)*f2*Pi^6*x2*If[1 >= f1 + f2, 1, 0]*
   ((-4*x2*(319 + 239*x2^2) + (533 + 2388*x2^2 + 370*x2^4)*Cos[y] - 
      4*x2*(113 + 178*x2^2)*Cos[2*y] - (139 - 138*x2^2 - 50*x2^4)*Cos[3*y] - 
      4*x2*(-13 + 3*x2^2)*Cos[4*y] - (-29 + 6*x2^2)*Cos[5*y] - 4*x2*Cos[6*y] - 
      3*Cos[7*y])*Log[1 + x2^2 - 2*x2*Cos[y]] + 
    x2*(2*(-101 - 31*x2^2 + 113*Cos[2*y] - 13*Cos[4*y] + Cos[6*y] + 420*Log[x2]) + 
      x2*(154*Cos[y] + 36*x2*Cos[2*y] - 161*Cos[3*y] + 26*x2*Cos[4*y] + 7*Cos[5*y] - 
        20*((126 + 37*x2^2)*Cos[y] + x2*(-74 - 52*Cos[2*y] + 5*x2*Cos[3*y]))*Log[x2]) - 
      (1/x2)*(12*(35 + 120*x2^2 + 18*x2^4 - 20*x2*(7 + 6*x2^2)*Cos[y] + 
         2*x2^2*(45 + 8*x2^2)*Cos[2*y] - 20*x2^3*Cos[3*y] + x2^4*Cos[4*y])*Csc[y]*
        (y*(Log[1 - x2/E^(I*y)] + Log[1 - E^(I*y)*x2]) + I*PolyLog[2, 1 - x2/E^(I*y)] - 
         I*PolyLog[2, 1 - E^(I*y)*x2]))))*Sin[y])/
  ((-1 + E^(2*I*y))^8*(f1 - f1^2 + f2*x2^2 - f2^2*x2^2 - 2*f1*f2*x2*Cos[y])^2)

the integral i want to perform is

NIntegrate[int, {f1, 0, 1}, {f2, 0, 1}, {x2, 0, \[Infinity]}, {y, 
  0, \[Pi]}]

The oscillations occur for the variable y near 0 and infinity. If you plot int (setting f1,f2,x2 to some random variables, they dont matter), you will see that the function is highly oscillating in these ranges. for instance

Plot[int /. f1 -> 1/2 /. f2 -> 1/3 /. x2 -> 1, {y, 3, \[Pi]}]
Plot[int /. f1 -> 1/2 /. f2 -> 1/3 /. x2 -> 1, {y, 0, .1}]

Answer 1

You can raise WorkingPrecision to handle the round-off error in the integrand. A 4D integral can be hard. Here calculating the integrand at high precision seems to take some time. We can see that raising PrecisionGoal adds considerably to the time taken, but note that the first three digits of the PrecisionGoal -> 3 result agreed with the PrecisionGoal -> 4 result. The result is probably trustworthy.

NIntegrate[
  int, {f1, 0, 1}, {f2, 0, 1}, {x2, 0, ∞}, {y, 0, π}, 
  PrecisionGoal -> 3, WorkingPrecision -> 30, 
  Method -> "MultidimensionalRule"] // AbsoluteTiming
(*  {320.151, 383544.798661976140977660944301}  *)

NIntegrate[
  int, {f1, 0, 1}, {f2, 0, 1}, {x2, 0, ∞}, {y, 0, π}, 
  PrecisionGoal -> 4, WorkingPrecision -> 30, 
  Method -> "MultidimensionalRule"] // AbsoluteTiming
(*  {1600.31, 383444.873108674521479252524670}  *)

---

Reply to comments:

@esches "GlobalAdaptive" is a strategy and "MultidimensionalRule" is a rule. Using Method -> rule automatically uses the global adaptive strategy. IIRC, I wanted to override the choice of "LevinRule" in NIntegrate above. (To specify both a strategy and rule with suboptions, use Method -> {strategy, Method -> {rule, ruleoptions}, strategyoptions}.) A discussion of precision/accuracy can be found here.

In general higher WorkingPrecision is needed when the problem has poor conditioning; in your case that comes from the integrand (in particular PolyLog near 0 and π). If you want a more accurate result above, raise PrecisionGoal; however, you will probably have to raise WorkingPrecision, too, since int gets worse at the endpoints and at some point you will need accurate values in those parts of the domain.

Sometimes raising "MaxErrorIncreases" solves the problem. But it's important to realize that you can get that error when round-off error is bad. You can see why from your graph. When you increase the sampling in the wild parts of the graph, the value of the integral and its error estimate will also vary wildly. It's virtually impossible for the integral to converge (and it seems to, it won't converge to the correct value). So before raising "MaxErrorIncreases", you should investigate round-off error. You can do that like I did in the last paragraph of this answer. Often the problem is in a small part of the domain, but luckily you found it. For instance, the following shows a loss of over 75 digits of precision (provided int is nonzero):

Block[{f1, f2, x2, y},
 {f1, f2, x2, y} = SetPrecision[{1/2, 1/3, 1, Pi - 1*^-8}, 100];
 100 - Precision[int]
 ]
(*  75.3115  *)

Of course, that is very close to the boundary of the domain. If you want to check the average over the interior, you can use Table, but for a 4D domain, extensive checking will be time-consuming:

(* scaled sampling points in interval {0, 1} *)
(* (open nonuniform nodes to get nearer the edges of the domain) *)
N[Sin[Subdivide[0, Pi/2, 8][[2 ;; -2]]]^2]
(*  {0.0380602, 0.146447, 0.308658, 0.5, 0.691342, 0.853553, 0.96194}  *)

With[{s = N[Sin[Subdivide[0, Pi/2, 8][[2 ;; -2]]]^2, 100]},
  100 - Precision /@ Flatten@
     DeleteCases[
      Table[int,
       {f1, s}, {f2, s},
       {x2, 1000^s},                       (* scaled to {1, 1000} ~= {1, ∞} *)
       {y, Rescale[s, {0, 1}, {0, Pi}]}],  (* scaled to {0, Pi} *)
      x_ /; x == 0, Infinity]
  ] // ListPlot

Machine precision has almost 16 digits of precision, so the points above the level 16 are useless for estimating the integral and the ones above 10 or 12 limit the accuracy that can be achieved. Raising "MaxErrorIncreases" would be futile in this case.

But suppose round-off error does not seem to be a problem. Then you might want to check whether the integrand is oscillatory or has singularities that might need special handling. Often I raise "MaxErrorIncreases" to 10000 and see things improve, if it takes only a few seconds. I will look at the integral and error estimates that the message provides and see how bad the estimate is. Sometimes it's close to the precision goal and the result is acceptable (you can lower PrecisionGoal to eliminate the message, if you want). Only rarely does raising "MaxErrorIncreases" solve a problem before I've figured out a better way to do the integral. (And, of course, sometimes I give up. There will always be integrals that cannot be done within given time limit.)