The graph of my 1-D function is well-behaved, but NIntegrate gives impossible results and lots of errors

Question

I have a one-dimensional integral that should converge, but Mathematica can't seem to handle using the options I have tried. Here's the setup:

The wavefunctions for a Morse Oscillator are given by:

psi[v_, r_, m_, De_, a_, re_] = (((v! Sqrt[a (2 lambda - 2 v - 1)/
        (Gamma[v + 1] Gamma[2 lambda - v])]) 
        z^(lambda - v - 1/2) E^(-z/2) LaguerreL[v, 2 lambda - 2 v - 1, z])
        /. z -> (2 lambda E^(-a (r - re)))) 
        /. lambda -> Sqrt[2 m De]/(a hb);

Here:

• hb is a universal constant (a positive real number)
• m is a positive real constant that is the same for all states
• De, a, and re are positive real numbers that are constant for a given state
• v is whole number that specifies a particular state
• r is the variable to integrate over, from -infinity to infinity

Typical values are:

mass = 1.2 10^(-26);
hb = 2.37 10^(-13);
De1 = 52000;
De2 = 12000;
a1 = 2.30;
a2 = 2.17;
re1 = 1.21;
re2 = 1.60;

What I need to do is integrate the product of two versions of psi, one with the "1" parameters above and the other with the "2" parameters above; the first one will have a series of different values of v ranging from 0 to 32 (for this set of parameters), and the second will always have v=0. We can look at the graphs of the functions and their product this way:

Manipulate[Grid[{{Plot[{psi[v, r, mass, De1, a1, re1], 
                        psi[0, r, mass, De2, a2, re2]}, 
                        {r, 0, 3.5}, PlotRange -> Full, ImageSize -> Large], 
                  Plot[psi[v, r, mass, De1, a1, re1] psi[0, r, mass, De2, a2, re2], 
                        {r, 1, 2.2}, PlotRange -> Full, ImageSize -> Large]}}], 
           {v, 0, 32, 1}]

Here is what it looks like for v=32, the most oscillatory of the functions, and as you can see the function is nonetheless well-behaved. So I try to execute this for various values of v:

integration[v_] := NIntegrate[psi[v, r, mass, De1, a1, re1] 
                              psi[0, r, mass, De2, a2, re2], 
                              {r, -Infinity, Infinity}]

For this set of parameters, the integration works fine up through v=14. Above that, I start to get failure to converge errors and then later that integration is converging too slowly. The result correctly increases from v=0 to v=13 to a maximum of about 0.39, and then decreases, but it should approach zero as v continues to increase. Instead the result blows up, first negative (to -18) and then positive (to over 80,000).

After monitoring the progress using Sow/Reap, I tried a bunch of things that didn't help. I tried increasing MaxRecursion to 30 and WorkingPrecision to 60, and most of the error messages went away (ignoring the warning that the precision of the argument function isn't that high), but the answers that come out are still wrong in the same way. I have tried changing the integration method to LocalAdaptive, DoubleExponential, Trapezoidal, MonteCarlo, QuasiMonteCarlo, DoubleExponentialOscillatory, and ExtrapolatingOscillatory to no effect.

Most strikingly, I tried limiting the integration range to -10 < r < 10, and I still get this behavior. We can see visually from the Manipulate graph above that for v=32, the function is quite well-behaved over this region with a magnitude that never exceeds ±3, which means that it shouldn't be possible for the integral to be larger than 3*20=60. And yet NIntegrate gives me over 177,000.

I am running out of ideas. Help?

Answer 1

These tweaks give a result without errors:

integration[v_] := 
 NIntegrate[
  psi[v, r, mass, De1, a1, re1] psi[0, r, mass, De2, a2, 
    re2], {r, -Infinity, Infinity}, MaxRecursion -> 20, 
  WorkingPrecision -> 32, PrecisionGoal -> 8]

Block[{mass = 1.2`100 10^(-26),
  hb = 2.37`100 10^(-13),
  De1 = 52000,
  De2 = 12000,
  a1 = 2.30`100,
  a2 = 2.17`100,
  re1 = 1.21`100,
  re2 = 1.60`100},
 integration[32]
 ]

(*  1.3768685099994004181924204909584*10^-8  *)