# Wrong numerical results from LegendreP

{Cos[Pi/180] // N, LegendreP[46, 0.9998476951563913], LegendreP[46, Cos[Pi/180]] // N}


give

{0.999848, 0.842007, 4.0625, 0.842007} (* Here is a typing error, the last element is not exist*)
{0.999848, 0.842007, 4.0625} (*Right version*)


And my version number of Mathematica is 10.3.0.0, platform is Microsoft Windows (64-bit).

Is it a foolish bug? The third element of the result list should be 0.842007, but my Mathematica gives a wrong result.

• bugs is reserved for when there is a consensus in the community that a bug exists (or confirmation from WRI). That said, something unusual is occurring, and it will require some extra investigation to see what it is. Also, why does your result have 4 elements in it? May 24, 2016 at 3:26
• Now, this is a question where the conclusions of this thread also apply. Have a look at N[List @@ LegendreP[46, Cos[Pi/180]]]; your problem is due to the fact that you are adding up a lot of big numbers to get a relatively tiny number, and that is not good numerics practice. LegendreP[46, N[Cos[Pi/180]]] avoids this since the numerical evaluation happens before LegendreP[] gets expanded. It is also well-known to those who know it that Legendre functions are difficult to numerically evaluate for arguments near $\pm 1$. May 24, 2016 at 4:16
• @rcollyer I give a list with 4 elements just for making comparisons. May 24, 2016 at 11:43
• @J.M. Is there some strategies to numerically evaluate Legendre functions for arguments near ±1? Can you give some advice? May 24, 2016 at 11:46
• What rcollyer was hinting at was that your input is a list of length 3, and your output is a list of length 4. Mathematica does not do that kind of magic. In any case, the brute-force method involves using arbitrary precision numbers and directly feeding such numbers into LegendreP[]. Clever methods are in the literature; search for them. May 24, 2016 at 11:54

From the LegendreP help page:

a bit of an extended comment, but i case anyone doesnt see the issue clearly, LegendreP[46, x] is a 46 order polynomial, with all even powers of x and alternating signs on the coefficients. We can separate out the positive and negative terms:

{neg, pos} = {
Total@MapIndexed[# x^(4  First@#2 - 4) &, #[[1 ;; ;; 2]]],
Total@MapIndexed[# x^(4 First@#2 - 2) &, #[[2 ;; ;; 2]]]} &@
CoefficientList[LegendreP[46, x], x][[1 ;; ;; 2]];

neg + pos == LegendreP[46, x] // Simplify


True

then compute the result at high precision:

 N[-neg /. x -> Cos[\[Pi]/180], 20]
N[pos /. x -> Cos[\[Pi]/180], 20]
Out[-1] - Out[-2]
`

1.5439246261069907522*10^16

1.5439246261069908364*10^16

0.842

You see the difference is in the 15th place, so we are on the hairy edge of machine precision yielding an approximate result depending on the actual order of the summation.

• That precisely is the "subtractive cancellation" I was alluding to. May 24, 2016 at 16:09