# 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
Commented Feb 26, 2012 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". Commented Feb 26, 2012 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
Commented Feb 26, 2012 at 20:39
• @acl: $k$ is $2 \pi$. I have added this now. Does the integration work for you? Commented Feb 26, 2012 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
Commented Feb 26, 2012 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). Commented Feb 27, 2012 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! Commented Feb 26, 2012 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! Commented Feb 26, 2012 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
Commented Feb 26, 2012 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). Commented Feb 26, 2012 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
Commented Feb 26, 2012 at 23:48