Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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[2] // 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.

share|improve this question
    
suggesting something to myself: may be use recurrence of LegendreP? –  chris Nov 15 '12 at 17:58
    
Is it possible at all to put in values for $\mu$ and $K$ and use NIntegrate? –  tkott Nov 15 '12 at 19:58
    
@tkott a symbolic solution was needed. Thanks for your interest. –  chris Nov 16 '12 at 7:31

1 Answer 1

up vote 13 down vote accepted

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[8]:= mapint[e_Plus] := 
  Integrate[#, {t, 0, Pi}, {p, 0, 2 Pi}] & /@ e;
mapint[e_] := Integrate[e, {t, 0, Pi}, {p, 0, 2 Pi}]

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

Out[10]= {25.681514, {0, 0, 0, 0}}

Here are timings:

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

Out[45]= {{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[14]:= AbsoluteTiming[r20 = sumint[20, ka, mu]; ]

Out[14]= {120.752106, Null}

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

Out[15]= {0.148191, Null}

In[16]:= LeafCount[r20s]

Out[16]= 385

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

share|improve this answer
    
Brilliant! Thanks a lot. –  chris Nov 15 '12 at 20:34
    
I tried my method on a fast machine, it ran all night for n=18. –  chris Nov 16 '12 at 7:36
1  
what does Method->"Procedural" do in the sum? Didn't find it in the documentation. –  chris Nov 16 '12 at 7:40
2  
@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$. –  Sasha Nov 16 '12 at 12:47

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.