Partial sum of Legendre polynomials and numerical error

Question

I am trying to compute the partial sums of the Legendre polynomials: $\sum_{l \geq 0} P_l(cos \gamma)$ of the first kind. The full sum are known to converge and I wanted to check that the error goes like $1/\sqrt{l}$. I wrote down three ways to compute the first 71 terms in the error. Unexpectedly, one of those (myerror2) takes a long time to compute and furthermore goes completely wild beyond $l \sim 50$. It actually blows up extremely rapidly for larger $l$. The other ones (myerror and myerror3) behave well and expected with the theory. What is going on?

exact[\[Gamma]_] := 1/2 Csc[\[Gamma]/2];
(*Sum[LegendreP[l,Cos[\[Gamma]]],{l,0,Infinity}]*)
my\[Gamma] = Pi/4.2;
myerror[\[Gamma]_] := 
  Table[Sum[LegendreP[l, Cos[\[Gamma]]], {l, 0, lmax}], {lmax, 1, 
     70}] - exact[\[Gamma]];
myerror2         = 
  Table[Sum[LegendreP[l, Cos[\[Gamma]]], {l, 0, lmax}], {lmax, 1, 
      70}] - exact[\[Gamma]] /. \[Gamma] -> my\[Gamma];
myerror3         = 
  Table[Sum[LegendreP[l, Cos[my\[Gamma]]], {l, 0, lmax}], {lmax, 1, 
     70}] - exact[my\[Gamma]];
Show[ListLogLogPlot[Abs[myerror2], PlotStyle -> Red, PlotMarkers -> "X"],
 ListLogLogPlot[Abs[myerror3], PlotStyle -> Black, PlotMarkers -> "<>"],
 ListLogLogPlot[Abs[myerror[my\[Gamma]]]],
 LogLogPlot[1.5/lmax^(1/2), {lmax, 1, 70}]]

As user Bob Hanlon suggested, using arbitrary-precision, instead of machine precision:

my\[Gamma] = N[10 Pi/42, 15];

allows for myerror2 to coincide with the other result, provided that the precision set is large enough. However, one does not seem to need to do that for myerror and myerror3. Furthermore, the computation time for myerror2 is still noticeably longer than the other two.

Answer 1

The point is, that in myerror2 Sum first evaluates the sum analytically ( that's why it takes longer than for myerror or myerror3) for the parameter gamma and then substitutes this gamma by mygamma, which has only MachinePrecision ( on my PC about 16).

Since the difference of the values of the summands (ranging from 0.09 to ...10^18) are greater than MachinePrecision, you get total wrong sum.

If you prevent Sum from immediate evaluation with Hold, ReplaceAll is done first (like in myerror and myerro3) and you get the wright result.

Of course you can solve the problem using higher Precision for mygamma.

mygamma = Pi/4.2;

l70 = LegendreP[70, Cos[gamma]]

(*  result not shown   *)

Sum[LegendreP[l, Cos[gamma]], {l, 70, 70}] /. gamma -> mygamma

(*   896.   *)

Hold[Sum][LegendreP[l, Cos[gamma]], {l, 70, 70}] /. 
   gamma -> mygamma // ReleaseHold

(*   -0.0126817   *)

ll = (List @@ Expand[l70]) /. gamma -> mygamma // N

(*   {-0.0950255, 126.892, -28206.9, 2.50097*10^6,  -1.18266*10^8, 
3.45865*10^9, -6.84288*10^10, 9.72619*10^11, -1.03659*10^13, 
8.55202*10^13, -5.59691*10^14, 2.96201*10^15, -1.28719*10^16, 
4.65028*10^16, -1.41075*10^17, 3.62314*10^17, -7.92911*10^17, 
1.48635*10^18, -2.39614*10^18, 3.33165*10^18, -4.00296*10^18, 
4.15969*10^18, -3.73806*10^18, 2.90146*10^18, -1.94063*10^18, 
1.11434*10^18, -5.46422*10^17, 2.27146*10^17, -7.926*10^16, 
2.29062*10^16, -5.38257*10^15, 1.00186*10^15, -1.42069*10^14, 
1.44144*10^13, -9.31671*10^11, 2.88157*10^10}   *)

Total@ll

(*   896.   *)