WorkingPrecision causes issue in the NIntegrate

Question

I really can't figure out why my code sometimes is not working. My integrals involve two variables (k and kz). The integration range for both of them are from zero to infinity. I found out that in some cases (not often) when the working precision is in some specific range, Mathematica gives me an error message saying that "the integrand has evaluated to non-numerical values...". However, if I change WorkingPrecision to higher or lower value, it works just fine. I am not sure why. Hope you can help me with this.

I copy my code here. Thanks a lot!

Clear[sofindmassz];
sofindmassz[λ_, T_, μ_, Δ_, opt : OptionsPattern[{workprec -> 16}]] := 
 Module[{wp, ϵ, ξ, E, dϵ, ddϵ, fBCSE, fBCSξ, f1d, mInth, mInt}, 
  wp = OptionValue[workprec]; ϵ[h_] := k^2 + kz^2 + 2 λ h k; ξ[h_] := ϵ[h] - μ; 
  E[h_] := √(ξ[h]^2 + Δ^2); 
  dϵ[h_] := 2 kz; ddϵ = 2; 
  f1d[h_] := -Exp[ξ[h]/T]/(T (Exp[ξ[h]/T] + 1)^2); 
  fBCSE[h_] := Tanh[E[h]/(2 T)]; 
  fBCSξ[h_] := Tanh[ξ[h]/(2 T)]; 
  mInth[h_] := k ((2 f1d[h] + E[h]/Δ^2 ((1 + ξ[h]^2/E[h]^2) (fBCSE[h] - fBCSξ[h]) +
      (1 - ξ[h]/E[h])^2 fBCSξ[h])) (dϵ[h])^2 - 1/2 (fBCSξ[h] - ξ[h]/E[h] fBCSE[h]) ddϵ); 
  mInt = mInth[1] + mInth[-1]; 
  Return[-Quiet[
      NIntegrate[
       NIntegrate[mInt, {k, 0, kz}, 
        Method -> {Automatic, "SymbolicProcessing" -> 0}, 
        WorkingPrecision -> wp], {kz, 0, ∞}, 
       Method -> {Automatic, "SymbolicProcessing" -> 0}, 
       WorkingPrecision -> wp]] - 
    Quiet[NIntegrate[
      NIntegrate[mInt, {kz, 0, k}, 
       Method -> {Automatic, "SymbolicProcessing" -> 0}, 
       WorkingPrecision -> wp], {k, 0, ∞}, 
      Method -> {Automatic, "SymbolicProcessing" -> 0}, 
      WorkingPrecision -> wp]]]]

For this code, if I run the following values,

sofindmassz[0, 0.14921714620005236`, 0.07455393513003296`, 1.016522853606922`, workprec -> 15]

the error messages pops up.

If I change workprec to 14. It works! Don't know what's wrong.

Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) 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! — Michael E2, Dec 25 '14 at 02:04
The problem relates to the nesting of NIntegrate. See 69638. Why it singles out workprec -> 15 is unclear to me, although it may be related to shifting from double precision to higher. — bbgodfrey, Dec 25 '14 at 03:22
Other values of WorkingPrecision produce warning messages too, but Quiet suppresses them. — bbgodfrey, Dec 25 '14 at 12:48
@Chien-TeWu, if bbgodfrey's answer was helpful, you should consider upvoting it. And accept an answer, if appropriate, too. (See the welcome comment above.) — Michael E2, Jan 24 '15 at 18:33
@MichaelE2 I'd like to upvote it. However, I don't have enough reputation. Could you let me know if there is any other way to do that? Thanks! — Chien-Te Wu, Jan 25 '15 at 21:08
I see. Too bad for now. I guess you need four more points. (I already upvoted your question, so I can't do anything.) I thought maybe you could upvote answers to your question. You certainly should be able to accept an answer regardless of reputation, but that will give you only two points. — Michael E2, Jan 25 '15 at 21:18

score 2 · Answer 1 · answered Jan 24 '15 at 18:31

One way to address the example is to provide input of the same or higher (arbitrary) precision,

sofindmassz[0, 0.14921714620005236`15, 0.07455393513003296`15, 1.016522853606922`15, 
  workprec -> 15]
(*  0.643229985259241  *)

One can get overflow, probably from

f1d[h_] := -Exp[ξ[h]/T]/(T (Exp[ξ[h]/T] + 1)^2)

because the size of ξ[h] can be quite large ~ k^2 when integrating to infinity. The alternative

f1d[h_] := -Exp[-ξ[h]/T]/(T (Exp[-ξ[h]/T] + 1)^2)

leads to underflow and other problems. Depending on how I played with it, it might lead to a kernel crash from running out of memory or take so long I aborted it. It you plot the integrand, it looks good and smooth up to k == 1000; but as k exceeds 10000 it starts to get bumpy and the bumps grow larger.

Here are some sample values of the integrand:

Grid@Table[mInt, {k, 1000000, 10000000, 1000000}, {kz, 0, k, 1000000}]

Mathematica graphics

I get the same results no matter what precision I use for the parameters in sofindmassz. While they are not the same as the sample points used by NIntegrate, the occasional extremely large values (compared to the estimated value of the integral) are disturbing.

It seems likely that numerical precision, underflow and/or overflow play a role in the evaluation of this integral. I've been unable to chase down the clear reasons for the difficulty. As a result I cannot comment on the accuracy of the result above, other than to say it was the most common result (up to the first 8 or digits).

Using a double integral

NIntegrate[mInt, {k, 0, ∞}, {kz, 0, k}, 
 Method -> {Automatic, "SymbolicProcessing" -> 0}, 
 WorkingPrecision -> wp]

instead of nested integrals did not seem to help. (It behaves somewhat differently, but it does not behave better.)

By the way, I am not a fan of using Quiet, especially to suppress numerics-related messages, unless I know why the warnings and errors are occurring and in fact expect them. Usually one needs to understand why Mathematica is warning of a possible numerical error. (If you get a lot, then you can set the preferences to log them in the messages window.)

score 1 · Answer 2 · answered Dec 25 '14 at 13:17

Replace your Return line (the final line in your code by

fkz[z_?NumericQ] := NIntegrate[mInt, {k, 0, z}, WorkingPrecision -> wp];
fk[z_?NumericQ] := NIntegrate[mInt, {kz, 0, z}, WorkingPrecision -> wp];
Quiet[-NIntegrate[fkz[kz], {kz, 0, \[Infinity]}, WorkingPrecision -> wp] - 
  NIntegrate[fk[k], {k, 0, \[Infinity]}, WorkingPrecision -> wp]]]

to eliminate the unwanted warning message. Creating separate functions fkz and fk with argument z_?NumericQ is the solution. The other changes merely delete code that no longer is necessary. Note, however, that Quiet hides still other warning messages.

WorkingPrecision causes issue in the NIntegrate

2 Answers2