# Numerical integration of oscillatory function

I am trying to evaluate the following oscillatory integration numerically but the answers are wrong or not accurate enough. I don't know what the problem is.

B = 2.5*10^(-6); l = 0.02; V = 0.01;
X[t_] := -((4 l t)/(B V)) - (4 I l (Exp[-V t] - 1))/(V (Exp[I B V] - 1))
k[w_] := Re[NIntegrate[Exp[I w t] Exp[X[t]], {t, 0, \[Infinity]}]]


I know from theoretical calculation that $k(w)/k(-w) =\exp(B w)$ but as you see below it is not true for the Mathematica answers for large $w$.

{k[2*10^6]/k[-2*10^6], Exp[B 2*10^6]}

{2.05121, 148.413}


I have tried changing the Accuracygoal and WorkingPrecision but it didn't help.

• With more precision, I can't make Mathematica think the answer is anything other than 9.571…. (Changing all the inexact numbers to exact numbers means Mathematica is free to use greater precision than machine precision.) – Patrick Stevens Sep 27 '15 at 10:39
• And with $w=3$ I obtain $k(w)/k(-w)$ being very close to 1, not $e^5$. – Patrick Stevens Sep 27 '15 at 10:41
• Are you sure that $\exp(Bw)$ is the theoretical convergence for large $w$ rather than small $w$? – Taiki Sep 28 '15 at 12:32
• For higher precision and accuracy goals computations it is better to use definitions with exact numbers. I.e. B = 25/10 * 10^6; l = 2/100; V=1/100 . – Anton Antonov Sep 28 '15 at 13:50

The integral can be solved exactly:

k[w_] = Integrate[Exp[I w t] Exp[X[t]], {t, 0, ∞}]
(*
ConditionalExpression[(
4^((-4 l + I B V w)/(B V^2)) * E^((4 I l)/((-1 + E^(I B V)) V)) *
((I l)/((-1 + E^(I B V)) V))^((-4 l + I B V w)/(B V^2)) *
(Gamma[(4 l - I B V w)/(B V^2)] -
Gamma[(4 l - I B V w)/(B V^2), (4 I l)/((-1 + E^(I B V)) V)])
)/V,
Re[V] > 0]
*)


The plot seems to agree with bbgodfrey's:

Block[{B = 25*10^(-7), l = 2/100, V = 1/100},
Table[N[{w, k[w]/k[-w]}, 3],
{w, 1, 10000000, 100000}]
] // ListLinePlot


As suggested by Patrick Stevens, it is highly desirable to use exact numbers.

B = 25*10^(-7); l = 2/100; V = 1/100;


Further, X[t]is entirely negligible for t > 1/10.

LogPlot[Evaluate[ReIm[Exp[X[t]]]], {t, 0, 1/10}, PlotRange -> All]


allowing the upper limit on the integral to be set to 1/10, which reduces running time.

k[w_] := Re[NIntegrate[Exp[I w t] Exp[X[t]], {t, 0, 1/10}, WorkingPrecision -> 30]]


High WorkingPrecision is essential to obtain reasonable accuracy at large w. The ratio, k[w]/k[-w], then can be computed at reasonable speed over a wide range of w.

Quiet@Plot[k[w]/k[-w], {w, 0, 10000000}]


In response to a question by Patrick Stevens abount truncating the range of integration, I point out that Re[X[t]], the exponential of which is the magnitude of the integrand, is given by

FullSimplify[Re[X[t]] // ComplexExpand, t > 0]
(* E^(4 (-800000 t + E^(-t/100) (-1 + E^(t/100)) Cot[1/80000000])) *)


which is strictly negative and decreases linearly for large t.

Plot[Evaluate[ReIm[X[t]]], {t, 0, 200}, PlotRange -> All]


Thus, the integrand decreases exponentially for large t. That the truncation also makes sense at small t can be shown by truncating the integral at 1 instead of 1/10. The computation, although somewhat slower, gives results that are indistinguishable to the eye.

• I'm afraid of "is entirely negligible" when you're truncating infinite regions - $1/n$ is entirely negligible for $n > 10^{100}$, and yet the integral of $1/x$ is still divergent… – Patrick Stevens Sep 30 '15 at 17:53
• @PatrickStevens Perhaps, I should point out in the answer that it is decreasing exponentially at large t and already is down to 10^-300 at t = 0.2.. – bbgodfrey Sep 30 '15 at 18:05
• The exponential decrease is enough for me, but I'd prefer that to be proved or at least indicated :P – Patrick Stevens Sep 30 '15 at 18:09
• @PatrickStevens Done – bbgodfrey Sep 30 '15 at 19:01