3
$\begingroup$

I am trying to integrate a function, and the error I get is greater than the result.

So I need to calculate

NIntegrate[osc[n,L,u,x]*f1[u, x]/f2[u, x], {u, 0, 1}, {x, 1, inf}]

where

f1[u_,x_] := u(1-u) ( 1-Cos[x]-x Sin[x] )

f2[u_,x_] := (1-u)^2 x^2 + 2u(1-u) x Sin[x] + 2( 1-Cos[x] ) u^2

and 

osc[n_,L_,u_,x_]:= Exp[-I u x(n+L)+2I n ArcTan[(u(1-Cos[x]))/(u Sin[x]+x(1-u))]]

Here, n and L are parameters: n is an integer (1, 2, 3, …), and L is real > 0

The best I got was by using Method -> "LevinRule":

For example, by using

nn=1; LL = 1;

I get

int[n_,L_,u_,x_] := osc[n,L,u,x]*f1[u, x]/f2[u, x];
inf = Infinity;
NIntegrate[ int[nn, LL, u, x], {u, 0, 1}, {x, 1, inf}, Method -> 
   {"LevinRule", "Kernel" -> osc[nn, LL, u, x]}]

(* -0.0741726 + 0.0851585 I *)

with an error of 0.18027814

Note that the error is greater than the result.

|improve this question
$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$ – bbgodfrey Feb 20 '15 at 20:49
  • $\begingroup$ The option WorkingPrecision might help. $\endgroup$ – Jinxed Feb 20 '15 at 22:25
  • 1
    $\begingroup$ What are int[] and inf? $\endgroup$ – Dr. belisarius Feb 20 '15 at 22:39
  • $\begingroup$ Sorry, int[n,L,u,x]=osc[n,L,u,x]*f1[u, x]/f2[u, x], and inf= infinity $\endgroup$ – oaroldan Feb 21 '15 at 14:48
  • $\begingroup$ @Jinxed Thanks Jinxed: increasing WorkingPrecision just give me more digit numbers but doesn’t change the results while spending more time. $\endgroup$ – oaroldan Feb 21 '15 at 19:23
2
$\begingroup$

It seems to me that infinitely many, semi periodic singularities makes for a very difficult integral. If you could find a substitution that would tame all the singularities at once, you could probably succeed in truly taming this integral.

I agree with the OP that each singularity is integrable, but one has also to wonder whether the sum of the infinitely many singularities exists. My feeling is that it does. Even though the magnitude of the integrand looks like it falls off as $1/x$, the oscillatory nature of the function does not rule out convergence, and the numerical calculations I get stay around 0.08(6) + 0.08(1) I, not too far from the OP's result.

My thought was that since the singularities occur as $u \rightarrow 1^-$, $x \rightarrow 2 \pi n^-$, $n$ an integer, we might divide the infinite strip {u, 0, 1}, {x, 0, Infinity} into infinitely many rectangles {u, 0, 1}, {x, 2 Pi n, 2 Pi + 2 Pi n} and sum the integrals. (One has to subtract the integral over {u, 0, 1}, {x, 0, 1} to get the OP's integral, of course.) A singularity occurs only in one corner of each rectangle. This suggests that using Duffy's coordinates might be useful. The substitution transforms a singularity at a corner to one along an edge. If the singularity is of the right type, it can be quite effective. In this case, it helps, but the singularity persists to be a nuisance. However, the IMT singularity handler helps, at least to the extent that with a higher working precision, the integral can be computed with the desired precision.

Summing the integrals over the rectangles can be done with NSum. This means approximating the sum using approximate terms. I'm not real sure what the compounding of the error implies. Given that all terms were computed with a relative precision of MachinePrecision/2 digits and the real and imaginary parts of the terms all have the same sign (except the imaginary part of the first term), the sum should be nearly as good as NSum usually calculates, I suppose.

Another strategy I used was to isolate the singularity in the corner. Away from the corner NIntegrate can accurately compute the integral with the "LevinRule" method. Luckily the singularity can be predicted to lie in the rectangle {u, 1 - 1/Max[n, 1]/4, 1}, {x, 3 Pi/2 + 2 Pi n, 2 Pi + 2 Pi n}. We just have to figure out how to integrate this singularity.

For some reason Method -> {"DuffyCoordinates", "Corners" -> {{1, 1}}} was insufficient to compute the integrals with the desired precision. But the NIntegrate tutorial linked above includes code for manually performing the substitution. We can then call NIntegrate specifying the singularity handler "SingularityHandler" -> "IMT", which finally seems to conquer the problem.

So finally, here is the code that implements the above strategy. The function DuffyCoordinates is given in this section of the NItegrate tutorial. I would include it, but I'm not sure about copyright restrictions. Warning: It takes a long time, about a minute or so per rectangle; hence the results are memoized. Print statements are included for the impatient. The singular integral always gives a singularity warning NIntegrate::slwcon, so I supressed it with Quiet.

ClearAll[intpart7, f1, f2, osc, int];

(* OP's definitions *)
f1[u_, x_] := u (1 - u) (1 - Cos[x] - x Sin[x])
f2[u_, x_] := (1 - u)^2 x^2 + 2 u (1 - u) x Sin[x] + 2 (1 - Cos[x]) u^2
osc[n_, L_, u_, x_] := 
 Exp[-I u x (n + L) + 2 I n ArcTan[(u (1 - Cos[x]))/(u Sin[x] + x (1 - u))]]
nn = 1; LL = 1;
int[u_, x_] := osc[nn, LL, u, x]*f1[u, x]/f2[u, x]

(* integrates over the n-th rectangle *)
intpart7[n_Integer] := Block[{x, u},
   Print[n];
   intpart7[n] = (* memoize for reuse *)
    NIntegrate[            (* nonsingular integral *)
      Boole[0 <= u <= 1 - 0.25/Max[n, 1] || 2 Pi n <= x <= 2 Pi n + 3 Pi/2]*int[u, x],
      {u, 0, 1}, {x, 2 Pi n, 2 Pi + 2 Pi n},
      MaxRecursion -> 25, 
      Method -> {"GlobalAdaptive", 
        Method -> {"LevinRule", "Kernel" -> osc[nn, LL, u, x]}}] +
     Quiet[Print["."];
      NIntegrate @@ Join[  (* singular part *)
        DuffyCoordinates[
         int[u, x], {u, 1 - 1/Max[n, 1]/4, 1},
         {x, 3 Pi/2 + 2 Pi n, 2 Pi + 2 Pi n}, {1, 1}],
        {PrecisionGoal -> 8, WorkingPrecision -> 20, 
         MaxRecursion -> 25, 
         Method -> {"GlobalAdaptive", "SingularityHandler" -> "IMT"}}
        ], NIntegrate::slwcon]
   ];

Grand result

intzero = NIntegrate[int[u, x], {u, 0, 1}, {x, 0, 1}]
(sum = NSum[intpart7[n], {n, 0, Infinity}]) // AbsoluteTiming
sum - intzero
(*
  -0.0751088 + 0.0187646 I                (* integral 0 < x < 1 *)
  {1891.108293, -0.161177 + 0.0992955 I}  (* {time, sum} *)

  -0.0860682 + 0.0805309 I                (* final integral *)
*)

No error estimate warnings!

Remark: The computed terms may be inspected with ? intpart7 or DownValues[intpart7]. The indices n computed are

Cases[DownValues[intpart7], HoldPattern[intpart7[n_Integer]] :> n, Infinity]
(*
  {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
   19, 20, 21, 22, 23, 24, 25, 26, 30, 93305, 93306, 93307}
*)

Accuracy of results

We can see that using more terms in the estimation by NSum tends to increase the magnitude of the result, but by less and less. I let this run, while I was doing other stuff. It probably takes a total of two hours, including the time above. NSum does not seem to estimate the tail of the sum very accurately. It's possible that it diverges very slowly, but the OP seems to think it converges. I have not thought very much about it.

seq = Table[NSum[intpart7[n], {n, 0, Infinity}, NSumTerms -> n] - intzero, {n, 15, 100, 10}]
(*
{-0.0860682 + 0.0805309 I,
 -0.0861501 + 0.0805336 I,
 -0.0860646 + 0.0805342 I,
 -0.0861001 + 0.0805403 I,
 -0.0861384 + 0.080535 I,
 -0.0861612 + 0.0805338 I,
 -0.0861782 + 0.080533 I,
 -0.0861947 + 0.080534 I,
 -0.0862018 + 0.0805343 I}
*)

ListLinePlot[Abs@Differences@seq, DataRange -> {15, 95}]

Mathematica graphics

Visualization of the Duffy transformation

The plots below show the real and imaginary parts of the integrand in the singular corner of the domain for n = 10 before and after the substitution of Duffy's coordinates. One can see that the singularity is moved to the edge as required by the "IMT" singularity handler.

Block[{n = 10, style, fn},
 {style[Re], style[Im]} = 
  Take["DefaultPlotStyle" /. (Method /. Charting`ResolvePlotTheme["Default", Plot3D]), 2];
 fn["int"] = {##} &;
 fn["Duffy"] = {First@DuffyCoordinates[##], {u, 0, 1/n}, {x, 1 - 2/n, 1}} &;
 GraphicsGrid@
  Table[Plot3D[part[#1], #2, #3, PlotRange -> 5, MaxRecursion -> 3, 
      PlotPoints -> 35, PlotStyle -> style[part], 
      AxesLabel -> Automatic, PlotLabel -> part[f]] & @@ 
    fn[f] @@ {int[u, x], {u, 1 - 1/Max[n, 1]/4, 1}, {x, 3 Pi/2 + 2 Pi n, 2 Pi + 2 Pi n}, {1, 1}},
   {part, {Re, Im}}, {f, {"int", "Duffy"}}
   ]
 ]

Mathematica graphics

|improve this answer
$\endgroup$
  • $\begingroup$ Michael Thanks a lot! Your help and commitment really motives me to also help other people. I have tested your program and it works perfectly. I will use it to compute a set of integrals that have the same behavior. When finished I'd like to let you now (if you don't bother). $\endgroup$ – oaroldan Feb 27 '15 at 3:12
  • $\begingroup$ @OmarAlbertoRoldanGarcia You're welcome. Best of luck with the rest of your project. Feel free to let me know how it goes. $\endgroup$ – Michael E2 Feb 27 '15 at 4:26
1
$\begingroup$

If you are happy with the findings so far, and just need help on dynamically creating the exclusions, you might want to proceed with the following index-generation code:

DeleteDuplicates@{#1,#2,Sequence@@Range[Ceiling[#2,2\[Pi]],Floor[#3,2\[Pi]],2\[Pi]],#3}&

which can be called with the variable and its lower and upper bounds like so, e.g.:

DeleteDuplicates@{#1,#2,Sequence@@Range[Ceiling[#2,2\[Pi]],Floor[#3,2\[Pi]],2\[Pi]],#3}&[x,1,10]

So, your NIntegrate becomes:

NIntegrate[int[nn, LL, u, x], {u, 0, 1},
  DeleteDuplicates@{#1,#2,
  Sequence@@Range[Ceiling[#2,2\[Pi]],Floor[#3,2\[Pi]],2\[Pi]],#3}&[x,xmin,xmax],
  Method->{"LevinRule", "Kernel"->osc[nn,LL,u,x]}]
|improve this answer
$\endgroup$
  • $\begingroup$ Thanks Jinxed, your tips have been very useful! I took me a little long to answer because I'm on a trip. The errors have indeed diminished. $\endgroup$ – oaroldan Feb 24 '15 at 22:52
  • $\begingroup$ @OmarAlbertoRoldanGarcia: You are very welcome! $\endgroup$ – Jinxed Feb 25 '15 at 14:51

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.