NIntegrate: how to obtain a more precise plot

Question

I want to plot a function involving the integral

a = 20;
d = 10^(-1);

int[t_]:= NIntegrate[E^(-w/a) Sin[(w-d)t]/(w-d)^2,{w,0,d,∞},Method->"PrincipalValue"]

Plot[int[t],{t,0,1000]

Ideally the plot should show damped perfect oscillations. But because of the numerical method of integration, the oscillations are not symmetric:

How can I obtain a more precise plot having perfect osccilations?

Answer 1

here is how I did that integral analytically, with a bit of manual manipulation to take care of the principal value:

a = 20;
d = 10^(-1);
With[{t = 100}, 
     NIntegrate[E^(-w/a) Sin[(w - d) t]/(w - d)^2, {w, 0, d, \[Infinity]},
       Method -> "PrincipalValue"]]

-1.05077

There are a bunch of steps to get here, in a nutshell break the integral into three parts {0,d} , {d, 2d} , {2 d,Inf} , do a variable substitution to make the second interval also {0,d} and add the first two integrands together, which eliminates the singularity.

    intinf =
     Integrate[E^(-w/a) Sin[(w - d) t]/(w - d)^2, {w, 2 d, \[Infinity]}, 
        Assumptions -> {t > 0, d > 0, a > 0}] +
     Integrate[
        Sin[w t]/w^2 (E^(-(w + d)/a) - E^((w - d)/a) ), {w, 0, d}, 
        Assumptions -> {t > 0, d > 0, a > 0}] // FullSimplify 
    % /. t -> 100 // N

(-2 (I + 20 t) ExpIntegralEi[1/200 - (I t)/10] + 800 E^(1/200) Sin[t/10] - 2 (-I + 20 t) (CoshIntegral[1/200 + (I t)/10] + SinhIntegral[1/200 + (I t)/10]))/(80 E^(1/200))

-1.05077 + 0. I

note I expect this should work as well with a and d symbolic but it will take much longer to evaluate the integral.

Plot[Chop@int, {t, 0, 1000}, PlotRange -> {-3, 6}]

Answer 2

I would turn your integral into an ODE. Here is an inactive form of your integral:

a = 20;
d = 1/10;
integral = Inactive[Integrate][
    Exp[-w/a] Sin[(w-d)t]/(w-d)^2,
    {w, 0, ∞},
    GenerateConditions->False
]

Inactive[Integrate][(E^(-w/20) Sin[t (-(1/10) + w)])/(-(1/10) + w)^2, {w, 0, ∞}, GenerateConditions -> False]

Let's take 2 derivatives with respect to t, and activate:

Activate @ D[integral, {t, 2}]

(20 (-20 t Cos[t/10] + Sin[t/10]))/(1 + 400 t^2)

That's enough information to create the differential equation, now we need to come up with initial conditions:

Activate[integral /. t->0]
Activate[D[integral, t] /. t->0]

0

-((I π + ExpIntegralEi[1/200])/E^(1/200))

We see that the first derivative is complex. However, this is because integral is divergent, and needs regularization to be evaluated. So, it is actually the real part that is needed. Now, we can solve the ODE:

if = NDSolveValue[
    {
    f''[t] == Activate @ D[integral, {t, 2}],
    f[0] == 0,
    f'[0] == Re @ Activate[D[integral, t] /. t->0]
    },
    f,
    {t, 0, 1000}
];

Plotting this solution is fast:

Plot[if[t], {t, 0, 1000}, PlotRange->All]

Let's compare it to plotting your int function:

p1 = Plot[int[t], {t, 0, 100}, PlotRange->All]; //AbsoluteTiming
p2 = Plot[if[t], {t, 0, 100}, PlotRange->All]; //AbsoluteTiming

p1
p2

{19.7421, Null}

{0.054623, Null}

Finally, as @george2079 says, it is possible to find the integral analytically. One method is to use DSolveValue:

DSolveValue[
    {
    f''[t] == Activate@D[integral, {t, 2}],
    f[0] == 0,
    f'[0] == Re @ Activate[D[integral, t] /. t->0]
    },
    f,
    t
]

Function[{t}, (1/( 200 E^(1/200)))(-200 E^(1/200) t + 200 t Cosh[1/200] + 200 t CosIntegral[I/200] + 5 I CosIntegral[I/200 - t/10] - 100 t CosIntegral[I/200 - t/10] - 5 I CosIntegral[I/200 + t/10] - 100 t CosIntegral[I/200 + t/10] - 200 t ExpIntegralEi[1/200] + 2000 E^(1/200) Sin[t/10] + t Sinc[I/200] + 200 t SinhIntegral[1/200] + 5 SinIntegral[I/200 - t/10] + 100 I t SinIntegral[I/200 - t/10] - 5 SinIntegral[I/200 + t/10] + 100 I t SinIntegral[I/200 + t/10])]