Compute integral symbolically or numerically

Question

I want to compute the integral of the following integrand

((0.000245587 + 0.0000651908 I) E^(-I kx)
   kx ky (kz^2 - 4 \[Pi]^2))/((kx^2 + ky^2) kz) - ((4.00033*10^-17 - 
    3.51178*10^-18 I) E^(-I kx) (kx^3 kz + kx ky^2 kz + kx^2 kz^2 + 
    4 ky^2 \[Pi]^2))/((kx^2 + ky^2) kz)

Where $k = 2 \pi$.

For this I have the following code:

SetAttributes[withTransforms, HoldAll];
withTransforms[code_] := 
 Block[{kx = k*Sin[a]*Cos[b], 
   ky = k*Sin[a]*Sin[b], kz = k*Cos[a]}, 
  code]

And the symbolic integration:

IntoS[int_] := 
 withTransforms[
  Integrate[
   int*k^2*Sin[a]*Cos[a], {a, 0, Pi/4}, {b, 0, 2*Pi}]]

And finally, the numerical:

IntoN[int_] := 
 withTransforms[
  NIntegrate[
   int*k^2*Sin[a]*Cos[a], {a, 0, Pi/4}, {b, 0, 2*Pi}]]

Now, IntoS, the symbolic integration doesn't give a result within 5 minutes (after which I have aborted the evaluation) on my core 2 duo. I need to compute many of these integrals for a plot and if it takes several minutes already to evaluate one point that won't be very useful.

Furthermore, IntoN throws

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

at me. How do I integrate this function?

Edit Scaling coefficients does not really work. I always get that some coefficients are quite large and others very small. This way, the integration never works. I can delete the "small" terms manually but that it not really automatic.

Edit 2: Mathematica is unable to compute the following integral numerically:

(E^(I kx) kx ky (-1 + kz^2))/(Sqrt[2] (kx^2 + ky^2) kz)

Symbolically this integral is $0$. However, I do need the numerical computation as the symbolic computation can take quite long and I need to plot the end-result. This is a special case of the term

(E^(I (kx x + ky y + kz z)) kx ky (-1 + kz^2))/(Sqrt[2] (kx^2 + 
   ky^2) kz)

And I want to plot in $x$ (with $y, z = 0$ for example) and this integral is sometimes impossible to compute symbolically.

Edit 3: I have "fixed" it by writing my own numerical method:

Trapezoidal2D[f_, {x_, a0_, b0_}, {y_, c0_, d0_}, m0_] :=

  Module[{h = (b0 - a0)/m0, k = (d0 - c0)/m0, a = N[a0], b = N[b0], 
    c = N[c0], d = N[d0], i, j, m = m0, X, Y},
   Subscript[X, k_] = a + i h;
   Subscript[Y, j_] = c + j k;
   Return[
    1/4 h k (Function[{x, y}, f][a, c] + Function[{x, y}, f][b, d] + 
       Function[{x, y}, f][a, d] + Function[{x, y}, f][b, d] + 
       2 Sum[Function[{x, y}, f][Subscript[X, i], c], {i, 1, 
          m - 1}] + 
       2 Sum[Function[{x, y}, f][Subscript[X, i], d], {i, 1, 
          m - 1}] + 
       2 Sum[Function[{x, y}, f][a, Subscript[Y, j]], {j, 1, 
          m - 1}] + 
       2 Sum[Function[{x, y}, f][b, Subscript[Y, j]], {j, 1, 
          m - 1}] + 
       4 Sum[Function[{x, y}, f][Subscript[X, i], Subscript[Y, 
          j]], {i, 1, m - 1}, {j, 1, m - 1}])];];

Answer 1

For the symbolic attempt I'd first rationalize to get exact input.

e1 = ((0.000245587 + 0.0000651908 I) E^(-I kx) kx ky (kz^2 - 
        4 \[Pi]^2))/((kx^2 + ky^2) kz) - ((4.00033*10^-17 - 
        3.51178*10^-18 I) E^(-I kx) (kx^3 kz + kx ky^2 kz + 
        kx^2 kz^2 + 4 ky^2 \[Pi]^2))/((kx^2 + ky^2) kz);
e2 = Rationalize[e1, 0]

Out[5]= ((245587/1000000000 + (162977 I)/2500000000) E^(-I kx)
   kx ky (kz^2 - 4 \[Pi]^2))/((kx^2 + ky^2) kz) - (1/((kx^2 + 
    ky^2) kz))(1/24997937670142213 - I/
    284755878785117533) E^(-I kx) (kx^3 kz + kx ky^2 kz + kx^2 kz^2 + 
    4 ky^2 \[Pi]^2)

I believe this does the transformation you do in a different way.

e3 = 
 e2*k^2*Sin[a]*Cos[a] /. {kx -> k*Sin[a]*Cos[b], 
    ky -> k*Sin[a]*Sin[b], kz -> k*Cos[a]} /. k -> 2*Pi

With this, we can now do as below.

In[26]:= Timing[ii = Integrate[e3, {a, 0, Pi/4}, {b, 0, 2*Pi}]]

Out[26]= {158.87, 
 Integrate[(8/284755878785117533 + (8*I)/24997937670142213)*Pi^3*
       (-2*I*Pi*BesselJ[2, 2*Pi*Sin[a]]*Cos[a]^2*Sin[a] + 

     BesselJ[1, 2*Pi*Sin[a]]*(I*(1 + Cos[a]^2) + Pi*Sin[a]*Sin[2*a])), 
     {a, 0, Pi/4}]}

In[28]:= N[ii]

Out[28]= 7.85307*10^-16 + 1.11299*10^-15 I

Whether this is accurate will depend on cancellation error and other vagaries of the quadrature.

Answer 2

If I define

expr = Function[{a, b},
   Evaluate@withTransforms[intgnd]];

(with intgnd the expression given at the top of the question), I get tiny differences between points symmetric around $b=\pi$:

ListLogPlot[Abs[expr[.3, #] + expr[.3, 2 Pi - #] & /@ Range[0, 2 Pi, .01]]]

so even if your coefficients are accurate, you will have lots of trouble obtaining good results unless you greatly increase working precision (and force the precision of your input numbers to be higher than machine precision).

If you can, it's probably better to rescale your variables at some earlier point to avoid such extreme numbers.

One may also see what is going on as follows:

t={
 NIntegrate[expr[a, b], {a, 0, Pi/4}, {b, 0, Pi}], 
 NIntegrate[expr[a, b], {a, 0, Pi/4}, {b, Pi, 2 Pi}]
 }

(*{-0.0000297965 + 0.000112249 I, 0.0000297965- 0.000112249 I}*)

Plus@@t
(*-2.23895*10^-16 + 1.1891*10^-16 I*)

ie, the two halves are very close to being equal in magnitude.