Unstable NIntegrate result vs plot

Question

Here I have a $\Psi$ function with a few arguments:

\[CapitalPsi][r1_, r2_, Nsize_, c_, eta_] := \[CapitalPsi][r2, r1, Nsize, c, eta] = \!\( \*UnderoverscriptBox[\(\[Sum]\), \(n2 = 1\), \(Nsize\)]\(( \*UnderoverscriptBox[\(\[Sum]\), \(n1 = 
        1\), \(n2\)]\(c[\([n1]\)]\)[\([n2]\)]  \*SqrtBox[ FractionBox[\(64  \*SuperscriptBox[\((eta)\), \(6\)]\), \(n1 \((n1 + 1)\) n2 \((n2 + 
            1)\)\)]]  \*SuperscriptBox[\(E\), \(\(-eta\)*\((r1 + r2)\)\)] \((LaguerreL[
           n1 - 1, 2, 2*eta*r1] LaguerreL[n2 - 1, 2, 2*eta*r2] + 
         LaguerreL[n2 - 1, 2, 2*eta*r1] LaguerreL[n1 - 1, 2, 
           2*eta*r2])\))\)\);

$c$ is the coefficient that is stored as a 47x47 matrix and can be download here

Here comes the interesting thing:

First get rid of the dummy dimension

c = Import["~/Downloads/c.mat"];
c = ArrayReshape[c, Dimensions[c]~DeleteCases~1];

Then plot vs integrate:

Plot[\[CapitalPsi][r1, 0.01, 47, c, 0.675314001427934]^2*r1*r1, {r1,40,50}]

NIntegrate[\[CapitalPsi][r1, 0.01, 47, c, 0.675314001427934]^2*r1*
  r1, {r1, 40, 50}]

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in r1 near {r1} = {44.8631}. NIntegrate obtained 2420.258935677104and 347.10349104798587 for the integral and error estimates. >>

= 2420.26

But how could it be? The plot seems smooth but it seems there is a divergence in the integral that the singular point at 44.8631 contribute to the majority of the integral (2420.258935677104 out of 2420.26)

Update:

To make it easier for people to get the coefficient matrix $c$, you can also do it in this way:

URLSave["https://drive.google.com/uc?export=download&id=\1lZeEfZ1I2UmVaHJ_hRAZqtgkipI8P6su", "c.mat"]
c = Import["c.mat"]
c = ArrayReshape[c, Dimensions[c]~DeleteCases~1];

Answer 1

This is an interesting question, even though it has the feel of an edge-case. The problem arises because the accuracy of the function evaluation depends on whether the function

c = ToExpression@Import["https://pastebin.com/raw/jMWZasFN", "Text"]; (* another source *)
Ψ[r1, 1/100, 47, SetPrecision[c, 100], 0.675314001427934`100]^2*r1*r1

is evaluated with a symbolic r1 or a numeric r1. With a numeric r1 only a few digits are lost due to rounding error. However if the function is evaluated with a symbolic r1 and then a numeric value is substituted for r1, more digits are lost than one has at machine precision (or double precision, which is about 16 digits).

Plot evaluates its function after r1 has been assigned a value (effectively with Block according to the docs). Here we evaluate the function with a symbolic r1 before passing it to plot, which results in a catastrophic loss of accuracy. The plot is very different than the OP's.

Block[{f},
 f[r1_?NumericQ] = Ψ[r1, 0.01, 47, c, 0.675314001427934]^2*r1*r1;
 Plot[f[r1], {r1, 40, 50}, MaxRecursion -> 3]
 ]

NIntegrate on the other hand evaluates the integrand before constructing a NumericFunction to use in computing the integral, which is disastrous in the OP's case. Below we prevent the symbolic evaluation of the integrand with ?NumericQ protection and SetDelayed. NIntegrate evaluates accurately without problems.

Block[{f},
 f[r1_?NumericQ] := Ψ[r1, 0.01, 47, c, 0.675314001427934]^2*r1*r1;
 NIntegrate[f[r1], {r1, 40, 50}]
 ]
(*  0.0338579  *)

We can examine the precision loss by computing the function two ways with 100-digit arbitrary-precision numbers and calculating the precision loss. Evaluating with symbolic r1 results in almost 22 digits lost, much more than machine precision. The results are completely useless. However, evaluating with a numeric r1 results in under three digits lost, which leaves around 13 accurate digits, sufficient to compute a numerical integral to 8-digit accuracy.

(* evaluate with symbolic r1 first *)
(Ψ[r1, 1/100, 47, SetPrecision[c, 100], 0.675314001427934`100]^2*r1*r1 /. r1 -> 45.`100 // 
   Precision) - 100
(*  -21.9455  *)

(* assign r1 a numeric value before evaluation *)
Block[{r1 = 45.`100},
  Ψ[r1, 1/100, 47, SetPrecision[c, 100], 0.675314001427934`100]^2*r1*r1 // Precision
  ] - 100
(*  -2.75681  *)

Answer 2

Looking at your function to be integrated, it has multiplications of coefficients of high value by others of very, very low one. Arithmetic problems could arise. Maybe some expert could add light to this.

I propose a workaround based on FunctionInterpolation of your integrand in the range of interest:

if = FunctionInterpolation[\[CapitalPsi][r1, 0.01, 47, c, 0.675314001427934]^2*r1*r1, 
     {r1, 40, 50}];
Plot[\[CapitalPsi][r1, 0.01, 47, c, 0.675314001427934]^2*r1*r1 - if[r1], {r1, 40, 50}]

The difference with the original is low. So integrate it:

NIntegrate[if[r1], {r1, 40, 50}]

(* 0.0338579 *)

which seems to be a reasonable result.