4
$\begingroup$

Consider the following example:

(* This in the first cell *)
ACf[z_] = Sqrt[2 + M - (M z^2)/2];
ABf[z_] = Sqrt[4 + 2 M - M z^2]/Sqrt[2 + M];
J[h_, z_] = h/Sqrt[2 M + 4 - 2 h^2 - M z^2];


(* This in the second cell *)
Block[{M = 10},
 NIntegrate[
  J[(Min[z, ACf[z]] - ABf[z]) t + ABf[z], z], {z, Sqrt[2 + M]/Sqrt[
   1 + M], (Sqrt[2] Sqrt[2 + M])/Sqrt[M]}, {t, 0, 1}, 
  WorkingPrecision -> 20, Method -> "TrapezoidalRule"]
 ]

With fresh kernel I evaluate the contents of the first cell. Then after the evaluation is finished, I evaluate the second one. And the result is:

NIntegrate::nlim: z = Sqrt[2.0000000000000000000+M]/Sqrt[1.0000000000000000000+M] is not a valid limit of integration.

NIntegrate[ J[(Min[z, ACf[z]] - ABf[z]) t + ABf[z], z], {z, Sqrt[2 + M]/Sqrt[ 1 + M], (Sqrt[2] Sqrt[2 + M])/Sqrt[M]}, {t, 0, 1}, WorkingPrecision -> 20, Method -> "TrapezoidalRule"]

In some cases, instead the kernel just quits without having printed anything.

I reproduce this in Mathematica 11.0, but not in 9.0.

What's happening here? In Mathematica 9.0 I get just NIntegrate::slwcon, not such strange errors as in 11.0.

$\endgroup$
1
  • $\begingroup$ I think there is a bug with NIntegrate here $\endgroup$
    – chuy
    Commented Nov 11, 2016 at 18:32

2 Answers 2

3
$\begingroup$

It has to do with the fact that when specifying Method -> "TrapezoidalRule", this integral simply won't evaluate at all. What it doesn't, I do not know. But it explains the behaviour you see.

Let's drop the Block and evaluate

M = 10;

NIntegrate[J[(Min[z, ACf[z]] - ABf[z]) t + ABf[z], z],
 {z, Sqrt[2 + M]/Sqrt[1 + M], (Sqrt[2] Sqrt[2 + M])/Sqrt[M]}, {t, 0, 
  1},
    WorkingPrecision -> 20, Method -> "TrapezoidalRule"]
(* NIntegrate[
 J[(Min[z, ACf[z]] - ABf[z]) t + ABf[z], z], {z, Sqrt[2 + M]/Sqrt[
  1 + M], (Sqrt[2] Sqrt[2 + M])/Sqrt[M]}, {t, 0, 1}, 
 WorkingPrecision -> 20, Method -> "TrapezoidalRule"] *)

It's returned as is.

What happens if we M =., then

Block[{M=10},
   NIntegrate[...]
]

?

Well, it will first evaluate with M having a value within a block. The result is the same NIntegrate, with no change (i.e. no evaluation). This is what is returned from the Block. But when it's returned, M loses its value, so the expression evaluates again, a second time. This time M has no value and you see the error you mentioned.

Removing this Method specification allows the evaluation to finish.

There is no mystery related to Block. The real question is: why doesn't this integral evaluate at all? I don't remember ever seeing NIntegrate return as entered ... Integrate does that often when it can't compute the result. But I have not seen it with NIntegrate.

$\endgroup$
4
  • $\begingroup$ But when NIntegrate returns unevaluated result, why is the evaluation then restarted? $\endgroup$
    – Ruslan
    Commented Nov 11, 2016 at 16:11
  • $\begingroup$ @Ruslan I don't have a full understanding of how this works, but it seems clear enough that if it didn't work this way, that would mess things up. Consider f[x_?NumericQ] := x and x=1; Then try Block[{x}, Echo@f[x]]. It would be quite weird if the evaluation were not restarted after f[x] is returned from Block. What exactly triggers the restart? I don't know. Maybe just a simple value change (x changed). $\endgroup$
    – Szabolcs
    Commented Nov 11, 2016 at 16:32
  • $\begingroup$ @Ruslan See also Update. It deal exactly with situations when the evaluation should be restarted but the system doesn't do this on its own. I am not experienced with Update. $\endgroup$
    – Szabolcs
    Commented Nov 11, 2016 at 16:34
  • 1
    $\begingroup$ Another factor is NIntegrate is HoldAll, so that its arguments are returned with M unevaluated. -- I agree with your comment that leaving Block resets M and triggers a new evaluation cycle. This is explicitly mentioned at the end of the tutorial Blocks and Local Values. $\endgroup$
    – Michael E2
    Commented Nov 12, 2016 at 13:34
2
$\begingroup$

Too long for a comment: What I present below is just a side issue to what's going on with Block, which Szabolcs has already explained. The integrand is not a good candidate for the trapezoidal rule, since it has an infinite singularity at t == 1 for z > Sqrt[2]. The integral is convergent, but the trapezoidal rule has trouble approximating its value.

The following is a simpler example of the numerical issue with NIntegrate and the "TrapezoidalRule". Having a function like Min, which is decomposed into piecewise components by NIntegrate, seems to be a key factor in reproducing the problem.

NIntegrate[Min[y, 1/2]/Sqrt[1 - x],
 {y, 0, 1}, {x, 0, 1},
 WorkingPrecision -> 20, Method -> "TrapezoidalRule"]

It also crashes the kernel sometimes. It seems always to crash if you try to Trace[] it.

Note that like the OP's integrand, this has an infinite singularity, and while the integral is convergent, it returns unevaluated. Also note the single integral

NIntegrate[1/Sqrt[1 - x], {x, 0, 1}, WorkingPrecision -> 20, Method -> "TrapezoidalRule"]

evaluates to a number, albeit with convergence warnings.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.