5
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ Try this:reference.wolfram.com/language/tutorial/… $\endgroup$ – Mariusz Iwaniuk Dec 30 '20 at 14:03
  • 1
    $\begingroup$ What are the magnitudes of the "large s"? E.g $s \in [10,100]$ or $s \in [10^5,10^6]$ ? $\endgroup$ – Anton Antonov Dec 30 '20 at 14:06
  • $\begingroup$ @AntonAntonov maybe up to 10^{4,5,6,7} but I dunno for sure depends on what I'll do later with this integral. $\endgroup$ – VN23 Dec 30 '20 at 14:10
6
$\begingroup$

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}
*)
$\endgroup$
8
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Cool. This is a better problem than I first thought when I saw the OP. $\endgroup$ – Michael E2 Dec 31 '20 at 6:13
  • $\begingroup$ Very nice. It's a pity AsymptoticIntegrate[] does not work on the OP's problem... $\endgroup$ – J. M.'s ennui Dec 31 '20 at 9:37
4
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ 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] $\endgroup$ – Michael E2 Dec 31 '20 at 2:28
  • $\begingroup$ 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. $\endgroup$ – J. M.'s ennui Dec 31 '20 at 3:05
0
$\begingroup$

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
$\endgroup$
1
  • 3
    $\begingroup$ 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 !. $\endgroup$ – Mariusz Iwaniuk Dec 30 '20 at 13:57

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.