3
$\begingroup$

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?

$\endgroup$
1
  • 2
    $\begingroup$ 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. $\endgroup$ Jan 16 at 12:14
2
$\begingroup$

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][0] == 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]]
     ][[2]]; (* 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"
    , "AdditiveTerm" -> 0
    , "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.

$\endgroup$
2
$\begingroup$

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

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.