# Unstable NIntegrate result vs plot

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];

• We could better help you if you posted code the we could evaluate. That is, if you gave use an actual value for c, not an expression we can't evaluate. – m_goldberg Feb 9 '18 at 4:46
• @m_goldberg The c value can be downloaded from the link in the question. Maybe I should have made it more noticeable. – James Feb 9 '18 at 5:07
• It would be nice if there were included an Import[] command that downloaded the data. Not sure if ufile.io supports that. Here are some discussions on Meta about ways to make large amounts of data or code available on SE: (1351), (1520), (2145). (I like pastebin, because I can inspect the code easily in a browser before deciding whether to download it into Mathematica.) – Michael E2 Feb 9 '18 at 14:30
• I used Google Drive so it will never expire, which I think is better than other choices. Please let me know if you have any trouble downloading the data. – James Feb 9 '18 at 20:53
• The URLSave did not work for me, so I copied the data here: c = ToExpression@Import["https://pastebin.com/raw/jMWZasFN", "Text"];, which also doesn't expire. – Michael E2 Feb 10 '18 at 13:56

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.675314001427934100]^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.675314001427934100]^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.675314001427934100]^2*r1*r1 // Precision
] - 100
(*  -2.75681  *)

• Happy to learn from your insight +1 – José Antonio Díaz Navas Feb 10 '18 at 11:19
• @JoséAntonioDíazNavas Your function interpolation made me wonder how the value of NIntegrate could become so large, so +1 to you, too. (It turns out that FunctionInterpolation evaluates its argument the same way as Plot.) – Michael E2 Feb 10 '18 at 13:39

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.

• Thanks! This is amazing because I assumed interpolation would not be a good approximation at all, but I was wrong! – James Feb 9 '18 at 20:05
• Well, your function is quite smooth in the range of integration, so why not to make a try? :)) – José Antonio Díaz Navas Feb 9 '18 at 20:07