3
$\begingroup$

Why two the same integrals give different values?

Psi[r_, n_] := (-1)^n *Exp[-1/2 *r^2] *Sqrt[2 *n!/Gamma[n + 3/2]]* 
   LaguerreL[n, 1/2, r^2];
KA[n1_, n2_] := 
 0.0032433*Integrate[Psi[r, n2]*r^2*Psi[r, n1]*r^2, {r, 0, Infinity}]
KA[87, 12]

Out[1356]= 0.

Psi[r_, n_] := (-1)^n *Exp[-1/2 *r^2] *Sqrt[2 *n!/Gamma[n + 3/2]]* 
   LaguerreL[n, 1/2, r^2];
KA[n1_, n2_] := 
 Integrate[Psi[r, n2]*0.0032433*r^2*Psi[r, n1]*r^2, {r, 0, Infinity}]
KA[87, 12]

Out[1352]= 6.14098*10^21
```
Improve this question
$\endgroup$
5
  • 3
    $\begingroup$ What happens if you use the exact rational 32433/10000000? $\endgroup$
    – Ghoster
    Mar 1 at 5:59
  • $\begingroup$ @Ghoster, thanks! $\endgroup$
    – Mam Mam
    Mar 1 at 6:08
  • 3
    $\begingroup$ The first integral is evaluated symbolically. In the second integral, there are floating point numbers in the integrand. Then Mathematica switches silently to NIntegrate, which uses a numerical approximation. Not sure why it fails so drastically. It is very oscillatory, but even on intervals far to the right I get absurdly large integrals. You might want to consider to send this example to Wolfram Support. $\endgroup$
    – Henrik Schumacher
    Mar 1 at 6:08
  • 1
    $\begingroup$ For such cases I'd highly recommend using FindSequenceFunction to discover exact expressions. In your case, ka[n_, n_] = 2 n + 3/2; ka[n1_, n2_] /; Abs[n1 - n2] == 1 = Sqrt[(n1 + n2 + 1) (n1 + n2 + 2)]/2; ka[_, _] = 0; (I've left out the numerical prefactor of $0.0032433$). I know this wasn't your question; but very generally, finding such regularities is more useful than actually integrating. $\endgroup$
    – Roman
    Mar 1 at 9:22
  • $\begingroup$ @MamMam you got effectively a good answer on the comments. If the people commenting do not offer an actual answer post outside the comments section, you should feel free to transform the comments into your own answer. Please add the linking details and outputs that explain the solution. This way, you will learn more from this problem, your question and answer will be more useful for future visitors and you will gain more reputation on the site. Win win win. You could also add the response from Wolfram Support if you get one. $\endgroup$
    – rhermans
    Mar 1 at 9:43

1 Answer 1

Reset to default
4
$\begingroup$

Evaluating high-order polynomials is numerically unstable.

Example: exact evaluation followed by numericalization is stable,

Psi[7, 87] // N
(*    -0.0271578    *)

Calling LaguerreL with numerical parameters is also stable,

Psi[7., 87]
(*    -0.0271578    *)

Expanding the Laguerre polynomial and inserting numbers, however, is unstable:

Psi[r, 87] /. r -> 7.
(*    5.64566*10^18    *)

The reason is that LaguerreL[87, 1/2, r^2] is a polynomial of order 174 in r, with (for $r=7$) contributions of order $10^{45}$ cancelling each other out to give a result of order $10^9$. This level of numerical cancellation is too much for 64-bit double-precision numbers to handle.

Illustrating this problem with a some graphics:

Plot[Psi[r, 87], {r, 0, 10}, PlotRange -> All]

enter image description here

Plot[Psi[r, 87] // Evaluate, {r, 0, 10}, PlotRange -> All]

enter image description here

Plot[Psi[r, 87] // Evaluate, {r, 0, 10}, PlotRange -> All,
     WorkingPrecision -> 100]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ thank you very much for these useful comments! $\endgroup$
    – Mam Mam
    Mar 3 at 14:56
  • $\begingroup$ You’re welcome. If you think that this (or another) answer is good enough, please consider using the check mark button ✓ to guide future visitors of this question. $\endgroup$
    – Roman
    Mar 4 at 10:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.