# why the integral without singularity does not converge

Here is an integral

NIntegrate[1/(px^2 + py^2 + 1)/((px - qx)^2 + (px - qy)^2 +
1)/(qx^2 + qy^2 + 1), {px, 0, 10}, {py, 0, 10}, {qx, 0, 10}, {qy, 0, 10}]


the integrant are all sufficiently small, but still MMA give a result and say the integral coverage too slow. I tried different method, none of them works. I could not understand what causes it? How to modify the code?

• Often one can ignore NIntegrate::slwcon if there are no other errors. (So you left out an important component of the question. You also left out the error message names, which are helpful to include because then people with similar problems can search the site for them.) – Michael E2 Jul 21 '17 at 17:45

This works well.

NIntegrate[1/(px^2 + py^2 + 1)/((px - qx)^2 + (px - qy)^2 + 1)/(qx^2 + qy^2 + 1),

{px, 0, 10}, {py, 0, 10}, {qx, 0, 10}, {qy, 0, 10},Method -> "CartesianRule"] // Timing


{58.875, 2.84714}

• Thanks for the answer, what is the reason this method work? Do you have some advice to speed up the integral? – p.s Jul 21 '17 at 5:22
• @p.s. See reference.wolfram.com/language/tutorial/… for more details. No, I don't have any advice to speed up the integral. In my opinion, one minute doesn't last too long. – user64494 Jul 21 '17 at 5:30

Lowering the precision goal slightly from 6 to 5 (perhaps advisable in a high-dimensional integral) and increasing MinRecursion avoids both the NIntegrate::slwcon and NIntegrate::eincr warnings:

NIntegrate[
1/(px^2 + py^2 + 1)/((px - qx)^2 + (px - qy)^2 + 1)/(qx^2 + qy^2 + 1),
{px, 0, 10}, {py, 0, 10}, {qx, 0, 10}, {qy, 0, 10},
PrecisionGoal -> 5,
MinRecursion -> 1] // AbsoluteTiming
(*  {5.563214, 2.8471463223061835}  *)


Update: Alternative

FWIW, you can do an expensive (i.e., takes time) symbolic integration on one or two of the variables, followed by a numerical integration. E.g, integrating over {py, 0, 10} symbolically first takes about as long as the above and needs no other special numerical processing. The numerical advantage here is reducing the dimension of the integral. (In fact, if you integrate over qx and then py, then the two-dimensional integral is quite fast and unproblematic; however, the symbolic integration takes about 17 seconds.)

NIntegrate[
Evaluate@Integrate[
1/(px^2 + py^2 + 1)/((px - qx)^2 + (px - qy)^2 + 1)/(qx^2 + qy^2 + 1),
{py, 0, 10},
Assumptions -> 0 < px < 10 && 0 < py < 10 && 0 < qx < 10 && 0 < qy < 10]
, {px, 0, 10}, {qx, 0, 10}, {qy, 0, 10}] // AbsoluteTiming
(*  {4.173041, 2.8471449687848587}  *)


We also get the extra digit of accuracy @user64494 got.

Update 2: Third try, inspired by Mariusz Iwaniuk

There is a slight "spike" in the value of the integrand near the subspace px == qx == qy. This might be why NIntegrate has trouble getting the error estimate to satisfy precision goals of 6 or more. Some weak evidence is that feeding the subspace as an "Exclusion" with the "LocalAdaptive" method to NIntegrate gives a fairly quick solution agreeing with the other solutions to slightly more than six digits.

NIntegrate[
1/(px^2 + py^2 + 1)/((px - qx)^2 + (px - qy)^2 + 1)/(qx^2 + qy^2 + 1),
{px, 0, 10}, {py, 0, 10}, {qx, 0, 10}, {qy, 0, 10},
Method -> "LocalAdaptive", Exclusions -> {px == qx == qy}] // AbsoluteTiming
(*  {3.497589, 2.8471576198024686}  *)

• A professional answer. Can you explain why AccuracyGoal -> 5 doesn't work well here? Thanks in advance. – user64494 Jul 21 '17 at 6:00
• @user64494 Thanks. AccuracyGoal sets a goal for the number of digits after the decimal point. When the result is 1, then Precision and Accuracy are the same. The corresponding AccuracyGoal for this integral would be about Accuracy[2.84714632230618355] or 4.54559. (The precision and accuracy goals combine, so that the absolute error estimate should be less than the sum of the errors implied by the two goals, and the correspondence is not quite as simple as I describe.) – Michael E2 Jul 21 '17 at 14:51
• @Michael E2 Thanks. You guys are really helpful. – p.s Jul 21 '17 at 20:48

NIntegrate[
1/(px^2 + py^2 + 1)/((px - qx)^2 + (px - qy)^2 + 1)/(qx^2 + qy^2 + 1),
{px, 0, 10}, {py, 0, 10}, {qx, 0, 10}, {qy, 0, 10},
Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}] // Timing

(*{1.45313, 2.84718}*)


or:

NIntegrate[1/(px^2 + py^2 + 1)/((px - qx)^2 + (px - qy)^2 + 1)/(qx^2 + qy^2 +
1), {px, 0, 10}, {py, 0, 10}, {qx, 0, 10}, {qy, 0, 10},
PrecisionGoal -> 5, MinRecursion -> 1, Method -> {"LocalAdaptive",
"SymbolicProcessing" -> 0}] // AbsoluteTiming

(*{0.687526, 2.84731}*)

• (+1) Nice. Seems to underestimate the error, though. Have to bump up the PrecisionGoal` to get a comparable quality result, and then it's not really faster (or so it seems to me). – Michael E2 Jul 21 '17 at 17:42
• @Mariusz Iwaniuk Thanks. – p.s Jul 21 '17 at 20:48