Re: SphericalHarmonics strange behavior
- To: mathgroup at smc.vnet.net
- Subject: [mg100985] Re: [mg100961] SphericalHarmonics strange behavior
- From: Luke Bloy <lbloy at seas.upenn.edu>
- Date: Fri, 19 Jun 2009 20:46:19 -0400 (EDT)
- References: <20090619072555.HPEMD.390879.imail@eastrmwml39>
Hi Bob, With some additional investigation yesterday I tracked the problem to the LegendreP and while i agree there is also a precision issue(which is discussed in some other bugs online), there seems to be an actual bug in the legendreP code dealing with negative m values. This bug is exposed when working with machine precision numbers. For example. LegendreP[1, 1, Cos[Pi/7]] // N Returns -0.433884 LegendreP[1, -1, Cos[Pi/7]] // N should be -1/2 LegendreP[1, 1, Cos[Pi/7]] = 0.216942 This is what is returned by LegendreP[1, -1, Cos[Pi/7]] // N However, LegendreP[1, -1, Cos[Pi/7.]] // N returns -0.216942 So it seems like the (-1)^m sign convention is not present in the machine precision routine that is run. -Luke Bob Hanlon wrote: > It's a precision issue > > diffY[l_, m_, t_, p_] := > Module[{val1, val2}, > val1 = SphericalHarmonicY[l, m, t, p]; > val2 = Sqrt[(2*l + 1) (l - m)!/(4*Pi (l + m)!)]* > LegendreP[l, m, Cos[t]]*E^(I*m*p); > Abs[val1 - val2]]; > > th1 = 1.324; > ph1 = 5.231; > > {diffY[1, 1, th1, ph1], diffY[1, -1, th1, ph1]} // > Chop > > {0,0.670051} > > th1 = 1.324`20; > ph1 = 5.231`20; > > {diffY[1, 1, th1, ph1], diffY[1, -1, th1, ph1]} // > Chop > > {0,0} > > Your definition would need to include an option for setting the working precision. > > > Bob Hanlon > > ---- "lbloy at seas.upenn.edu" <lbloy at seas.upenn.edu> wrote: > > ============= > Hi all, > > working with some spherical harmonics and have some questions about the code as implemented in mathematica. > > based on http://mathworld.wolfram.com/SphericalHarmonic.html > > I put together the following code to see what is implemented > > th1 = 1.324; > ph1 = 5.231; > > diffY[l_, m_, t_, p_] := > Module[{val1, val2}, val1 = SphericalHarmonicY[l, m, t, p]; > val2 = > Sqrt[((2*l + 1)*(l - m)!)/((4*Pi)*(l + m)!)]* > LegendreP[l, m, Cos[t]]*E^(I*m*p); > Print[{val1, val2, Abs[val1 - val2]}]; > Return[Abs[val1 - val2]]; ]; > > diffY[1, 1, th1, ph1] > > diffY[1, -1, th1, ph1] > > Clearly there is some sign error somewhere I was wondering if anyone had any suggestions. > > thanks > Luke > > >