0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ "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. $\endgroup$ – J. M. is away Jul 27 '17 at 11:39
  • $\begingroup$ Thank you. I just tried to convert the all value to rational fraction form but it still not work. $\endgroup$ – JOJO Jul 27 '17 at 11:46
  • $\begingroup$ 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.) $\endgroup$ – Michael E2 Jul 27 '17 at 14:04
  • $\begingroup$ The order is around 17.4. The output inform for LevinIntegrandReduce[] is so strange. I will study about that. $\endgroup$ – JOJO Jul 27 '17 at 15:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.