NIntegral failing to converge

Question

I have the following integrals that I am trying to calculate:

a[q_, E_, K_, y_] = 1 + q/(2*K)*Log[(q^2 - K*q - E)/(q^2 + K*q - E)];

f[Epol_?NumericQ, K_?NumericQ, y_?NumericQ] := NIntegrate[a[q, Epol, K, y],{q, y, Infinity}, Method -> "GlobalAdaptive", WorkingPrecision -> 20, MaxPoints -> 100];

rhs[x_?NumericQ, y_?NumericQ, Epol_?NumericQ] := 4/Pi *NIntegrate[K^2/(x - 2/Pi *y + K^2/4 + Epol - 2/Pi*f[Epol, K, y]), {K, 0,y},Method -> "GlobalAdaptive", MaxPoints -> 100];

Then I want to use that to solve an equation as the following:

FindRoot[E == rhs[0.2, 1, E], {E, -0.01}]

I'm getting always the convergence error in the integral:

NIntegrate::maxp: The integral failed to converge after 121 integrand evaluations. NIntegrate obtained 0.0818604 +0.0435709 I and 0.0005506179310998352` for the integral and error estimates.

I tried to change the maximum number of evaluation points allowed, the working precision, and the accuracy/precision goals, nothing has worked out well for values of x more than 0.1, any other solutions to try are welcome.

Even the f function really struggles to return a value when I plug in some random input. Are you confident that these integrals are convergent? Can you provide a set of values for the parameters for which f returns without errors or warnings? — MarcoB, Dec 02 '19 at 23:24
Avoid using capital letters for variables, especially ones used by Mathematica for other purposes, such as E (base of the natural logarithm) and K (dummy summation/integration variable). — Michael E2, Dec 03 '19 at 14:11

score 1 · Answer 1 · answered Dec 03 '19 at 12:53

1

Multiple problems with your code. Here is better code:

Clear[a, f, rhs]

a[q_, E_, K_, y_] := 1 + q/(2*K)*Log[(q^2 - K*q - E)/(q^2 + K*q - E)];

f[Epol_?NumericQ, K_?NumericQ, y_?NumericQ] := 
  NIntegrate[a[q, Epol, K, y], {q, y, Infinity}, 
   Method -> "GlobalAdaptive", WorkingPrecision -> 20, 
   PrecisionGoal -> 2];

rhs[x_?NumericQ, y_?NumericQ, Epol_?NumericQ] := 
  4/Pi*NIntegrate[
    K^2/(x - 2/Pi*y + K^2/4 + Epol - 2/Pi*f[Epol, K, y]), {K, 0, y}, 
    Method -> "GlobalAdaptive", PrecisionGoal -> 2];

FindRoot[var == rhs[0.2, 1, var], {var, -0.01}] 

(* Lots of messages.. *)
(* {var -> -0.010069} *)

[...] nothing has worked out well for values of x more than 0.1 [...]

Is this something you expect:

FindRoot[var == rhs[12, 1, var], {var, -0.01}]

(* {var -> 0.036143} *)

answered Dec 03 '19 at 12:53

Anton Antonov

37,787
3
100
178

I don't clearly see what changes have you made except for reducing the precision goal, which I have already tried. I think I have figured out the problem, it has to do with how Mathematica searches for solutions in FindRoot. A possible workaround is to put a starting point very close to the result to be able to avoid these convergence problems. For instance FindRoot[var == rhs[0.2, 1, var], {var, -0.44}]gives a value which is logical (based on what I know of the behavior of this equation) without any error messages. – Dec 03 '19 at 13:03
But this is not a trustworthy solution since it means the equation has a solution around this value which leads to a convergent integral. Also, I might not have always an intuition of the logical outcome, which means predicting the starting point becomes harder. Anyway, I will try to write a code in python since I'm better at debugging there, I don't think I will get it to work with mathematica. – Dec 03 '19 at 13:09
@RagheedAlhyder "I don't clearly see what changes have you made except for reducing the precision goal [...] " -- look more carefully. – Anton Antonov Dec 03 '19 at 13:19
Sorry I don't see it, is it the removal of maximum points? Anyway the performance has stayed the same. – Dec 03 '19 at 13:33

score 1 · Answer 2 · answered Dec 07 '19 at 11:44

Do as much of integration you can analytically.

a[q_, EE_, K_, y_] = 
  1 + q/(2*K)*Log[(q^2 - K*q - EE)/(q^2 + K*q - EE)] // 
FullSimplify[#, 0 < K < y < 1 && q > y && Element[{EE}, Reals]] &

int = Integrate[a[q, EE, K, y], {q, y, \[Infinity]}, 
        Assumptions -> 0 < K < y < 1 && q > y && Element[{EE}, Reals]]

{*   ConditionalExpression[(1/(
   4 K))(K Sqrt[-4 EE - K^2] \[Pi] + 
  K Sqrt[-4 EE - 
 K^2] (ArcTan[(K - 2 y)/Sqrt[-4 EE - K^2]] - 
  ArcTan[(K + 2 y)/Sqrt[-4 EE - K^2]]) + (2 EE + K^2) ArcTanh[(
 K y)/(EE - y^2)] - 
 y (2 K + y Log[1 + (2 K y)/(EE - y (K + y))])), 4 EE + K^2 < 0] 
     *}

.

ff[EE_, K_, y_] = int[[1]];

rhs[x_?NumericQ, y_?NumericQ, Epol_?NumericQ] := 
   4/Pi*NIntegrate[
   K^2/(x - 2/Pi*y + K^2/4 + Epol - 2/Pi*ff[Epol, K, y]), {K, 0, y}, 
  MaxRecursion -> 50]

EEfr[x_, y_] := EE /. First@FindRoot[EE == rhs[x, y, EE], {EE, -.1}]

Here a plot of Epol for negative x. for positive x there are more convergence problems may be due to low working precision. Further the condition with analytical integral 4 EE + K^2 < 0 has to be regarded.

pl = Plot3D[EEfr[x, y], {x, -5, 0}, {y, 0, 1}, PlotRange -> All, PlotPoints -> 10]

NIntegral failing to converge

2 Answers2