0
{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.

tanghe2014
  • 195
  • 8
  • [tag: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? – rcollyer May 24 '16 at 03:26
  • 3
    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$. – J. M.'s missing motivation May 24 '16 at 04:16
  • @rcollyer I give a list with 4 elements just for making comparisons. – tanghe2014 May 24 '16 at 11:43
  • @J.M. Is there some strategies to numerically evaluate Legendre functions for arguments near ±1? Can you give some advice? – tanghe2014 May 24 '16 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. – J. M.'s missing motivation May 24 '16 at 11:54
  • Thank you all very much! Now I know it is not a bug, but a bad numeric practice. – tanghe2014 May 24 '16 at 11:56
  • Oh, sorry, the last element in the list is not exist. It is a typing error. – tanghe2014 May 24 '16 at 12:00
  • N[LegendreP[46, Cos[Pi/180]], $MachinePrecision] seems to give the correct answer. – xslittlegrass May 24 '16 at 14:34
  • @xslittlegrass, yes, the arbitrary precision evaluation seems to give mostly sound results, despite the threat of subtractive cancellation. But we are not always so lucky. – J. M.'s missing motivation May 24 '16 at 14:53

2 Answers2

3

From the LegendreP help page:

LegendreP

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Squigglyteeth
  • 106
  • 6
2

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.

george2079
  • 38,913
  • 1
  • 43
  • 110