9
$\begingroup$

The integral

Integrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x}]
(* 512/5355 *)

can be solved analytically.

Trying to apply NIntegrate

NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x}]

fails (although with Method -> "PrincipalValue")

How can I force NIntegrate to calculate? Thanks!

$\endgroup$
12
$\begingroup$

One should indicate the singular line (as it is described in the documentation) to calculate numerically the improper integral under consideration:

 NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x},Exclusions -> {t == x}]

0.0956116

$\endgroup$
  • $\begingroup$ Thanks, how simple... $\endgroup$ – Ulrich Neumann May 23 at 9:49
18
$\begingroup$

How can I force NIntegrate to calculate?

Below are listed "forcing" answers. They should apply in a wide range of situations with minimal understanding of the integrands.

Diagnosing

First let us look at the messages given by NIntegrate:

NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x}]

(* During evaluation of In[18]:= NIntegrate::zeroregion: Integration region {{0.5,1},{1.,0.999999999999999999999999999999975153439150570957241015732418974750}} cannot be further subdivided at the specified working precision. NIntegrate assumes zero integral there and on any further indivisible regions. *)

(* During evaluation of In[18]:= NIntegrate::inumri: The integrand (t^4 x^3)/Sqrt[-t+x] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0.5,1},{0.999999999999999999999999999999975153439150570957241015732418974750,0.999999999999999999990527764909997233148869688962838439242343509680}}. *)

(* NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x}] *)

The messages NIntegrate::zeroregion and NIntegrate::inumri are issued because of the application of the default singularity handler "IMT".

Approaches

1. Using Exclusions is one alternative.

2. Another alternative is to prevent the singularity handler application (and increase MaxRecursion.)

NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x}, 
 Method -> {"GlobalAdaptive", "SingularityHandler" -> None, 
   MaxRecursion -> 120}]

(* 0.0956116 *)

3. A third alternative is to use the tuning parameters for "IMT", if you think "IMT" is beneficial. (Described in NIntegrate's advanced documentation.)

NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x}, 
 Method -> {"GlobalAdaptive", 
   "SingularityHandler" -> {"IMT", "TuningParameters" -> 2}}]

 (* During evaluation of In[16]:= NIntegrate::zeroregion: Integration region {{0.75,1},{1.,0.999999999999999999999999999925933445985376189112492843112221898520}} cannot be further subdivided at the specified working precision. NIntegrate assumes zero integral there and on any further indivisible regions. *)

(* 0.0956116 *)

4. Use more sampling points per integration region.

5. Related to 4, use higher MinRecursion:

NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x}, 
 MinRecursion -> 4]

(* 0.0956116 *)

6. Use higher precision.

7. Switch the order of integration.

$\endgroup$
13
$\begingroup$

I think the problem is that the error estimation at the singularity drives the recursive subdivision too far. In addition to the other methods presented, here are some more.

Use a different rule (with a different error estimator):

NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x}, 
 Method -> "GaussKronrodRule"]
(*  0.0956116  *)

Switch the order of integration:

NIntegrate[(t^4 x^3)/Sqrt[-t + x], {t, 0, 1}, {x, t, 1}]
(*  0.0956116  *)

Use a higher working precision:

NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x}, 
 WorkingPrecision -> 16]
(*  0.09561157754126271  *)

Addendum

I feel NIntegrate should handle the OP's integral without user intervention. The singularity should be easy to identify automatically and easy to handle computationally. I think the problem is that for some unknown reason, the singularity is mishandled and that it could possibly be a bug. Here are three "fixes" for which there is absolutely no mathematical or computational grounding that I can imagine:

NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x}, 
 Exclusions -> x == 100]  (* x == 100 is way outside the interation region *)
(*  0.0956116  *)

NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1}, {t, 0, x}, 
 Exclusions -> t == 100]  (* ditto *)
(*  0.0956116  *)

(* Specify an ordinary point as a singularity in the `x` interval *)
NIntegrate[(t^4 x^3)/Sqrt[-t + x], {x, 0, 1/2, 1}, {t, 0, x}] 
(*  0.0956116  *)

NIntegrate seems to apply "UnitCubeRescaling", which is similar to the following substitution, which I left earlier in a comment:

NIntegrate[
 ((t^4 x^3)/Sqrt[-t + x] /. t -> t x) * Abs@ Det@ D[{x, t x}, {{x, t}}],
 {x, 0, 1}, {t, 0, 1}]
(*  0.0956116  *)

One can partially see into the workings of NIntegrate using IntegrationMonitor:

ireg = NIntegrate[(t^4 x^3)/Sqrt[-t + x],
  {x, 0, 1}, {t, 0, x},
  IntegrationMonitor :> (Return[#, NIntegrate] &)]

Mathematica graphics

If we compare the integrands from my substitution and from the transformation done by NIntegrate, we will see that they are equivalent, although symbolically they are different expressions:

First[ireg]["NumericalFunction"]["FunctionExpression"]
((t^4 x^3)/Sqrt[-t + x] /. t -> t x) Abs@Det@D[{x, t x}, {{x, t}}]
(*
  (t^4 x^8)/Sqrt[x - t x]
  (t^4 x^7 Abs[x])/Sqrt[x - t x]
*)

The only difference is that I wrapped the Jacobian determinant in Abs[]. Since 0 <= x <= 1, there's no significant difference between x and Abs[x]. Or is there?:

NIntegrate[(t^4 x^8)/Sqrt[x - t x], {x, 0, 1}, {t, 0, 1}]

NIntegrate::zeroregion: Integration region...cannot be further subdivided ....

NIntegrate::inumri: The integrand (t^4 x^8)/Sqrt[x-t x] has evaluated to Overflow....

(*  NIntegrate[(t^4 x^8)/Sqrt[x - t x], {x, 0, 1}, {t, 0, 1}]  *)
NIntegrate[(t^4 x^7 Abs[x])/Sqrt[x - t x], {x, 0, 1}, {t, 0, 1}]
(*  0.0956116  *)

In the NIntegrate rescaling, we get the same errors as the OP. In the manual one, with Abs[x], it works without a hitch. It would seem that Abs[x] triggers a different handling of the singularity.

Again, I would repeat that I can see no justification for why the OP's code shouldn't just simply work.

$\endgroup$
  • $\begingroup$ Referenced this answer in mine. (Three times :) $\endgroup$ – Anton Antonov May 23 at 12:18
  • $\begingroup$ @AntonAntonov Thanks. I added something about why I think the behavior of the OP's integral is strange. I know it's been a while since you worked on NIntegrate, but maybe you would have an insight. $\endgroup$ – Michael E2 May 23 at 17:19
  • $\begingroup$ Ok, I will try to investigate/comment in more detail in the next few days. A few of preliminary comments. 1) I had to implement and utilize "UnitCubeRescaling" for variety of reasons. Just the conceptual elegance would have been sufficient, though. Of course NIntegrate had non-symbolic way of handling functional boundaries. 2) The default IMT singularity handler is fairly aggressive in flattening the singularity. If we have arbitrary precision that is kind of fine. But that extra precision hunger has to be curbed. Hence using $MaxExtraPrecision. (cont.) $\endgroup$ – Anton Antonov May 23 at 23:58
  • $\begingroup$ (cont.) 3) Can you repeat your analysis using "SymbolicProcessing" -> 0 ? You might find some answers for the observed manual and automatic transformations differences... 4) To investigate I would combine the explanations from NIntegrate's Advanced Documentation for IMT and "UnitCubeRescaling". I will very likely use the functions of the context "NIntegrateUtilities". $\endgroup$ – Anton Antonov May 24 at 0:04
7
$\begingroup$

You can do a linear variable substitution $y = x - t$, so that the singularity becomes more manageable:

Integrate[((x - y)^4 x^3)/Sqrt[y], {x, 0, 1}, {y, 0, x}]
(* 512/5355 *)

NIntegrate[((x - y)^4 x^3)/Sqrt[y], {x, 0, 1}, {y, 0, x}]
(* 0.0956116 *)

Or even eliminate the singularity completely by substituting $z = \sqrt{x-t}$:

Integrate[2 (x - z^2)^4 x^3, {x, 0, 1}, {z, 0, Sqrt[x]}]
(* 512/5355 *)

NIntegrate[2 (x - z^2)^4 x^3, {x, 0, 1}, {z, 0, Sqrt[x]}]
(* 0.0956116 *)

In my experience this way of proceeding is often much more fruitful than addressing the technical/methodical difficulties of NIntegrate.

$\endgroup$
  • $\begingroup$ I was just about to add your new, second method to my answer. +1 $\endgroup$ – Michael E2 May 23 at 12:37
  • 1
    $\begingroup$ I wish someone would explain the downvote. The IMT singularity handler works better when the singularity aligns with one of the coordinate axes. The change of variables in the first example here does that. I was going to add to my answer the example NIntegrate[((t^4 x^3)/Sqrt[-t + x] /. t -> u x) Abs@ Det@ D[{x, u x}, {{x, u}}], {x, 0, 1}, {u, 0, 1}], which is another change of variables, {x, t} -> {x, u x}, that accomplishes a similar alignment and is essentially the same idea. $\endgroup$ – Michael E2 May 23 at 12:43
  • $\begingroup$ @MichaelE2 I downvoted the original version because it was too short of explanations, just proposing a substitute. The new version discusses the singularity elimination. (IMT "just" does singularity flattening.) Nevertheless, one of the ideas behind "big functions" like NIntegrate (and NDSolve, NMinimize, etc.) is that we should not think that much when using them, their frameworks should support "simple user" usage. Meaning 1) automatic method (options) selection, or 2) effective tweaking with minimal mathematics knowledge and/or understanding. $\endgroup$ – Anton Antonov May 23 at 13:09
  • 1
    $\begingroup$ @AntonAntonov I actually agree with your downvote, and believe the automatic methods should be promoted instead. For those readers who have enough analysis skill to get a substitution going, what I wrote is probably rather trivial. $\endgroup$ – Roman May 23 at 15:16
  • 1
    $\begingroup$ @AntonAntonov Thanks for the explanation. I think it helps the site if shortcomings and potential improvements are pointed out. $\endgroup$ – Michael E2 May 23 at 16:18

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.