# Triple Summation Calculation Which Includes an Integral

I am a beginner in Mathematica. I have to calculate the following summation:$$C_3(\alpha)=\frac{1}{3}\sum_{j,k,l=3}^{\infty}(2j+1)(2k+1)(2l+1)(C_jC_k+C_kC_l+C_lC_j)f_{jkl}$$ where $$f_{jkl}=\frac{1}{4\pi}\int_{\theta=0}^\pi\int_{\phi=0} ^{2\pi}sin{\theta}d{\theta}d{\phi}P_j(cos\theta)P_k(cos(\theta-\alpha)P_l[cot\alpha(sin\theta(1-cos\alpha)+cos\theta sin\alpha)]$$, $$C_i=\frac{1}{i(i+1)}$$ and $$\alpha$$ runs from 0 to 120 degrees. Here $$P_j's$$ are Legendre polynomials. I have to calculate $$C_3({\alpha})$$ for values of $$\alpha$$ in the above-mentioned range. The summation runs up to infinity from the theoretical formula but I guess after j,k,l=10 maybe the value of the summation will saturate. I have written the code to calculate the summation for $$\alpha=20$$ degrees:

c[j_] := 1/(j*(j + 1))
c[l_] := 1/(l*(l + 1))
c[k_] := 1/(k*(k + 1))
f[j_, k_, l_, \[Alpha]_] :=
Assuming[\[Theta] \[Element] Reals && \[Phi] \[Element] Reals,
Integrate[
LegendreP[j, Cos[\[Theta]]]*
LegendreP[k, Cos[\[Theta] - \[Alpha] Degree]]*
LegendreP[l,
Cot[\[Alpha] Degree]*(Sin[\[Theta]]*(1 - Cos[\[Alpha] Degree]) +
Cos[\[Theta]]*Sin[\[Alpha] Degree])]*Sin[\[Theta]], {\[Theta],
0, \[Pi]}, {\[Phi], 0, 2 \[Pi]}]]
N[Sum[(2j+1)*(2k+1)*(2l+1)*(C_j*C_k+C_k*C_l+C_l*C_j)*f[j, k, l, 20 Degree]/(12*\[Pi]), {j, 3, 10}, {k, 3, 10}, {l, 3, 10}]]


But I am unable to get the value of the summation the code never gets executed and runs for a long time.

• What geometric meaning has the $\cot\alpha(\sin\theta(1-\cos\alpha)+\cos\theta \sin\alpha)$ angle? Commented Nov 25, 2023 at 11:38
• The argument of the Legendre polynomial is cotα(sinθ(1−cosα)+cosθsinα) which I found after some calculation. Actually, there are 4 unit vectors (\hat{n},\hat{a},\hat{b},\hat{c})in polar coordinates, and the arguments of Legendre polynomial is simply the dot product of \hat{n} with the other three unit vectors. From dot product cotα(sinθ(1−cosα)+cosθsinα) term arises. Commented Nov 25, 2023 at 12:00
• I meant to say the arguments of Legendre polynomial cos(θ) comes from the dot product of \hat{n} and \hat{a}, cos(θ-α) comes from the dot product of \hat{n} and \hat{b} and cotα(sinθ(1−cosα)+cosθsinα) comes from the dot product of \hat{n} and \hat{c}. \hat{n} makes polar angle θ and azimuthal angle ϕ in spherical polar coordinates. Commented Nov 25, 2023 at 12:08
• Something like this math.stackexchange.com/q/4807073/435814 ? Commented Nov 25, 2023 at 12:10
• Yes! Actually, that question you mentioned was posted by me. I assumed \hat{a} along the z-axis and found the modified expression as mentioned here. Commented Nov 25, 2023 at 12:14

There were typos in the code you posted as I noted in the comment above. You had terms like C_j, etc., which should be c[j], etc. Since you're interested in a numerical result, I replaced Integrate with NIntegrate, and N@Sum with NSum. You had three c[_] definitions; one will suffice since they're basically the same formula. I also removed the Degree unit while tinkering, and just divided alpha by 180 in your f definition; you can revert this if you wish. I just wanted to declutter as much as possible while trying to get this to work.

After making the changes, it ran for a while without returning an answer, so I just ran f[__] by itself inside Table for a few values. I received a number of messages indicating convergence difficulties for the default NIntegrate option GlobalAdaptive, so I tried LocalAdaptive to see if it would handle your integrals more gracefully, and it did. I didn't see any more messages when I ran my test. I set up a function sumjlk to allow me to try varying the limits on your sum, to determine how long it might take. It ran acceptable fast for a small case, so I bumped up the limits to match your original problem.

Here's the modified code:

c[jkl_] := 1/(jkl*(jkl + 1))
f[j_, k_, l_, \[Alpha]_] :=
NIntegrate[
LegendreP[j, Cos[\[Theta]]]*
LegendreP[k, Cos[\[Theta] - \[Alpha]]]*
LegendreP[l,
Cot[\[Alpha]]*(Sin[\[Theta]]*(1 - Cos[\[Alpha]]) +
Cos[\[Theta]]*Sin[\[Alpha]])]*Sin[\[Theta]], {\[Theta],
0, \[Pi]}, {\[Phi], 0, 2 \[Pi]},Method->"LocalAdaptive"];

sumjlk[jmax_,kmax_,lmax_] := NSum[(2j+1)*(2k+1)*(2l+1)*(c[j] c[k]+c[k]*c[l]+c[l]*c[j])*f[j, k, l, 20 Degree],
{j, 3, jmax}, {k, 3, kmax}, {l, 3, lmax}]/(12*\[Pi]);


And this is the result I obtained after roughly one minute:

sumjlk[10,10,10]


EDIT: I left out a factor of Pi in my earlier post when converting degrees to radians. After fixing and rerunning, this is the result.

(* 1.30704*)