4
$\begingroup$

Executing the following code, NDSolve begins to solve a system of two coupled ordinary differential equations for the dependent variables $\rho_0(t)$ and $\lambda_k(t)$. Around $t=-2$ NIntegrate warns that

NIntegrate::inumri: The integrand [...] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0.,4.64782*10^14}}.

causing NDSolve to abort. Could some NIntegrate-whiz tell me what's wrong with my integrand? I looked at it in various 3d plots and it appears completely regular. I tried rerunning the calculation with different values for WorkingPrecision and MaxRecursion in case NIntegrate was having difficulties with small numbers, but without success. Any other suggestions how to diagnose this further and get NDSolve to solve my ODEs would be much appreciated!

sqrt[x_?NumericQ,y_?NumericQ]=Piecewise[{{I Sqrt[-x],Re[x]<0&&Re[y]>=0},{-I Sqrt[-x],Re[x]<0&&Re[y]<0}},Sqrt[x]];
Derivative[1,0][sqrt][x_,y_]=1/(2sqrt[x,y]);
Derivative[0,1][sqrt][x_,y_]=0;

αk[1,1][m2_,g2_,z_,c_]:=1/2 (1/c+m2/z-I g2/z)+I sqrt[1/(c z)-1/4 (1/c-m2/z+I g2/z)^2,-(1/c-m2/z)];
αk[1,2][m2_,g2_,z_,c_]:=1/2 (1/c+m2/z+I g2/z)-I sqrt[1/(c z)-1/4 (1/c-m2/z-I g2/z)^2,(1/c-m2/z)];
αk[2,1][m2_,g2_,z_,c_]:=1/2 (1/c+m2/z-I g2/z)-I sqrt[1/(c z)-1/4 (1/c-m2/z+I g2/z)^2,-(1/c-m2/z)];
αk[2,2][m2_,g2_,z_,c_]:=1/2 (1/c+m2/z+I g2/z)+I sqrt[1/(c z)-1/4 (1/c-m2/z-I g2/z)^2,(1/c-m2/z)];

w={1,1,1,1,-1,-1,-1,-1,-2,2};
μ=Join[Flatten@Table[αk[i,j][m2/s^2,g2/s^2,z,c]s^2,{s,{k,Λ}},{i,2},{j,2}],{k^2/c,Λ^2/c}];

ΔkerT=T Sum[w[[j]]Log[Exp[1/T Sqrt[p^2+μ[[j]]]]-1],{j,10}];
kerT[0][m2_,g2_,z_,c_,T_,p_]=k D[ΔkerT,k]/.k->1;
kerT[j_][m2_,g2_,z_,c_,T_,p_]:=(-1)^j Derivative[j,0,0,0,0,0][kerT[0]][m2,g2,z,c,T,p];

tfI[j_][m2_,g2_,z_,T_,c_,d_,pmax_,ops___]:=(2π^(d/2))/((2π)^d Gamma[d/2]) NIntegrate[kerT[j][m2,g2,z,c,T,p] p^(d-1),{p,0,pmax},ops];

ρ0Flow[ρ0_?NumericQ,λk_?NumericQ,ηk_?NumericQ,g2_?NumericQ,Z1_?NumericQ,T_?NumericQ,c_?NumericQ,d_?NumericQ,N_?NumericQ,pmax_,ops___]:=-(2+ηk)ρ0+(3/2 tfI[1][2ρ0 λk,g2,Z1,T,c,d,pmax,ops]+(N-1)/2tfI[1][0,0,1,T,c,d,pmax,ops])
λkFlow[ρ0_?NumericQ,λk_?NumericQ,ηk_?NumericQ,g2_?NumericQ,Z1_?NumericQ,T_?NumericQ,c_?NumericQ,d_?NumericQ,N_?NumericQ,pmax_,ops___]:=2ηk λk+λk^2 (9/2 tfI[2][2ρ0 λk,g2,Z1,T,c,d,pmax,ops]+(N-1)/2tfI[2][0,0,1,T,c,d,pmax,ops])

AbsoluteTiming@Module[{runner = 0, counter = 0},
  With[{ηk = 0, g2 = 0, Z1 = 1, T = 0.1 E^t, c = 1, d = 3, N = 2,
     pmax = ∞}, 
   run[1] = 
    NDSolve[{ρ0'[
        t] == ρ0Flow[ρ0[t], λk[t], ηk, g2, Z1, 
        T, c, d, N, pmax], λk'[
        t] == λkFlow[ρ0[t], λk[t], ηk, g2, 
        Z1, T, c, d, N, pmax], ρ0[0] == 0.02, λk[0] == 
       0.5}, {ρ0, λk}, {t, -10, 0}, 
     StepMonitor :> 
      counter++ If[Abs[t] > runner, 
        Print@Chop[{counter, Round[t, 1], 
           E^(2 t) ρ0[t], λk[t]}]; runner++]]]]
$\endgroup$
  • $\begingroup$ @MichaelE2 That's odd. I reach $t = -2$ in about 40 seconds using MMA 11.2 on macOS 10.13. What's your system? $\endgroup$ – Casimir Nov 18 '17 at 20:09
  • 1
    $\begingroup$ I thought I was using 11.2, but it turns out I pasted it into the beta...Well, that's something to look forward to, hmm? $\endgroup$ – Michael E2 Nov 18 '17 at 20:28
3
+25
$\begingroup$

NIntegrate fails when the arguments of exponentials in the integrand become too large, causing overflows. This can be seen by removing unnecessary code from AbsoluteTiming@NDSolve ... and enclosing it in EvaluationData as follows.

With[{ηk = 0, g2 = 0, Z1 = 1, T = 0.1 E^t, c = 1, d = 3, n = 2, pmax = ∞},  
    EvaluationData[NDSolve[{ρ0'[t] == ρ0Flow[ρ0[t], λk[t], ηk, g2, Z1, T, c, d, n, pmax], 
    λk'[t] == λkFlow[ρ0[t], λk[t], ηk, g2, Z1, T, c, d, n, pmax], 
    ρ0[0] == 0.02, λk[0] == 0.5}, {ρ0, λk}, {t, -10, 0}]]];

Which produces an enormous expression and many error messages. (Note that N has been replaced by n everywhere in the code in the question, because N is a Mathematica function. Using this reserved symbol is not, however, the source of the problem in the question.) The NIntegrate integrand that caused the error then is obtained from

Cases[%, Hold[Message[NIntegrate::inumri, z__]] -> z, Infinity, 1]

which also is very large but contains 18 terms of the general form

-(((29.0983 + 50.3998 I) E^(116.393 Sqrt[1/2 - (I Sqrt[3])/2 + p^2]))/
    ((-1 + E^(116.393 Sqrt[1/2 - (I Sqrt[3])/2 + p^2])) (1/2 - (I Sqrt[3])/2 + p^2)^(5/2)))

Moreover, the limits of integration over p are {0., 4.64782*10^14}. Thus, the exponentials overflow. Probably, the solution is to restructure the integrand so that it consists instead of terms of the general form

-(((29.0983 + 50.3998 I)/
   ((-E^(-116.393 Sqrt[1/2 - (I Sqrt[3])/2 + p^2]) + 1) (1/2 - (I Sqrt[3])/2 + p^2)^(5/2)))

which should avoid the overflows. I hope that this forensic analysis it helpful and regret that I do not have time to pursue the problem futher.

$\endgroup$
  • $\begingroup$ Thanks! This already clears up quite a few things. Cases[%, Hold[Message[...]]] is a great way to pick out problematic terms out of an unwieldy error-throwing expression (+1 for that alone). Do you know if there is any way to tell Mathematica to do the restructuring that avoids overflow automatically? $\endgroup$ – Casimir Nov 20 '17 at 10:07
  • 1
    $\begingroup$ @Casimir You might process each of the 10 terms in ΔkerT separately, and only combine them in the formula for tfI. This will allow you to identify the offending exponentials more readily and recast them to avoid overflows, perhaps by using Together to put the each term over a common denominator, multiply itsNumerator and Denomiator by Exp[-...], Simplify`, and only then recombine the numerator and denominator. Some experimentation may be required to determine where in the processing to best do this. $\endgroup$ – bbgodfrey Nov 20 '17 at 18:33

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.