0
$\begingroup$

While integrating the following function,

me = 511000/29979200(*eV*);
re = 2.81794*10^-13(*cm*);
h = 4.135667516 * 10^-15;(*eV*)
c = 29979200;(*cm/s*)

Zg = 79;
Ag = 196.96655;
\[Rho] = 19.320(*g/\[Mu]m^3*);

\[Xi][\[Nu]_?NumericQ, T_?NumericQ] := (100 me c^2 h \[Nu])/(
 T (T - h \[Nu] ) Zg^(1/3));(*Per Bremmstrahlung*)
\[Phi]1[T_?NumericQ, \[Nu]_?NumericQ] := 
 20.863 - 2 Log[1 + (0.55846 \[Xi][T, \[Nu]])^2] - 
  4 (1 - 0.6 Exp[-0.9 \[Xi][T, \[Nu]]] - 
     0.4 Exp[-105 \[Xi][T, \[Nu]]]);(*Per Bremmstrahlung*)
\[Phi]2[T_?NumericQ, \[Nu]_?NumericQ] := \[Phi]1[T, \[Nu]] - 
  2/3 1/(1 + 6.5 \[Xi][T, \[Nu]] + 
    6 \[Xi][T, \[Nu]]^2);(*Per Bremmstrahlung*)
f[Z_?NumericQ ] := (Z/137)^2 (1/(1 - (Z/137)^2) + 0.20206 - 
    0.0369 (Z/137)^2 + 0.0083 (Z/137)^4 - 
    0.002 (Z/137)^6);(*Per Bremmstrahlung*)

dEdx2[T_?NumericQ] := 
  6.022*10^23*\[Rho]/
   Ag NIntegrate[(4 Zg^2 re/
        137 1/\[Nu] (1 + (T/(T - h \[Nu]))^2) (\[Phi]1[T, \[Nu]]/4 - 
          1/3 Log[Zg] - f[Zg]) - 
       2/3 T/(T - 
         h \[Nu]) (\[Phi]2[T, \[Nu]] - 1/3 Log@Zg - f[Zg])), {\[Nu], 
      0, T/h}];

ParallelTable[{j, dEdx2[j]}, {j, 10000, 2*10^6, 10000}]

I get a very long series of errors, the first of which is

NIntegrate::zeroregion Integration region {{2.41799*10^18,2.41798934786516787199999999993411568411875634969888072969436258221*10^18}} cannot be further subdivided at the specified working precision. NIntegrate assumes zero integral there and on any further indivisible regions.

What does it mean, and how do I get rid of it?

$\endgroup$
  • 2
    $\begingroup$ You have constants that differ by dozens of orders of magnitude; numerical results will be catastrophically bad. Please, work in natural units, where all parameters are of order one. $\endgroup$ – AccidentalFourierTransform Jul 7 '18 at 15:03
  • $\begingroup$ @AccidentalFourierTransform I understand your point, but if you look at integral, all of the constant except for h are multiplicative. $\endgroup$ – mattiav27 Jul 7 '18 at 16:15
  • 1
    $\begingroup$ So what? If most constants can be factored out, do so manually. The error message is telling you explicitly that the integral is numerically unstable. The very first step is to make it stable by working in natural units. Being a phd student in, I presume, physics/chemistry/engineering, you should know better by now: always use natural units, and especially when doing numerics! $\endgroup$ – AccidentalFourierTransform Jul 7 '18 at 16:18
3
$\begingroup$

Obtaining a zero length integration region means that most likely the singularity handler is not aggressive enough.

So, one way to deal with your problem is to try to change the singularity handler and/or the precision goal. Here is an answer doing that for a similar question.

Further, you can rationalize all the numerical entities and use higher precision. (Not done here.)

First let us redefine dEdx2:

ClearAll[dEdx2]
Options[dEdx2] = Options[NIntegrate];
dEdx2[T_?NumericQ, opts : OptionsPattern[]] := 
  6.022*10^23*\[Rho]/
    Ag NIntegrate[(4 Zg^2 re/
        137 1/\[Nu] (1 + (T/(T - h \[Nu]))^2) (\[Phi]1[T, \[Nu]]/4 - 
         1/3 Log[Zg] - f[Zg]) - 
      2/3 T/(T - h \[Nu]) (\[Phi]2[T, \[Nu]] - 1/3 Log@Zg - 
         f[Zg])), {\[Nu], 0, T/h}, opts];

This takes less than 40s on my laptop (without any messages):

AbsoluteTiming[
 res3 = 
   Table[
     {j, 
      dEdx2[j, 
       Method -> {"GlobalAdaptive", "SingularityHandler" -> "DoubleExponential"}, 
       MaxRecursion -> 20, PrecisionGoal -> 3]}, 
     {j, 10000, 2*10^6, 10000}];
 ]
Short[res3]

(* {37.456, Null} *)

(* {{10000, -6.33661*10^43}, {20000, -1.26732*10^44}, {30000, \
-1.92223*10^44}, {40000, -2.53464*10^44}, {50000, -3.18431*10^44}, \
{60000, -3.84447*10^44} ... } *)

The computation above with PrecisionGoal->4 takes ~15 times longer; the results are assigned to res4. Here are the relative errors:

Through[{Min, Max, Mean, Median}[
  Abs[(res4[[All, 2]] - res3[[All, 2]])/res4[[All, 2]]]]]

(* {0., 0.0134694, 0.00151681, 0.000671349} *)
| improve this answer | |
$\endgroup$
  • $\begingroup$ What is res4 ? $\endgroup$ – Mariusz Iwaniuk Jul 8 '18 at 7:08
  • $\begingroup$ What result you assigned? $\endgroup$ – Mariusz Iwaniuk Jul 8 '18 at 11:08

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.