# How to increase the evaluation speed of this numerical integral?

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

• 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? – Anton Antonov Jul 13 '16 at 11:57
• @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. – an offer can't refuse Jul 13 '16 at 12:04
• Yeah, I noticed the discrepancies in the integral estimates. Thanks. – Anton Antonov Jul 13 '16 at 12:07
• 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... – Anton Antonov Jul 13 '16 at 12:40
• can you post it? seperate is fine – an offer can't refuse Jul 13 '16 at 12:42

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:= psi1[1/100, 1/100, z0/2, t0/2, Method -> {"GlobalAdaptive"},
MinRecursion -> 2, PrecisionGoal -> 3] // AbsoluteTiming

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

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

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


### PrecisionGoal->5

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

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

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

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


### PrecisionGoal->6

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

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

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

Out= {7.79748, -5.95563*10^-12 + 7.27904*10^-12 I}

• Is parallel possible for this problem? I have a multi core machine... – an offer can't refuse Jul 13 '16 at 12:56
• I do not think parallelization would help much in this case. What kind of precision and accuracy goals are you looking for? – Anton Antonov Jul 13 '16 at 12:58
• 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? – J. M.'s technical difficulties Jul 13 '16 at 12:59
• @AntonAntonov accuracy goal about 3 is fine. – an offer can't refuse Jul 13 '16 at 13:00
• @buzhidao, any reason why you're not considering (quasi-)Monte Carlo, if you do not need so much accuracy anyway? – J. M.'s technical difficulties Jul 13 '16 at 13:02