# Symbolic vs numerical integration

This is related to a question that I asked before but didn't resolve quite well.

I want to obtain the solution of the following expression analytically:

Assuming[a>0 && s<0, Integrate[(w E^(-w/a)Sin[((w - s)^2)/2])/((w - s)^2)/2, {w, 0, ∞}]]


Mathematica gives an analytical solution as:

(1/(4 a))E^(-(s/
a)) (2 a E^((-I + 1/a) s) + 2 a E^((I + 1/a) s) - 4 a E^(s/a) +
2 (a - s + I a s) ExpIntegralEi[(-I + 1/a) s] +
2 a ExpIntegralEi[(I + 1/a) s] - 2 s ExpIntegralEi[(I + 1/a) s] -
2 I a s ExpIntegralEi[(I + 1/a) s] - 4 a ExpIntegralEi[s/a] +
4 s ExpIntegralEi[s/a] - 2 a Log[-1 - I a] + 2 s Log[-1 - I a] +
2 I a s Log[-1 - I a] + a Log[1 - I a] - s Log[1 - I a] +
I a s Log[1 - I a] - 2 a Log[-1 + I a] + 2 s Log[-1 + I a] -
2 I a s Log[-1 + I a] + a Log[1 + I a] - s Log[1 + I a] -
I a s Log[1 + I a] - a Log[-(I/(-I + a))] + s Log[-(I/(-I + a))] +
I a s Log[-(I/(-I + a))] + (a - s - I a s) Floor[Arg[1 + I a]/(
2 \[Pi])] (Log[I + 1/a] + Log[-((I a)/(-I + a))]) -
a Log[I/(I + a)] + s Log[I/(I + a)] - I a s Log[I/(I + a)])


For an appropriate choice of parameters:

a = 10^-6;
s = -10^(-1) Sqrt[10^(-3)^2 + b^2];


The plot of the solution in terms of the parameter b would be as which shows values around 10^250.

The function Sin[((w - s)^2)/2])/((w - s)^2)/2 can be approximated by the Dirac delta function 𝛿(w-s). It seems that Mathematica uses the sifting property of the delta function (Integrate[w E^(-w/a) 𝛿(w-s), {w, 0, ∞}] = -s E^(-s/a)) to approximate the result. If so, since s is negative, the result would be so large.

But if we approximate the Sin function using the first terms of its Taylor expansion, the function converges so quickly to around 10^-12 which seems reasonable. Numerically, this approximation would suffice but I need exact analytical solutions which as I explained are problematic. What is the problem?

• @MarcoB Thanks for the edit. Feb 5 '16 at 19:03
• Since s<0 the quantity (w-s) is strictly positive on [0,Inf] , so I'm not following your delta approximation result. In any case I do not believe Integrate does any such approximation. Feb 5 '16 at 19:17
• @george2079 I am not sure it needs to be negative to define delta function see here. One problem is that the sifting property needs to be defined over a symmetric domain (-Inf,Inf), which for my case clearly is not (0,Inf). Feb 5 '16 at 19:39
• Thats my point. delta is zero over your entire domain. Your analytic result is going to be difficult if not impossible to evaluate numerically with any accuracy due to working precision issues. Feb 5 '16 at 19:51

Evaluating the integral analytically seems to take far too long on my machine, so I aborted it. Since you mention that you studied the integrated expression, perhaps you could consider showing result you got from the integration in your question.

Having said that, numerical integration suggests that the values of the integral are very small over a wide range of $b$; the values are, however, highly oscillatory:

Clear[integral]
integral[b_?NumericQ] :=
NIntegrate[(w E^(-w/a) Sin[((w - s)^2)/2])/((w - s)^2)/2, {w, 0, Infinity}]

Plot[integral[b], {b, 0, 120}, PlotPoints -> 300] Since numerical results seem to be accessible relatively easily, and your symbolic results are dodgy, you may want to reassess your need for an explicit symbolic integral in this case.

• Thanks for the suggestion. It is really helpful. But you didn't answer my question. I included the analytical solution. If you plot its real part, you will see a very different result. Feb 5 '16 at 19:29
• By the way, how can I insert the plot here? I uploaded the plot in a hosting site and I followed the instructions, but it didn't showed up. Feb 5 '16 at 19:32
• @Farhad use the Insert Graphics button or use the Mathematica upload tool Feb 5 '16 at 20:19
• @Sjoerd C. de Vrie Thanks for the point. I included the plot. Feb 5 '16 at 20:31
• @MarcoB I put the exact code you presented in the above in Mathematica 7, but It says "NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections". What is the problem? Feb 6 '16 at 17:30

Perhaps useful I think this substitution provides a preferable form for the numerical integration:

 exp = ((w E^(-w/a) Sin[((w - s)^2)/2])/((w - s)^2)/2 /.
w -> -a Log[g]) D[-a Log[g], g] f[b_?NumericQ] :=
Block[{a = 10^-6, s = -10^(-1) Sqrt[10^(-3)^2 + b^2]},
- a^2 NIntegrate[exp/a^2, {g, 0, 1}, WorkingPrecision -> 30,
MaxRecursion -> 200]]
ListPlot[Table[{b, f[b]}, {b, 0, 40, 1/2}], Joined -> True] Essentially the same result as MarcoB, but I think it should give confidence in the accuracy.

Note I use ListPlot[Table[...]] instead of Plot so that f is always given an exact rational argument and so can work at high precision.

• Thanks. Can you explain what did you do exactly? Feb 6 '16 at 17:35