1

Why does the error occur for P[34, 14] whereas for P[33, 14] has no error? How to get rid of this error?

ClearAll["Global`*"]

Psi[r_, n_] := (-1)^n Exp[-1/2 r^2] Sqrt[2 n!/Gamma[n + 3/2]]* 
   LaguerreL[n, 1/2, r^2];


P[n1_, n2_] := 
  NIntegrate[
   Psi[r, n2]4.915/rExp[-r*1.65]Psi[r, n1]r^2, {r, 
    0, [Infinity]}];


P[33, 14]


Out[1008]= -0.500849


P[34, 14]


During evaluation of In[1009]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in r near {r} = {6.42115}. NIntegrate obtained 0.488132829059508and 9.895807414830366*^-7 for the integral and error estimates.


Out[1009]= 0.488133
Mam Mam
  • 1,843
  • 2
  • 9
  • I am not getting any errors when I run your code. v12 on a mac – bmf Dec 18 '22 at 15:01
  • @bmf, my version: "12.0.0 for Microsoft Windows (64-bit) (April 6, 2019)". Try to run P[41,63] – Mam Mam Dec 18 '22 at 15:06
  • see the screenshot for the value you quoted in the comment. Perhaps you should add the details of the operational system and the verion in the main question – bmf Dec 18 '22 at 15:08
  • @bmf, strange situation – Mam Mam Dec 18 '22 at 15:33
  • @user293787, could you please take a look at this question https://mathematica.stackexchange.com/questions/277927/%d0%a1an-this-integral-be-solved-exactly-using-integrate – Mam Mam Dec 29 '22 at 14:01

1 Answers1

5

OP is trying to numerically evaluate highly oscillating integrals. In this case, one can do the integration exactly, using high precision only at the end to add terms:

f0[a_,n_]=Integrate[r^n*Exp[-r^2]*Exp[-a*r],{r,0,Infinity},Assumptions->{a>0,n>0}];
f[a_,n_]:=f[a,n]=N[f0[a,n],100 (* increase if necessary *)];

coeff[n_]:=(-1)^nSqrt[2n!/Gamma[n+3/2]];


Px[n1_,n2_]:=Px[n1,n2]=If[n1>n2,Px[n2,n1],
  coeff[n1]coeff[n2]4915/1000Total[Map[f[165/100,#[[1,1]]+1]#[[2]]&,
  CoefficientRules[LaguerreL[n1,1/2,r^2]*LaguerreL[n2,1/2,r^2],r]]]//N];

It is also quite fast:

AbsoluteTiming[Table[Px[n1,n2],{n1,0,50},{n2,0,50}];]
(* about 3 seconds *)
user293787
  • 11,833
  • 10
  • 28
  • Thanks! Could you please describe this part in more detail, it's a bit difficult to understand Px[n1_,n2_]:=Px[n1,n2]=If[n1>n2,Px[n2,n1], coeff[n1]coeff[n2]4915/1000Total[Map[f[165/100,#[[1,1]]+1]#[[2]]&, CoefficientRules[LaguerreL[n1,1/2,r^2]*LaguerreL[n2,1/2,r^2],r]]]//N]; – Mam Mam Dec 24 '22 at 22:19
  • The Px[n1_,n2_]:=Px[n1,n2]= part is called memoization, see this tutorial. If you are asking about Total[Map[f[165/100,#[[1,1]]+1]*#[[2]]&,CoefficientRules[...,r]]]//N I propose you run the pieces and look at intermediate results, for example step1 = CoefficientRules[13*r^2-7*r+2,r] then step2 = Map[ftemp[165/100,#[[1,1]]+1]*#[[2]]&,step1] then step3 = step2/.ftemp->f then step4=Total[step3] then N[step4]. If there is a step you do not understand, let me know which step. – user293787 Dec 25 '22 at 05:43
  • The polynomial 13*r^2-7*r+2 in step1 is just an example, in the actual code this is replaced by the product of Laguerre polynomials. – user293787 Dec 25 '22 at 05:45
  • Thanks, this part Exp[-r^2] in f0[a_,n_] is this part Exp[-1/2 *r^2] in Psi[r_, n_]? I don't understand why there is no 0.5 – Mam Mam Dec 25 '22 at 23:10
  • In your integral there is a factor Exp[-1/2*r^2] from Psi[r,n1] and another factor Exp[-1/2*r^2] from Psi[r,n2], that gives Exp[-1/2*r^2]*Exp[-1/2*r^2] == Exp[-r^2]. – user293787 Dec 26 '22 at 04:59
  • Thank you very much! – Mam Mam Dec 26 '22 at 08:57