# The problem when numerical integrate (NIntegrate) a large order Bessel function

I am trying to solve this problem by using the NIntegrate function to evaluate the following code

cm = 1000*4200;
Q = 5*10^-4;
Kr = 0.64;
Kz = 0.64;
rw = 0.075;
Tin = 35 - 13.26;
n = 0.3;
b = 50;
K1 = 2.5; K2 = 2.5; b1 = 50; b2 = 50; \
\[Theta] = 5.5; \[Lambda]s = 1.6; rs = 10^-2; cs = 2000*800; c1 =
c2 = cm; g = 0.015; g1 = 0.025; g2 = 0.035; Tm0 = 13.26; RR = 10;
TMsa[r_] :=
Tm0 + r^((cm Q)/(4 b Kr n \[Pi]))*
NIntegrate[(-((2 (-2 (-1 + E^((a b Sqrt[Kr])/Sqrt[
Kz]))^2 K1 K2 Sqrt[Kz] +
a^2 b (1 + E^((2 a b Sqrt[Kr])/Sqrt[Kz])) (b2 K1 +
b1 K2) Kr Sqrt[Kz] +
a^3 b b1 b2 (-1 + E^((2 a b Sqrt[Kr])/Sqrt[Kz])) Kr^(3/2)
Kz + a (-1 + E^((2 a b Sqrt[Kr])/Sqrt[Kz])) Sqrt[
Kr] (b K1 K2 - (b2 K1 + b1 K2) Kz)) rw^(-((cm Q)/(
4 b Kr n \[Pi]))) Tin)/(a^3 b Sqrt[
Kr] (K1 ((-1 + E^((2 a b Sqrt[Kr])/Sqrt[Kz])) K2 +
a b2 (1 + E^((2 a b Sqrt[Kr])/Sqrt[Kz])) Sqrt[Kr]
Sqrt[Kz]) +
a b1 Sqrt[
Kr] ((1 + E^((2 a b Sqrt[Kr])/Sqrt[Kz])) K2 +
a b2 (-1 + E^((2 a b Sqrt[Kr])/Sqrt[Kz])) Sqrt[Kr]
Sqrt[Kz]) Sqrt[Kz]) \[Pi])))*((-BesselJ[(cm Q)/(
4 b Kr n \[Pi]), a rw] BesselY[(cm Q)/(4 b Kr n \[Pi]),
a r] + BesselJ[(cm Q)/(4 b Kr n \[Pi]), a r] BesselY[(
cm Q)/(4 b Kr n \[Pi]), a rw])/(BesselJ[(cm Q)/(
4 b Kr n \[Pi]), a rw]^2 +
BesselY[(cm Q)/(4 b Kr n \[Pi]), a rw]^2))*a, {a,
0, \[Infinity]},
Method -> {"LevinRule", "LevinFunctions" -> {BesselJ, BesselY},
"MaxOrder" -> 100}]


I found that as the order in the Bessel function (i.e., (cm Q)/(4 b Kr n \[Pi])) is large, it will show the following message

TMsa[3]
NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in a near {a} = {127.}. NIntegrate obtained -3.42963*10^9 and 1.658176851051496*^13 for the integral and error estimates. >>
-6.93058*10^17


The result is totally wrong. I still don't figure out where the trouble is happened. I were also trying to increase the WorkingPrecision or MaxRecursion but still get the wrong answer. Thank you for your time.

• "trying to increase the WorkingPrecision" does not work when you have numbers like 0.64 or 13.26 in your parameters. Convert them to exact rational fractions like 16/25 or 663/50 and try again. – J. M. is away Jul 27 '17 at 11:39
• Thank you. I just tried to convert the all value to rational fraction form but it still not work. – JOJO Jul 27 '17 at 11:46
• How large is the order? Have you tried LevinIntegrandReduce[]` and considered its output? (Sorry, I have no M this week. I remember this or a similar integral, but I forget how it worked out.) – Michael E2 Jul 27 '17 at 14:04
• The order is around 17.4. The output inform for LevinIntegrandReduce[] is so strange. I will study about that. – JOJO Jul 27 '17 at 15:38