# Compute integral symbolically or numerically

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}])];];

-
I haven't tried your code, but it seems like you meant to evaluate withTransforms[int] before passing to (N)Integrate... otherwise, it doesn't make sense to call that on the output of NIntegrate, which is just a number. It could possibly change how your integral behaves, but certainly will change your answer. –  rm -rf Feb 26 '12 at 19:03
@RM Thanks! At least for some integrands, that helps (I don't get that error anymore) but still Mathematica takes a very long time (at this moment it didn't finish yet) to compute that integral. If I use your modification on IntoN it gives for the same integrand the same result for the unmodified IntoS. I'm not sure what you mean with "change your answer". –  Jonas Teuwen Feb 26 '12 at 19:15
maybe I am just missing it, but I don't see what k is supposed to be (I changed NIntegrate to nintegrate to see what would the integrand would be, and a k is left in there). I also tried setting k=1 and plotting the real and imaginary parts of the integrand and see no obvious problems; but perhaps I am looking at the wrong integrand. –  acl Feb 26 '12 at 20:39
@acl: $k$ is $2 \pi$. I have added this now. Does the integration work for you? –  Jonas Teuwen Feb 26 '12 at 20:43
OK, I see. No, I have the same problem that you report, and I don't see why. I tried Plot3D[Im[withTransforms[intgnd]],{a, 0, Pi/4},{b, 0, 2*Pi}] (here intgnd is the first expression in your question), and see no problem (also, no problem is evident with its real part). –  acl Feb 26 '12 at 20:54

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.

-
Thanks, this works for certain integrands, but not for all unfortunately. (Like the last one I gave). –  Jonas Teuwen Feb 27 '12 at 21:35

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.

-
Okay, I was already thinking that it must be some kind of round-off problem. I will attempt to scale the coefficients. Thanks! –  Jonas Teuwen Feb 26 '12 at 21:14
This might be interesting for others, no? In that case I would leave this as an answer. Something else: If I scale the coefficients, some other coefficients become quite small and then those integrals fail to compute. Argh! –  Jonas Teuwen Feb 26 '12 at 21:43
@JonasTeuwen I will leave it then. I have also added a description of what happens if I split the domain of integration in two (I don't get any warning then, but the two parts are nearly equal in magnitude). It is strange that you get such difference in magnitudes of coeffs; is there some physical reason for this? Perhaps this can be prevented earlier in the calculation? –  acl Feb 26 '12 at 22:14
Well, I can't really find a physical reason for this, but the calculations sometimes give $\textrm{Large coefficient}$ and then again $\textrm{Large coefficient}^{-1}$. I have edited the start post again, maybe these integrals can be computed (I have already scaled the coefficients). –  Jonas Teuwen Feb 26 '12 at 22:29
I will try to look at this tomorrow if I can find some time. About automatically deleting small terms, have you looked at Chop? –  acl Feb 26 '12 at 23:48