NIntegrate Oscillating kernel

Question

I have a finite highly oscillating integral of the type:

NIntegrate[x^2 Exp[I f[x] s] , {x, 0, 1}]

where f is a simple function like Log[1 + x]/x. This integral becomes highly oscillatory for large s and there are numerical errors showing up. Mathematica's standard numerical integration methods are only for kernels Exp[I x s].

How do I deal with this? Do I need to manually implement a method? I think this should be possible because the oscillating function is still trigonometric and f varies more slowly than x.

Answer 1

Here are two ways, a simple change of variables and an assist to the Levin Rule.

Change of variables

Block[{f, s},
 integrand = ComplexExpand[
    x^2 Exp[I f[x] s], TargetFunctions -> {Re, Im}
    ] Dt[x];
 f[x_] := Log[1 + x]/x;
 s = 10^7;
 NIntegrate[integrand /. x -> 1/u /. Dt[u] -> 1, {u, Infinity, 1}]
 ]

(*  -2.48219*10^-9 + 5.17734*10^-7 I  *)

Explicit Levin Rule

For details on the Levin Rule, see this:

Block[{f, s},
 li = NIntegrate`LevinIntegrandReduce[x^2 Exp[I f[x] s], x];
 f[x_] := Log[1 + x]/x;
 s = 10^7;
 NIntegrate[li, {x, 0, 1}, Method -> "LevinRule"]
 ]

(*  -2.48219*10^-9 + 5.17734*10^-7 I  *)

Check consistency:

Block[{f, s},
 li = NIntegrate`LevinIntegrandReduce[x^2 Exp[I f[x] s], x];
 f[x_] := Log[1 + x]/x;
 s = 10^7;
 {NIntegrate[li, {x, 0, 1}, Method -> "LevinRule"],
  NIntegrate[li, {x, 0, 1}, Method -> "LevinRule", 
   WorkingPrecision -> 32],
  NIntegrate[li, {x, 0, 1}, Method -> {"LevinRule", "Points" -> 11}, 
   WorkingPrecision -> 24],
  NIntegrate[li, {x, 0, 1}, Method -> {"LevinRule", "Points" -> 21}, 
   WorkingPrecision -> 24]}
 ]
(*
{-2.48219*10^-9 + 5.17734*10^-7 I,
 -2.4821858856455486537866819639768*10^-9 + 
   5.1773393974054904734746810429446*10^-7 I, 
 -2.48218588564554865378666*10^-9 + 
   5.17733939740549047347468*10^-7 I,
 -2.48218588564554865378660*10^-9 +
   5.17733939740549047347468*10^-7 I}
*)

Answer 2

Asymptotics (integration by parts)

Since the integrand is highly oscillatory for large $s$, we can obtain an asymptotic expression by performing integration by parts.

Observe that

$$ \begin{align} & \int_0^1 x^2 \mathrm{e}^{\mathrm{i} f(x) s} \,\mathrm{d}x \\ &=\int_0^1 \frac{x^2}{\mathrm{i} f'(x) s} \,\mathrm{d}\left( \mathrm{e}^{\mathrm{i} f(x) s} \right) \\ &= \left[ \frac{x^2 \mathrm{e}^{\mathrm{i} f(x) s}}{\mathrm{i} f'(x) s} \right]_0^1 - \int_0^1 \mathrm{e}^{\mathrm{i} f(x) s} \,\mathrm{d} \left( \frac{x^2}{\mathrm{i} f'(x) s} \right). \end{align} $$

The first term (boundary term) is already small, since it is of order $1/s$. Assuming $f$ is sufficiently well behaved, the second term is even smaller, since we have a highly oscillatory integrand of order $1/s$ (I would expect the second term to be of order $1/s^2$).

So we merely need to evaluate the boundary term:

f[x_] := Log[1 + x] / x
boundaryTerm[s_][x_] := x^2 Exp[I f[x] s] / (I f'[x] s);
result = boundaryTerm[s][1] - Limit[boundaryTerm[s][x], x -> 0]

$$-\frac{\mathrm{i} \, 2^{\mathrm{i} s}}{s (1/2 - \log 2)} $$

Table[{s, result // N}, {s, 10^Range[4, 7]}] // TableForm

$s$	result
$10^4$	$-0.000465659 + 0.00022631 \,\mathrm{i}$
$10^5$	$0.0000508565 + 9.70392 \cdot 10^{-6} \,\mathrm{i}$
$10^6$	$4.92323 \cdot 10^{-6} + 1.60227 \cdot 10^{-6} \,\mathrm{i}$
$10^7$	$-2.48291 \cdot 10^{-9} + 5.17734 \cdot 10^{-7} \,\mathrm{i}$

This is consistent with Michael E2's result.

Addendum: more terms

We can obtain a full asymptotic expansion by performing integration by parts iteratively. We get $$ \int_0^1 x^2 \mathrm{e}^{\mathrm{i} f(x) s} \,\mathrm{d}x = \frac{A_1}{s} + \frac{A_2}{s^2} + \frac{A_3}{s^3} + \text{etc.}, $$ where $$ \begin{align} T_1 (x) &= \frac{x^2}{\mathrm{i} f'(x)}, \\ T_n (x) &= \frac{T_{n-1}' (x)}{\mathrm{i} f'(x)}, \end{align} $$ and $$ A_n = (-1)^{n + 1} \left[ T_n (x) \mathrm{e}^{\mathrm{i} f(x) s} \right]_0^1. $$

In the case of $s = 10^7$, taking 6 terms allows us to achieve all of the decimal places in Michael E2's WorkingPrecision -> 32 check:

f[x_] := Log[1 + x] / x

t[1] = x^2 / (I f'[x]);
t[n_] := t[n] = D[t[n - 1], x] / (I f'[x]);

a[n_] := a[n] =
  (-1)^(n + 1) * Subtract[
    Limit[t[n] Exp[I f[x] s], x -> 1],
    Limit[t[n] Exp[I f[x] s], x -> 0]
  ];

asympTerm[n_] := a[n] / s^n;

nMax = 6;
asympSum = Sum[asympTerm[n], {n, nMax}];
asympSum /. {s -> 10^7} // N[#, 34] &
(*
  -2.4821858856455486537866819639768*10^-9
  +5.177339397405490473474681042944648*10^-7 I
*)

(* Check size of omitted term *)
asympTerm[nMax + 1] /. {s -> 10^7} // N
(*
  4.98211*10^-47 - 4.15958*10^-42 I
*)

WARNING. This is an asymptotic expansion, NOT a convergent expansion. At some point the terms will start getting larger, and the series will NOT converge. Always check the size of the first omitted term.

Answer 3

Once again demonstrating the power of Hadamard's maxim, let us use a parabolic contour to evaluate this oscillatory integral:

parabolic[a_, x_] = Simplify[InterpolatingPolynomial[{{0, 0}, {1/2, -a}, {1, 0}}, x]];

With[{s = 1*^7, a = 1/2}, 
     NIntegrate[With[{x = x + I parabolic[a, x]}, 
                     x^2 Exp[I Log[1 + x]/x s]]
                (1 + I Derivative[0, 1][parabolic][a, x]), {x, 0, 1}, 
                MaxRecursion -> 12, WorkingPrecision -> 25]]
   -2.4821858856388198658695*10^-9 + 5.177339397405635636628810*10^-7 I

Answer 4

If you use the method "LocalAdaptive", the integration routine of MMA will put more sample points in the region where the integrand oscillates strongly. Here is an example:

NIntegrate[x^2 Exp[I Log[1 + x]/x 10000000], {x, 0, 1}, 
 Method -> "LocalAdaptive"]

This gives:

9.90376*10^-8 + 0.0000697064 I