# NIntegrate Oscillating kernel

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.

• Dec 30, 2020 at 14:03
• What are the magnitudes of the "large s"? E.g $s \in [10,100]$ or $s \in [10^5,10^6]$ ? Dec 30, 2020 at 14:06
• @AntonAntonov maybe up to 10^{4,5,6,7} but I dunno for sure depends on what I'll do later with this integral.
– VN23
Dec 30, 2020 at 14:10

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 = NIntegrateLevinIntegrandReduce[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 = NIntegrateLevinIntegrandReduce[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}
*)


# 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] - 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.

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 = 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.

• Cool. This is a better problem than I first thought when I saw the OP. Dec 31, 2020 at 6:13
• Very nice. It's a pity AsymptoticIntegrate[] does not work on the OP's problem... Dec 31, 2020 at 9:37

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

• Nice. One can use a polygonal path as well: NIntegrate[x^2 Exp[I Log[1 + x]/x 10^7], {x, 0, 1/2 - I/2, 1}, WorkingPrecision -> 25, MaxRecursion -> 12] Dec 31, 2020 at 2:28
• Yes, that works nicely too, @Michael. A serious effort (and a more complicated integrand) would probably have to delve into determining "steepest descent contours", but the simple approaches you and I give are already good enough for this situation. Dec 31, 2020 at 3:05

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

9.90376*10^-8 + 0.0000697064 I

• I wouldn't trust these results.Even NIntegrate[x^2 Exp[I Log[1 + x]/x 10000000], {x, 0, 1}, Method -> "DoubleExponentialOscillatory", WorkingPrecision -> 30, MaxRecursion -> 10000] error is about 0.035 !. Dec 30, 2020 at 13:57