1
$\begingroup$

How can I get this section of the code to solve this integration faster. I have rather large values but I don't need it to be too specific. However if I lower the precision it says the integration is 0. The function is only in a very small section of the x values according to the 3d plot. However if I analyze the smaller section it still takes a considerable amount of time. Any help is appreciated

timewj2 = AbsoluteTiming[
   Wjlambda5 = (1/8.952*10^8)*
      NIntegrate[
       Dot[Re[{0, -8.952*10^8*(15.675459922348419449)*(((-8.\
9652950534407867999303353030387276032035253437793*10^768 - 
                   2.1141452007443239580927869423477693294460547800052\
7*10^769*I)*
                 BesselJ[
                  1, (75062.4870217452581789128 + 
                    75062.4983440564559965704*I )*
                   x] + (2.1141452007443239580927869423477693294*10^\
769 - 8.965295053440786799930335303038727603*10^768*I)*
                 BesselY[
                  1, (75062.4870217452581789128 + 
                    75062.4983440564559965704*I )*x])*
              Exp[I*(1)*(41.2281675906100480252)*s]), 0}], 
         Re[{0, -8.952*10^8*(15.675459922348419449)*(((-8.\
9652950534407867999303353030387276032035253437793*10^768 - 
                   2.1141452007443239580927869423477693294460547800052\
7*10^769*I)*
                 BesselJ[
                  1, (75062.4870217452581789128 + 
                    75062.4983440564559965704*I )*
                   x] + (2.1141452007443239580927869423477693294*10^\
769 - 8.965295053440786799930335303038727603*10^768*I)*
                 BesselY[
                  1, (75062.4870217452581789128 + 
                    75062.4983440564559965704*I )*x])*
              Exp[I*(1)*(41.2281675906100480252)*s]), 0}]]*
        x, {x, (0.02357124714249428777766010023597686995), \
(0.04022197675031113473009371147375231401)}, {s, 
        0, (0.15240030480060962059241091992589645088)}, 
       WorkingPrecision -> (2000 - 1), MaxPoints -> 80];
   ];
Print["Time Wj5 (s) = ", N[timewj2[[1]],10]];

EDIT: Like This?

prec = 100;
σ5 = SetPrecision[500*(1.492*^6), prec];
vs5 = SetPrecision[15.675459, prec];
γ5 = SetPrecision[75062.487021745, prec];
e = SetPrecision[-8.965295053440786799*10^768 - 
    2.1141452007443239588*10^769*I, prec];
j = SetPrecision[
   2.1141452007443239580927*10^769 - 8.96529505344078679*10^768*I, 
   prec];
mm = SetPrecision[1, prec];
k = SetPrecision[41.2281675, prec];
r45 = SetPrecision[0.02357124714249428777766010, prec];
r56 = SetPrecision[0.0402219767503111347300937114, prec];
lam = SetPrecision[.1524, prec];
timewj2 = AbsoluteTiming[
   Wjl5 = (1/σ5)*
      NIntegrate[
       Dot[
         Re[{0, -σ5* vs5*((e*BesselJ[1, γ5*x] + j*BesselY[1, γ5*x])*Exp[I*mm*k*s]), 0}], 
         Re[{0, -σ5* vs5*((e*BesselJ[1, γ5*x] + j*BesselY[1, γ5*x])*Exp[I*mm*k*s]), 0}]
         ]*x, 
       {x, r45, r56}, 
       {s, 0, lam}, 
       WorkingPrecision -> prec - 1, 
       MaxPoints -> 80
       ];
   ];
Wjl5
Print["Time Wj5 (s) = ", N[timewj2, prec2]];

But my output is still wrong its getting 1.02*10^1542 and the answer is close to 0 but not 0.

Improve this question
$\endgroup$
13
  • 1
    $\begingroup$ You seem to mixing machine arithmetic with arbitrary precision arithmetic. That never works well. I suggest expressing all quantities at a modest specific precision and working with that. I also suggest exercising patience. Difficult computations often take some time. $\endgroup$
    – m_goldberg
    May 16, 2019 at 15:37
  • $\begingroup$ This is the shortest version of the numerical equivalent i could find. I set the precision of all calculations using SetPrecision and right now it is at 2000. The imaginary numbers being multiplied by the bessel function are from a resulting simultaneous equation solve that requires that level of precision. Furthermore the 3d plot will not show an accurate version unless set to a precision of say 1500. This is just one equation needed to solve for my calculations but it takes up the most time. All the other equations are done in 14 seconds or so. This one is 11405.1. Any way I can chip it down $\endgroup$
    – Rookey
    May 16, 2019 at 15:48
  • 1
    $\begingroup$ When you include machine numbers like 8.952*10^8 in you code, then all computations involving that number in any way are done at machine precision, what you set with SetPrecision notwithstanding. I advise doing what I suggested in my 1st comment: re-expressing all the arithmetic quantities at a modest precision, say 20, and working strictly at that precision in your calculations. Also, I don't know what you are referring to when you mention "3rd plot". You show no plots. Further, plotting doesn't require high precision; it's always done at machine precision. $\endgroup$
    – m_goldberg
    May 16, 2019 at 16:11
  • 1
    $\begingroup$ Moreover I observe that the integrand can be brought into the form $f_1(x) \cdot g_1(s) + f_2(x) \cdot g_2(s) $, so for the integral, we obtain $(\int f_1 (x) \, \mathrm{d} x) \, (\int g_1 (s) \, \mathrm{d} s) + (\int f_2 (x) \, \mathrm{d} x) \, (\int g_2 (s) \, \mathrm{d} s) $. And the 1D-integrals should be much easier to compute. $\endgroup$
    – Henrik Schumacher
    May 16, 2019 at 17:03
  • 2
    $\begingroup$ Dot[Re[{0,-σ5*vs5*((e*BesselJ[1,γ5*x]+j*BesselY[1,γ5*x])*Exp[I*mm*k*s]),0}], Re[{0,-σ5*vs5*((e*BesselJ[1,γ5*x]+j*BesselY[1,γ5*x])*Exp[I*mm*k*s]),0}]]*x == x*Re[-σ5*vs5*((e*BesselJ[1,γ5*x]+j*BesselY[1,γ5*x])*Exp[I*mm*k*s])]^2 and that might make your code smaller and a little faster. Check to make certain this is correct. $\endgroup$
    – Bill
    May 16, 2019 at 18:07

1 Answer 1

Reset to default
1
$\begingroup$

We introduce the notation

f = a*BesselJ[1, \[Gamma]5*x] + b*BesselY[1, \[Gamma]5*x];
g = Exp[I*k*s];

then the required integral is

Wjl5 = (1/\[Sigma]5)*
      NIntegrate[
       Dot[Re[{0, f*g, 0}], Re[{0, f*g, 0}]]*x, {x, r45, r56}, {s, 0, 
        lam}];

The multiplication operation Dot can be performed, 

Wjl5 = (1/\[Sigma]5)*
          NIntegrate[
           Re[f*g]^2*x, {x, r45, r56}, {s, 0, 
            lam}];

The integral over s is calculated analytically using Re[fg]^2*x = (Re[f]*Cos[k*s] - Im[f]*Sin[k*s])^2*x, we have

Integrate[(Re[f]*Cos[k*s] - Im[f]*Sin[k*s])^2, {s, 0, lam}]

(*Out[]= (4 f k lam Conjugate[f] - 
 8 Im[f] Re[f] Sin[k lam]^2 + (f^2 + Conjugate[f]^2) Sin[
   2 k lam])/(8 k)*)

define the parameters

\[Sigma]5 = 500*(1.492*^6);
vs5 = 15.675459;
\[Gamma]5 = 75062.4870217452581789128 + 75062.4983440564559965704*I;
e = -8.965295053440786799*10^768 - 2.1141452007443239588*10^769*I;
j = 2.1141452007443239580927*10^769 - 8.96529505344078679*10^768*I;
mm = 1;
k = 41.2281675;
r45 = 0.02357124714249428777766010;
r56 = 0.0402219767503111347300937114;
lam = .1524;
a = -\[Sigma]5*vs5*e;
b = -\[Sigma]5*vs5*j;

finally calculate the integral

With[{f = 
    a*BesselJ[1, \[Gamma]5*x] + b*BesselY[1, \[Gamma]5*x]}, 
  NIntegrate[(
    4 f k lam Conjugate[f] - 
     8 Im[f] Re[f] Sin[k lam]^2 + (f^2 + Conjugate[f]^2) Sin[
       2 k lam])/(8 k)*x, {x, r45, r56}, 
   Method -> "MonteCarlo"]] // AbsoluteTiming

(*During evaluation of In[]:= NIntegrate::maxp: The integral failed to converge after 50100 integrand evaluations. NIntegrate obtained 0. +0. I and 0.` for the integral and error estimates.*)

Out[]= {5.69314, 0. + 0. I}
Improve this answer
$\endgroup$

Your Answer

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

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