# Speeding up trigonometric integral

Context

On a possible non trivial toric topology for the Universe (nothing less!).

Problem

I would like to carry out the following integral for $\ell=2,4\cdots 20$.

$$\int _0^{\pi }\int _0^{2 \pi } \sin(\theta ) P_{\ell }\left(\frac{1}{2} \sin(\theta ) (\cos(\phi )-\sin(\phi ))\right) P_{\ell }\left(\frac{1}{2} \sin(\theta ) (\cos(\phi )+\sin (\phi ))\right) \times$$ $$P_{\ell }\left(\frac{\sin (\theta ) ((K+1) \cos (\phi )+(1-K) \sin (\phi ))}{2 \sqrt{K^2+2 \mu K+1}}\right)d\phi d\theta$$ where $P_\ell$ are Legendre Polynomials.

For $\ell$ larger than say 8, Mathematica takes forever and runs out of memory. I have found out a method to circumvent the memory problem (see attempt below) but it still takes very long time to carry out the integration for $\ell>10$.

Question

Is there a way to be smart about this class of integral? Another approach to the one below?

Attempt

I have defined the integrant as

integ[ℓ_] :=
LegendreP[ℓ, (Cos[p] - Sin[p]) Sin[t]/
2] LegendreP[ℓ, (Cos[p] + Sin[p]) Sin[t]/
2] LegendreP[ℓ, ((1 + K) Cos[p] + (1 - K) Sin[p]) Sin[
t]/2/Sqrt[1 + 2 K μ + K^2]] Sin[t]


And the following integration rule:

r1 = {Exp[ Complex[0, b_] t +  p Complex[0, c_] ] ->
Integrate[Exp[ I b  t + I c p], {t, 0, Pi}, {a, 0, 2 Pi}],
Exp[ Complex[0, b_] t ] ->
Integrate[Exp[ I b  t], {t, 0, Pi}, {a, 0, 2 Pi}],
Exp[   p Complex[0, c_] ] ->
Integrate[Exp[ I c p], {t, 0, Pi}, {a, 0, 2 Pi}]}


so that I expand my integrant into sines and cosines into complex exponentials, and carry out the integration via substitution:

   integ // TrigToExp // Expand //
Collect[#, Exp[_]] & // (# /. r1) & // Apart // Simplify

(*-((Pi*(43*K^2 + 84*K*μ + 43))/
(560*(K^2 + 2*K*μ + 1))) *)


A simple timing

Table[{i, integ[i] // TrigToExp // Expand //
Collect[#, Exp[_]] & // (# /. r1) & // Apart //
Simplify;//Timing}, {i, 2, 6,2}]


suggest a $n^{4.5}$ scaling.

• suggesting something to myself: may be use recurrence of LegendreP? Nov 15, 2012 at 17:58
• Is it possible at all to put in values for $\mu$ and $K$ and use NIntegrate? Nov 15, 2012 at 19:58
• @tkott a symbolic solution was needed. Thanks for your interest. Nov 16, 2012 at 7:31

Use the following representation of the Legendre polynomials: $$P_n(x) = 2^n \sum_{k=0}^n x^k \binom{n}{k} \binom{\frac{n+k-1}{n}}{n}$$ Note that the sum effectively is over $k \equiv n \bmod 2$.

Expand each Legendre polynomial into a sum. Integration with respect to $\theta$ is easy: $$\int_0^{\pi} \sin^{k_1+k_2+k_3+1} \theta \mathrm{d}\theta = \operatorname{Beta}\left(\frac{k_1+k_2+k_3+2}{2}, \frac{1}{2}\right) \tag{1}$$ Integration with respect to $\phi$ is more involved. We need the following result, for $p,q \in \mathbb{Z}_{\geqslant 0}$: $$\int_0^{2\pi} \cos^{p} \phi \sin^q \phi \mathrm{d}\phi = 2 \cos^2\left( \frac{\pi p}{2} \right) \cos^2\left( \frac{\pi q}{2} \right) \operatorname{Beta}\left(\frac{1+p}{2}, \frac{1+q}{2} \right) \tag{2}$$ notice that the right hand side vanishes whenever $p$ or $q$ are odd.

We now apply binomial theorem to the following and use eq. $(2)$ $$\int_0^{2\pi} \left(\cos \phi - \sin\phi\right)^{k_1} \left(\cos \phi + \sin\phi\right)^{k_2} \left((1+\kappa)\cos \phi + (1-\kappa)\sin\phi\right)^{k_3} \mathrm{d}\phi$$ Combining these results in the following Sum representation on your integral:

binoms[n_, k_] := Binomial[n, k] Binomial[(n + k - 1)/2, n];

sumint[n_Integer, kappa_, mu_] :=
Sum[If[EvenQ[k1 + k2 + k3],
binoms[n, k1] binoms[n, k2] binoms[n,
k3] Beta[(k1 + k2 + k3)/2 + 1, 1/2] 2^(3 n + 1 - k1 - k2 - k3)
1/(1 + 2 kappa mu + kappa^2)^(k3/2)
If[EvenQ[m1 + m2 + m3],
Binomial[k1, m1] Binomial[k2, m2] Binomial[k3, m3] (-1)^
m1 (1 + kappa)^(k3 - m3) (1 - kappa)^
m3 Beta[(1 + k1 - m1 + k2 - m2 + k3 - m3)/2, (1 + m1 + m2 + m3)/
2], 0], 0], {k1, Mod[n, 2], n, 2}, {k2, Mod[n, 2], n, 2}, {k3,
Mod[n, 1], n, 2}, {m1, 0, k1}, {m2, 0, k2}, {m3, 0, k3},
Method -> "Procedural"]


Verification:

integ[n_, kappa_, mu_] :=
LegendreP[n, 1/2 (Cos[p] - Sin[p]) Sin[t]] LegendreP[n,
1/2 (Cos[p] + Sin[p]) Sin[t]] LegendreP[
n, (((1 + kappa) Cos[p] + (1 - kappa) Sin[p]) Sin[t])/(
2 Sqrt[1 + 2 kappa mu + kappa^2])] Sin[t]


First check against Integrate for low values of n:

In:= mapint[e_Plus] :=
Integrate[#, {t, 0, Pi}, {p, 0, 2 Pi}] & /@ e;
mapint[e_] := Integrate[e, {t, 0, Pi}, {p, 0, 2 Pi}]

In:= Table[
sumint[n, ka, mu] -
mapint[Expand[integ[n, ka, mu], _Sin | _Cos]], {n, 0, 3}] //
Simplify // AbsoluteTiming

Out= {25.681514, {0, 0, 0, 0}}


Here are timings:

In:= Table[{n, First[AbsoluteTiming[sumint[n, ka, mu];]]}, {n, 1,
12}]

Out= {{1, 0.}, {2, 0.}, {3, 0.010000}, {4, 0.050001}, {5,
0.070001}, {6, 0.270005}, {7, 0.420008}, {8, 1.180024}, {9,
1.510030}, {10, 3.510070}, {11, 4.340087}, {12, 9.130183}}


Evaluation for $n=20$ takes 2 minutes on a real fast machine:

In:= AbsoluteTiming[r20 = sumint[20, ka, mu]; ]

Out= {120.752106, Null}

In:= AbsoluteTiming[r20s=Simplify[r20];]

Out= {0.148191, Null}

In:= LeafCount[r20s]

Out= 385


The expression is not fit for display here, as coefficients are some 20-digit integers. I can post the result if requested.

• Brilliant! Thanks a lot. Nov 15, 2012 at 20:34
• I tried my method on a fast machine, it ran all night for n=18. Nov 16, 2012 at 7:36
• what does Method->"Procedural" do in the sum? Didn't find it in the documentation. Nov 16, 2012 at 7:40
• @chris On the Sum's ref-page click on the 'More information' buttom. Find the table of methods. "Procedural" is listed there. It's effect is that Sum just adds terms, does not attempt any symbolic algorithms. Existence of this Method is motivated by the question, is it more efficient to sum $10^9$ terms of $\sum\limits_{k=1}^n k$, or compute the closed form formula and substitute $n=10^9$. Nov 16, 2012 at 12:47