# Levin Rule Integrate Power of Oscillating Function

I want to numerically integrate a function of the form:

    NIntegrate[Exp[I a x] x^-n (1+Exp[I b x])^n,{x,1,10^4},
Method -> {"LevinRule", "LevinFunctions" -> {"TrigRelated"}}]


BTW, Mathematica understands this better using Euler's formula exp = cos + i sin .

The Levin Rule wants a explicit oscillating kernel like Exp[I 1/50 x] . The obvious solution would be to expand the bracket using the binomial theorem and then there would be contributions like Exp[I (a+l b) x]  - oscillations with well defined frequency. But then I need to compute N integrals instead of 1 and there will be cancellations. Is there a smarter way of telling the Levin Rule what to do than the Binomial expansion?

• With: g[a_, b_, n_, Z_] := NIntegrate[Exp[I a x] x^-n (1 + Exp[I b x])^n, {x, 1, Z}, WorkingPrecision -> 30]; f[1/500, 1/500, 10, 10^4] standard method works fine. Jan 16 at 12:14

The automatic Levin-rule reduction of the integrand picks $$e^{iax}$$ for the kernel and $$e^{-ibx}/x^n$$ for the amplitude, which fails to capture the whole oscillatory component. You can rescue the Levin rule manually. See this tutorial for details.

We need to construct a basis $$\{w_j\}$$, the Levin "kernel," that satisfies a homogeneous linear ODE $$w_j'(x)=\sum_k A_{j,k}(x) w_k(x) \,.$$ For this problem, we may pick for the kernel the following: $$w_j(x) = e^{iax}\left(1+e^{ibx}\right)^j\,, \quad j=0,\dots,n \,.$$ A simple calculation (luckily) shows $$w_j'(x)=i a \, w_j(x)+i b j \left(w_j(x)-w_{j-1}(x)\right) \,,$$ which is programmatically valid for $$j=0$$ since $$w_{j-1}$$ is multiplied by $$j$$.

In the code below, only the kernel levinK[n] and the differential matrix $$A$$ levinDM[n] are needed to implement the Levin rule in NIntegrate. The others were written for checking the ODE above. While levinDM[] uses levinODE[] and levinVARS[], it could directly compute the matrix.

ClearAll[levinODE, levinICS];
levinODE[n_] :=
Table[w[j]'[x] == I*a*w[j][x] + I*b*j*(w[j][x] - w[j - 1][x]), {j, 0, n}];
levinICS[n_] := Table[w[j] == 2^j, {j, 0, n}];
levinVARS[n_] := Array[w, n + 1, 0];
levinK[n_] := Table[Exp[I a x] (1 + Exp[I b x])^j, {j, 0, n}];
levinDM[n_] := CoefficientArrays[
levinODE[n][[All, -1]],
Through[levinVARS[n][x]]
][]; (* Part 2 are the coefficients of the linear part *)


We have to specify the "Amplitude", or NIntegrate will choose $$w_0(x)$$ for determining it instead of $$w_n(x)$$.

Block[{n = 4, a = 3, b = 4},
NIntegrate[Exp[I a x] x^-n (1 + Exp[I b x])^n, {x, 1, 10^4},
Method -> {"LevinRule"
, "Amplitude" -> SparseArray[{{n + 1} -> x^-n}, {n + 1}]
, "Kernel" -> levinK[n]
, "DifferentialMatrix" -> levinDM[n]
, "Points" -> 20 + 10 n (* for a=3, b=4, n up to 8 or higher *)
}]
]

(*  -0.01205 + 0.0929632 I  *)


The number of collocation "Points" needed depends on how oscillatory the integrand is (a and b) and the dimension of the Levin system (n). I don't know off-hand a formula for it. For a and b on the order of $$10^{-3}$$, the formula 5 + 5 n works for small n.

I thought it was cool how easy it was to do the Levin rule by hand, but there's an easier way. @J.M.'s advice for the OP's previous oscillatory integral applies here, too. A complex integration path can dampen the oscillations, something like {x, 1, 1 + 2I, 10^4 + 2I, 10^4}. It's considerably faster than the Levin rule. Since the integrand is analytic in a simply connected, complex neighborhood of 1 <= x <= 10^4 (for example, Re[x] >= 0), a high-order Gauss rule will be beneficial.

Block[{n = 4, a = 3, b = 4},
NIntegrate[
Exp[I a x] x^-n (1 + Exp[I b x])^n,
{x, 1, 1 + 2 I, 10^4 + 2 I, 10^4},
Method -> {"GaussKronrodRule", "Points" -> 21}]
]

(*  -0.01205 + 0.0929632 I  *)