4
$\begingroup$

I am trying to numerically estimate a complex 3D integral. Here is its version considerably reduced for the purposes of the present discussion:

NIntegrate[(x + ξ)^4/(
 Sqrt[ξ] (y^2 + (x + ξ)^2)^3), {x, -100, 100}, {y, 0, 
  100}, {ξ, 0, 100}]

The integral has an integrable singularity. One can make it sure by passing to spherical coordinates. I believe that it should be possible to calculate it. However, it is this singularity that makes the integral difficult to estimate numerically. Regularization like 

...{y, 0.001,100}, {ξ, 0.001, 100}...

does not help.

I tried most (though yet not all) of methods and strategies for the singular integrals listed in the tutorial/NIntegrateIntegrationStrategies. This was so far unsuccessful.

One of the recommendations I received in the message reporting the divergence was to increase "MaxErrorIncreases". However, I failed to find any documentation showing how to apply this option. Therefore,

My first question: Do you know, how to use "MaxErrorIncreases"?

The second question: Do you have any idea on calculating an integral of this type?

Edit: To address the question of @Akku14 Yes, I managed to integrate it over ksi exactly: 

intt = Integrate[
  1/Sqrt[ξ]* ((x + ξ)^4)/(y^2 + (x + ξ)^2)^3, {ξ, 
   0, ∞}, Assumptions -> {y > 0, x ∈ Reals}]

which yields 

(* 1/64 π (10 (1/(x - I y)^(3/2) + 1/(x + I y)^(3/2)) - (
   12 I (1/Sqrt[x - I y] - 1/Sqrt[x + I y]))/y + 
   3 I (1/(x - I y)^(5/2) - 1/(x + I y)^(5/2)) y)  *)

On the one hand, this result seems suspicious. I am not quite sure that it is right. On the other hand, this does not really help. The integral

NIntegrate[intt, {x, -100, 100}, {y, 0, 100}]

still does not converge, at least, by itself.

Improve this question
$\endgroup$
13
  • $\begingroup$ Is actual integral in infinite limits? $\endgroup$
    – yarchik
    Mar 2, 2020 at 9:45
  • $\begingroup$ @yarchik No. The actual integrand has also a step-like function as a factor. The latter is zero beyond +-100 along x, +100 along y and ksi. This makes the limits finite. I omitted this factor, since the difficulty is related to the singularity, rather than to this factor, and it will make the question more obscure. On the other hand, like this, the limits given above are to extent arbitrary. One can instead +-100 take +-10 to try. $\endgroup$
    – Alexei Boulbitch
    Mar 2, 2020 at 9:50
  • $\begingroup$ The shown example is analytically integrable with respect to at least one variable. May be your real problem is also integrable with Integrate in x,y or xi. This makes things much easier. $\endgroup$
    – Akku14
    Mar 2, 2020 at 10:35
  • 1
    $\begingroup$ You could also play with the singularity handler. Usage example: Method -> {"GlobalAdaptive", "SingularityHandler" -> "DuffyCoordinates", "SingularityDepth" -> 20, "MaxErrorIncreases" -> 20000}. See reference.wolfram.com/language/tutorial/… for documentation of these options. Good luck. Nothing really worked for me yet. $\endgroup$
    – Michael E2
    Mar 2, 2020 at 13:31
  • 1
    $\begingroup$ For me the principal problem is calculating the integral, which I've failed to do, almost as bad as entering MaxErrorIncreases into the documentation search bar. My goodness, not a single relevant hit until the 17th (NIntegrate), which is the only relevant hit in all 33 hits. The tutorial I linked does not even come up in the whole list. Your first question should have been "easily found in the documentation," but it's not. $\endgroup$
    – Michael E2
    Mar 2, 2020 at 13:59

2 Answers 2

Reset to default
1
$\begingroup$

At present I developed a workaround for this problem. It is lengthily, and, generally, I do not like it very much. However, it seems to work. I will be grateful for all sorts of criticisms and suggestions.

So, let us first integrate over x and y, and find out this integral as the function of ξ. Note that I make a regularization of the denominator by the machine number 10^-16. 

    iter = Join[Table[10^-i, {i, 0, 16}] // N, Table[i, {i, 2, 90, 5}], 
    Table[i, {i, 91, 110, 0.1}], Table[i, {i, 115, 200, 5}]] // Sort;

Here iter is a non-homogeneous iterator.

Below int is the list of values of pairs {ξ, integralOverXandY} 

int = Table[{ξ, 
        NIntegrate[(x + ξ)^4/( (y^2 + (x + ξ)^2)^3 + 
          10^-16), {x, -100, 100}, {y, 0, 100}, 
         Method -> {"LocalAdaptive", 
           Method -> {"ClenshawCurtisRule", "Points" -> 10}}, 
         AccuracyGoal -> 3, PrecisionGoal -> 5]}, {ξ, iter}];

In the following I interpolate the results:

f = Interpolation[int, InterpolationOrder -> 1];
Show[{
  ListPlot[int, PlotRange -> All, 
   AxesLabel -> {Style["ξ", 16], Style["int", 16]}],
  Plot[f[ξ], {ξ, 0.0001, 300}, PlotStyle -> Red]
  }]

It looks as follows:

enter image description here

Now it can be integrated over ξ:

NIntegrate[f[ξ]/Sqrt[ξ], {ξ, 10^-16, 300}]

(*  261.063   *)

I also tried to calculate the integral over x and y with the options:

Method -> {"AdaptiveMonteCarlo", "MaxPoints" -> 1000000}, 
 AccuracyGoal -> 2, PrecisionGoal -> 4]

and 

Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule"}, 
 AccuracyGoal -> 2, PrecisionGoal -> 4]

This returns error messages but makes the calculation. The plot appears considerably worse, but the result of integration is astonishingly the same.

Improve this answer
$\endgroup$
2
  • $\begingroup$ With the help of Rubi I could integrate first over y and then over x analytically, The result has a removable singularity, where you have a jump. But the function looks quite different. This 1D function can then be integrated numerically. $\endgroup$
    – yarchik
    Mar 4, 2020 at 21:46
  • $\begingroup$ @yarchik Thank you. It is a very valuable observation. I will try it up. My preliminary trial shows that the result of the Rubi integration over x and y does not coincide with my numeric result, neither does it with my expectation (the latter, of course, may be wrong). I will look at it precisely. Why do not you write your comment as an answer? $\endgroup$
    – Alexei Boulbitch
    Mar 5, 2020 at 10:24
0
$\begingroup$

(Extended comment...)

It might be a good idea to run NIntegrate computational set-up with a few different option settings in order to verify the integral values. For example,

AbsoluteTiming[
 Association@
  Table[i -> 
    NIntegrate[(x + \[Xi])^4/(Sqrt[\[Xi]] (y^2 + (x + \[Xi])^2)^3), {x, -100, 100}, {y, 0, 100}, {\[Xi], 0, 100}, MinRecursion -> i, 
     MaxRecursion -> 100, 
     Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0, 
       "SingularityHandler" -> None, 
       Method -> {"ClenshawCurtisRule", "Points" -> 10}}, 
     AccuracyGoal -> 3, PrecisionGoal -> 5], {i, 1, 10}]
 ]

(* {72.8907, <|1 -> 142.097, 2 -> 158.287, 3 -> 174.651, 
  4 -> 174.651, 5 -> 174.651, 6 -> 174.651, 7 -> 174.651, 
  8 -> 174.651, 9 -> 174.651, 10 -> 174.651|>} *)
Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Thank you I am going to do this. If you have in mind some options that are more promising then others, I would be grateful to know it. $\endgroup$
    – Alexei Boulbitch
    Mar 3, 2020 at 15:03

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.