Message Boards Message Boards

How to speed up NIntegrate for an integral?

Posted 2 years ago

Hello, I am trying to evaluate numerically an integral of a fairly complicated function. Here is my code:

K1K0[x_]=BesselK[1,x]/BesselK[0,x];
J0J0Y0Y0[z_]:=BesselJ[0,z]*BesselJ[0,z]+BesselY[0,z]*BesselY[0,z];
J1J1Y1Y1[z_]:=BesselJ[1,z]*BesselJ[1,z]+BesselY[1,z]*BesselY[1,z];
J0J1Y0Y1[z_]:=BesselJ[0,z]*BesselJ[1,z]+BesselY[0,z]*BesselY[1,z];
G1[w_,kap_]:=(kap-w)*((K1K0[Sqrt[kap-w]])^2)*J0J0Y0Y0[Sqrt[w]];
G2[w_,kap_]:=Sqrt[w*(kap-w)]*K1K0[Sqrt[kap-w]]*J0J1Y0Y1[Sqrt[w]];
G3[w_,kap_]:=w*J1J1Y1Y1[Sqrt[w]];
BESSELTERM[w_,kap_]:=((kap-w)/w)*((K1K0[Sqrt[kap-w]])^2)/(G1[w,kap]+2*th*G2[w,kap]+th*th*G3[w,kap]);
u = 4;
th = Exp[u];
t=1;
k=1000;
result=(2/π^2)*NIntegrate[BESSELTERM[1/y,k]*Exp[-t/y]/y^2,{y,1/k,∞},Method->"GaussKronrodRule",MaxRecursion->40,WorkingPrecision->16];
CForm[N[result,30]]

The result of running this code is obtained in a few seconds and it is: 0.3924117267170341. However, the result is not very accurate. I need about 30-35 accurate digits. But if I set WorkingPrecision->35, then the code slows down enormously, and the result is obtained in about half an hour. In addition, I need to vary the parameter u up to u=50, and I observe that the computational time rapidly increases with u for u > 4. Parameters t and k may also vary in the intervals -6 <=log10(t)<=6 and -3<=log10(k)<=3, which often causes error messages from MATHEMATICA, to avoid which WorkingPrecision has to be increased, which seems to further increase the cost. Is there any way to speed up these calculations, perhaps by chosing other options of NIntegrate? I don't understand why these calculations are so slow. Function BESSELTERM[w,kap] is singular t w=0. Lesław

POSTED BY: Leslaw Bieniasz
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract