I have a function F which behaves like you can see on the plot shown below. I need to calculate

NIntegrate[F[r], {r, 0, Infinity}]

but it returns a huge result, something like 10^80. If I apply Method -> "LevinRule", it still returns an unreasonable number. How can I caluclate the integral?

Perhaps it will help if I tell you that the function is based on Laguerre polynomials.


Here is the code defining the function via some number of auxiliary functions and variables

Z = 88
n = 120
k = -2
m = 1
alpha = 7.2973525664*^-3
alpha = alpha*Z
nr = n - Abs[k]
gam = Sqrt[k^2 - (alpha)^2]
enk = m/Sqrt[1 + (alpha)^2/(nr + gam)^2]
cnk = Sqrt[m^2 - enk^2]
ap = If[nr == 0, 0,
  ap = Sqrt[(cnk (nr - 1)!  / 
    Gamma[nr + 2 gam + 1])*
     ( (nr + gam + gam/k*Sqrt[alpha^2 + (nr + gam)^2])/
        (2 gam^2/k^2 (alpha^2 + (nr + gam)^2)))]]
am = Sqrt[(cnk*nr!  / 
  Gamma[nr + 2 gam])*
   ( (nr + gam - gam/k*Sqrt[alpha^2 + (nr + gam)^2])/
      (2 gam^2/k^2 (alpha^2 + (nr + gam)^2)))]
up[r_] := (2 cnk*r)^gam*Exp[-cnk*r]*ap*2 cnk*r*
  LaguerreL[nr - 1, 2 gam + 1, 2 cnk*r]
um[r_] := (2 cnk*r)^gam*Exp[-cnk*r]*am*
  LaguerreL[nr, 2 gam - 1, 2 cnk*r]
g[r_] := 1/Sqrt[2 k (k - gam)]*(alpha*up[r] - (k - gam) um[r])
f[r_] := 1/Sqrt[2 k (k - gam)]*(alpha*um[r] - (k - gam) up[r])
rho[r_] := Abs[g[r]]^2 + Abs[f[r]]^2
NIntegrate[rho[r], {r, 0, Infinity}, Method -> "LevinRule"]
Plot[rho[r], {r, 0, 55000}, PlotRange -> All]
  • 6
    $\begingroup$ It is impossible!, because we do not know the F[r] function.Give it to see what we can do. $\endgroup$ Dec 19, 2016 at 19:10
  • 2
    $\begingroup$ does the function end like that at 50,000 or keep oscillating? If it keeps growing it obviously diverges. $\endgroup$
    – george2079
    Dec 19, 2016 at 20:26
  • 2
    $\begingroup$ If the integration range is finite you can find the zeroes of function (or extrema points) and integrate between them. See the answers in this discussion for integration between a given set of points. $\endgroup$ Dec 19, 2016 at 20:50
  • 1
    $\begingroup$ Have you tried a higher WorkingPrecision? $\endgroup$
    – Greg Hurst
    Dec 20, 2016 at 3:29
  • $\begingroup$ Dear friends, thank you for your replies. As asked, I have attached the code defining the function. I tried also different WorkingPrecision and integration by small pieces (though not between zeros). No progress. $\endgroup$ Dec 20, 2016 at 9:31

2 Answers 2


A bit of experimentation indicates that "LevinRule" does not work well, and that "MaxErrorIncreases" and MaxRecursion must be increased. Most importantly, WorkingPrecision must be increased substantially. To do so, first increase the precision of the only finite precision constant.

alpha = SetPrecision[7.2973525664*^-3, 65];

Then, set the options of NIntegrate as follows.

NIntegrate[rho[r], {r, 0, 55000}, Method -> {"GlobalAdaptive", 
    "SymbolicProcessing" -> 0, "MaxErrorIncreases" -> 2000}, MaxRecursion -> 20, 
    WorkingPrecision -> 60, PrecisionGoal -> 8, AccuracyGoal -> Infinity]
(* 0.999999999862974654189196482270003235458935716234694791163667 *)

as desired. AbsoluteTiming for the computation is 1.4 seconds.

  • $\begingroup$ I think the "ClenshawCurtisOscillatoryRule" is a red herring. I get your result merely with the increased precision and recursion (and the restriction to the effective support of the integrand). $\endgroup$
    – Michael E2
    Dec 21, 2016 at 14:36
  • $\begingroup$ @MichaelE2 Actually, my claim is that "ClenshawCurtisOscillatoryRule" works better than "LevinRule", I chose the former over other rules I tried, because it was slightly faster, at least for n == 30, for which I ran a number of tests. By the way, nice answer (+1). $\endgroup$
    – bbgodfrey
    Dec 21, 2016 at 16:02
  • $\begingroup$ I see (& thanks, and +1 to you, too). The choice of the Levin rule, btw, is inappropriate because there is no oscillatory component of the integrand (high-degree polynomials don't count for it). In fact, Method -> "LevinRule" with your other changes gives an accurate result (albeit by rejecting Levin and using Gauss-Kronrod). I don't know much about "ClenshawCurtisOscillatoryRule" beyond the docs, which indicate that the integrand is not of the appropriate form for it either; but mainly I was just suggesting your answer could be simplified. $\endgroup$
    – Michael E2
    Dec 21, 2016 at 16:26
  • $\begingroup$ Interesting side note: Looking at what "ClenshawCurtisOscillatoryRule" actually does in this case, I see it is effectively equivalent to Method -> {"GaussKronrodRule", "Points" -> 20}. Not sure if that's normal or not. $\endgroup$
    – Michael E2
    Dec 21, 2016 at 16:43

The source of the problem

The problem is the expansion of LaguerreL[], with degree nr == 118. If it is evaluated with a symbolic argument r, it will expand to power series form (the usual linear combination of powers of r); in this form, it is impossible to evaluate accurately at machine precision. A very high precision is needed, and in addition, the machine precision coefficients, which derive from the definition of alpha, need to have their precision increased. I set alpha = 7.2973525664`100*^-3 for the following. For example, compare the plots:

Plot[rho[r], {r, 0, 55000}, PlotRange -> All, 
 PlotLabel -> "Evaluated on a numeric argument"]

Plot[Evaluate@rho[r], {r, 0, 55000}, PlotRange -> All, 
 PlotLabel -> "Evaluated on a symbolic argument, with low precision"]

Plot[Evaluate@Rationalize[rho[r], 0], {r, 0, 55000}, PlotRange -> All,
  WorkingPrecision -> 50, 
 PlotLabel -> "Evaluated on a symbolic argument, with insufficient precision"]

Plot[Evaluate@Rationalize[rho[r], 0], {r, 0, 55000}, PlotRange -> All, 
 WorkingPrecision -> 60, 
 PlotLabel -> "Evaluated on a symbolic argument, with high precision"]

Here we can see why WorkingPrecision -> 60 was successful in @bbgodfrey's answer.

An easy fix to the problem

The trick is to evaluate LaguerreL[] only on numeric arguments using ?NumericQ to control the evaluation. The built-in routine will evaluate LaguerreL[] accurately on machine precision inputs. This may be used closest to the source of the problem in up[] and um[] or simply on the top-level function rho[]:

(* with alpha = 7.2973525664`*^-3, machine precision as in the OP *)
rho[r_?NumericQ] := Abs[g[r]]^2 + Abs[f[r]]^2; 
NIntegrate[rho[r], {r, 0, Infinity}, MinRecursion -> 3, MaxRecursion -> 20] // AbsoluteTiming
(*  {0.727109, 1.0000000000001033`}   *)

I wouldn't consider this an oscillatory integral, since it eventually stops oscillating. (The Levin rule expects an oscillatory function like Sin[] or BesselJ[] etc.) However it does have many oscillation, so a sufficiently fine sampling is necessary, which is the purpose of MinRecursion -> 3 above. This could also be achieved with a high number of sample points, such as Method -> {"GaussKronrodRule", "Points" -> 61}.

Splitting the interval up can speed things up in this case, since the function effectively has support over a finite interval and is analytic. Something like this takes only 0.32 sec.:

NIntegrate[rho[r], {r, 0, 50000, Infinity}, Method -> {"GaussKronrodRule", "Points" -> 21}]

One final remark: The integral can be evaluated symbolically; but it takes a really long time, and the numerics issue of the expanded polynomial persists. Below is a relatively fast way, since the integrand, after removing the unnecessary Abs, expands into a linear combination of terms of the form Exp[-cnt r] * r^s for various powers s.

alpha = 7.2973525664`200*^-3;
(* ... re-execute OP's defs ...*)

i[s_] = Integrate[E^(-2 cnk r) r^s, {r, 0, Infinity}, Assumptions -> s > 0];

List @@ (rho[r] /. Abs -> Identity // Expand) /. E^(-2 cnk r)*r^s_ :> i[s] // Total

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.