0
$\begingroup$

I have a three dimensional integration, which takes about 30s using LocalAdaptive and "GaussKronrodRule". Using GlobalAdaptive and other rules will take more time. I have tried many different combinations of integration strategies, it seems that Method -> {"LocalAdaptive", Method -> "GaussKronrodRule"} is the fastest or some others will report error. I am wonder how to make this integration excute faster, the speed is very important to me. Below is the integration code:

Clear["Global`*"];
k = Sqrt[kx^2 + ky^2 + kz^2]; 
k0 = 1/100;  \[CapitalDelta]k = k0/10;  t0 = -1000;   v = 1;  z0 = 
 v t0;   kInfinity = 10 \[CapitalDelta]k; 


psi[x_?NumericQ, y_?NumericQ, z_?NumericQ, t_?NumericQ] := 
  1/(2 Pi Sqrt[Pi] \[CapitalDelta]k)^(3/2)  NIntegrate[
    Exp[-kx^2/(2 \[CapitalDelta]k^2)] Exp[-ky^2/(2 \
\[CapitalDelta]k^2)] Exp[-(kz - k0)^2/(2 \[CapitalDelta]k^2)] Exp[
      I kx x] Exp[I ky y] Exp[
      I kz (z - z0)] Exp[-I k v (t - t0)] {Sqrt[(1 + kz/k)/
        2], (kx + I ky)/Sqrt[2 (k^2 + k kz)]}, {kx, -kInfinity, 
     kInfinity}, {ky, -kInfinity, kInfinity}, {kz, -kInfinity + k0, 
     kInfinity + k0}, 
    Method -> {"LocalAdaptive", Method -> "GaussKronrodRule"}];

Excuting using:

psi[0.01, 0.01, z0/2, t0/2] // Timing
$\endgroup$
  • $\begingroup$ I get five times faster result using the default options of NIntegrate. What speed-up are you looking for? (5 times, 100 times.) In what ranges of the parameters? $\endgroup$ – Anton Antonov Jul 13 '16 at 11:57
  • $\begingroup$ @AntonAntonov If it is possible, less 1s is acceptable, because later on I need to integrate w.r.t x, y, z. Note that the defaut option returns the result is not accurate, because I have an analytical expression. The result returned by my code is close to that expression. $\endgroup$ – an offer can't refuse Jul 13 '16 at 12:04
  • $\begingroup$ Yeah, I noticed the discrepancies in the integral estimates. Thanks. $\endgroup$ – Anton Antonov Jul 13 '16 at 12:07
  • $\begingroup$ I think it is better to separate your integral into two integrals -- you function psi calls NIntegrate with a vector of two elements. By separating into two integrals you can investigate and optimize the computations better. If I use the first term only, Sqrt[(1 + kz/k)/2], NIntegrate finishes in less than a second. With the second term I can get some speed-up using certain options values... $\endgroup$ – Anton Antonov Jul 13 '16 at 12:40
  • $\begingroup$ can you post it? seperate is fine $\endgroup$ – an offer can't refuse Jul 13 '16 at 12:42
2
$\begingroup$

I did manage to get a 6-7 times speed up by separating the definition of psi into two integrals and experimenting with the option values. (I am not sure how significant the speed-up is... OP mentioned in the comments that 1 second is the goal.)

Here are the two new functions:

Clear[psi1]
psi1[x_?NumericQ, y_?NumericQ, z_?NumericQ, t_?NumericQ, opts : OptionsPattern[]] := 
  1/(2 Pi Sqrt[Pi] \[CapitalDelta]k)^(3/2) NIntegrate[
    Exp[-kx^2/(2 \[CapitalDelta]k^2)] Exp[-ky^2/(2 \
\[CapitalDelta]k^2)] Exp[-(kz - k0)^2/(2 \[CapitalDelta]k^2)] Exp[
      I kx x] Exp[I ky y] Exp[
      I kz (z - z0)] Exp[-I k v (t - t0)] (Sqrt[(1 + kz/k)/
        2]), {kx, -kInfinity, kInfinity}, {ky, -kInfinity, 
     kInfinity}, {kz, -kInfinity + k0, kInfinity + k0}, opts];

Clear[psi2]
psi2[x_?NumericQ, y_?NumericQ, z_?NumericQ, t_?NumericQ, opts : OptionsPattern[]] := 
  1/(2 Pi Sqrt[Pi] \[CapitalDelta]k)^(3/2) NIntegrate[
    Exp[-kx^2/(2 \[CapitalDelta]k^2)] Exp[-ky^2/(2 \
\[CapitalDelta]k^2)] Exp[-(kz - k0)^2/(2 \[CapitalDelta]k^2)] Exp[
      I kx x] Exp[I ky y] Exp[
      I kz (z - z0)] Exp[-I k v (t - t0)] ((kx + I ky)/
       Sqrt[2 (k^2 + k kz)]), {kx, -kInfinity, 
     kInfinity}, {ky, -kInfinity, kInfinity}, {kz, -kInfinity + k0, 
     kInfinity + k0}, opts];

The functions are almost identical except for the terms coming from the vector {Sqrt[(1 + kz/k)/2], (kx + I ky)/Sqrt[2 (k^2 + k kz)]} in the original definition of psi.

Some computational results and timings follow. (Note that I have posted results with two different precision goals 5 and 6. With precision goal 3 the computations are even faster, and the obtained results reach the precision goal 6.)

Fastest computations adhering to communication with OP

In[127]:= psi1[1/100, 1/100, z0/2, t0/2, Method -> {"GlobalAdaptive"}, 
  MinRecursion -> 2, PrecisionGoal -> 3] // AbsoluteTiming

Out[127]= {0.070369, 0.0000133329 - 6.65116*10^-7 I}

In[126]:= psi2[1/100, 1/100, z0/2, t0/2, 
  Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
  PrecisionGoal -> 3] // Timing

Out[126]= {0.729304, -5.95563*10^-12 + 7.27904*10^-12 I}

PrecisionGoal->5

In[111]:= 
psi1[1/100, 1/100, z0/2, t0/2, 
  Method -> {"LocalAdaptive", Method -> "GaussKronrodRule"}, 
  PrecisionGoal -> 5] // AbsoluteTiming

Out[111]= {0.787942, 0.0000133343 - 6.65033*10^-7 I}


In[112]:= 
psi2[1/100, 1/100, z0/2, t0/2, 
  Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
  PrecisionGoal -> 5] // Timing

Out[112]= {3.16929, -5.95563*10^-12 + 7.27904*10^-12 I}

PrecisionGoal->6

In[109]:= 
psi1[1/100, 1/100, z0/2, t0/2, 
  Method -> {"LocalAdaptive", 
    Method -> "GaussKronrodRule"}, PrecisionGoal -> 6] // AbsoluteTiming

Out[109]= {1.2058, 0.0000133343 - 6.65033*10^-7 I}


In[110]:= 
psi2[1/100, 1/100, z0/2, t0/2, 
  Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
  PrecisionGoal -> 6] // Timing

Out[110]= {7.79748, -5.95563*10^-12 + 7.27904*10^-12 I}
| improve this answer | |
$\endgroup$
  • $\begingroup$ Is parallel possible for this problem? I have a multi core machine... $\endgroup$ – an offer can't refuse Jul 13 '16 at 12:56
  • $\begingroup$ I do not think parallelization would help much in this case. What kind of precision and accuracy goals are you looking for? $\endgroup$ – Anton Antonov Jul 13 '16 at 12:58
  • $\begingroup$ I must stress tho that cubature is in general more difficult to do than quadrature ("curse of dimensionality"). Here, the OP is lucky that some speed-ups were possible. Now that I think about it: what if one uses an actual multidimensional rule like Genz-Malik ("MultidimensionalRule") instead of taking a Cartesian product of one-dimensional rules? $\endgroup$ – J. M.'s technical difficulties Jul 13 '16 at 12:59
  • $\begingroup$ @AntonAntonov accuracy goal about 3 is fine. $\endgroup$ – an offer can't refuse Jul 13 '16 at 13:00
  • $\begingroup$ @buzhidao, any reason why you're not considering (quasi-)Monte Carlo, if you do not need so much accuracy anyway? $\endgroup$ – J. M.'s technical difficulties Jul 13 '16 at 13:02

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.