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.

  • 1
    $\begingroup$ 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.) $\endgroup$ – Patrick Stevens Sep 27 '15 at 10:39
  • $\begingroup$ And with $w=3$ I obtain $k(w)/k(-w)$ being very close to 1, not $e^5$. $\endgroup$ – Patrick Stevens Sep 27 '15 at 10:41
  • 1
    $\begingroup$ Are you sure that $\exp(Bw)$ is the theoretical convergence for large $w$ rather than small $w$? $\endgroup$ – Taiki Sep 28 '15 at 12:32
  • $\begingroup$ 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 . $\endgroup$ – 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, ∞}]
   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)])
   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

Mathematica graphics


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]

enter image description here

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

enter image description here


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]

enter image description here

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.

  • 1
    $\begingroup$ 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… $\endgroup$ – Patrick Stevens Sep 30 '15 at 17:53
  • $\begingroup$ @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.. $\endgroup$ – bbgodfrey Sep 30 '15 at 18:05
  • $\begingroup$ The exponential decrease is enough for me, but I'd prefer that to be proved or at least indicated :P $\endgroup$ – Patrick Stevens Sep 30 '15 at 18:09
  • 1
    $\begingroup$ @PatrickStevens Done $\endgroup$ – bbgodfrey Sep 30 '15 at 19:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.