0
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ What geometric meaning has the $\cot\alpha(\sin\theta(1-\cos\alpha)+\cos\theta \sin\alpha)$ angle? $\endgroup$
    – yarchik
    Commented Nov 25, 2023 at 11:38
  • $\begingroup$ 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. $\endgroup$
    – Rosstopher
    Commented Nov 25, 2023 at 12:00
  • $\begingroup$ 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. $\endgroup$
    – Rosstopher
    Commented Nov 25, 2023 at 12:08
  • $\begingroup$ Something like this math.stackexchange.com/q/4807073/435814 ? $\endgroup$
    – yarchik
    Commented Nov 25, 2023 at 12:10
  • $\begingroup$ 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. $\endgroup$
    – Rosstopher
    Commented Nov 25, 2023 at 12:14

1 Answer 1

0
$\begingroup$

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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