Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2005
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: LegendreP -- up to which order?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg53573] Re: LegendreP -- up to which order?
  • From: Maxim <ab_def at prontomail.com>
  • Date: Tue, 18 Jan 2005 05:08:30 -0500 (EST)
  • Organization: MTU-Intel ISP
  • References: <cscicl$5lp$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

On Sun, 16 Jan 2005 02:09:24 +0000 (UTC), <gert.lindner at snafu.de> wrote:

>
> I compute the surface potentials generated by dipols within a  
> homogeneous sphere
> by using the associated Legendre Polynomials of first kind.
>
> The symmetric results caused by the angle positions are good for an  
> order < 60
> (Mathematica 5.1, 64 Bit optimized). The order is indicated in the cyan
> parameter field by 'm' (over sigma).
> *******************
>  The question is: What do you think about 'm' - how big can (should) I  
> set the
>  value for 'm'? (over sigma)
> There are failures (outliers for some angles) in the case m>65.
> Please help me (can you post it in the group for me?)
> Thanks
>  Gert Lindner

LegendreP[50, 1, x] equals Sqrt[(-1 - x)/(-1 + x)] times a polynomial of  
degree 50. For your values of phi and theta this is a very rapidly  
oscillating function, for which truncation errors become significant:

In[1]:=
LegendreP[50, 1, Cos[32*Degree]] // N
LegendreP[50, 1, Cos[32*Degree]] // N[#, 20]&

Out[1]=
-0.36566706

Out[2]=
-5.8185531335990549641

One way is to use arbitrary-precision numbers in the input:

\[CurlyTheta] = {92, 92, 115, 92, 60, 46, 60, 92, 115, 71, 32, 32, 71,  
115, 92, 46, 0, 46, 92, 115, 71, 32, 32, 71, 115, 92, 60, 46, 60, 92, 115,  
92, 92}*Degree;

\[CurlyPhi] = {108, 72, 145, 144, 129, 90, 51, 36, 35, 159, 135, 45, 21,  
180, 180, 180, 0, 0, 0, 0, 201, 225, 315, 339, 216, 216, 231, 270, 309,  
324, 324, 252, 288}*Degree;

\[CapitalPhi][m_] := (-Subscript[P,  
x])*(Cos[\[CurlyPhi]]/(4*Pi*\[Sigma]*R^2))*Sum[((2*n + 1)/n)*(Subscript[z,  
0]/R)^(n - 1)*LegendreP[n, 1, Cos[\[CurlyTheta]]], {n, 1, m}];

Block[{m, R, \[Sigma], z, P},
   m = 70; R = 0.0912`20; \[Sigma] = 1.`20;
   z /: Subscript[z, 0] = 0.075`20;
   P /: Subscript[P, x] = 50.`20*10^-9;
   ListPlot[\[CapitalPhi][m],
     PlotRange -> All, PlotStyle -> Hue[0.7],
     AxesLabel -> {"", "\[CapitalPhi]"}, PlotJoined -> True]
]

However, if the arguments of LegendreP are floating-point numbers, then  
Mathematica uses a different and more stable method of evaluating  
LegendreP (I suppose that numerical rules are used before algebraic  
rules), and LegendreP[50, 1, Cos[32*Degree] // N] works fine as well:

In[3]:=
LegendreP[50, 1, Cos[32*Degree] // N]

Out[3]=
-5.8185531

So you can just convert \[CurlyTheta] and \[CurlyPhi] to machine numbers:

\[CurlyTheta] = {92, 92, 115, 92, 60, 46, 60, 92, 115, 71, 32, 32, 71,  
115, 92, 46, 0, 46, 92, 115, 71, 32, 32, 71, 115, 92, 60, 46, 60, 92, 115,  
92, 92}*Degree // N;

\[CurlyPhi] = {108, 72, 145, 144, 129, 90, 51, 36, 35, 159, 135, 45, 21,  
180, 180, 180, 0, 0, 0, 0, 201, 225, 315, 339, 216, 216, 231, 270, 309,  
324, 324, 252, 288}*Degree // N;

or replace LegendreP[n, 1, Cos[\[CurlyTheta]]] with LegendreP[n, 1,  
1.*Cos[\[CurlyTheta]]] in the definition of \[CapitalPhi]:

\[CapitalPhi][m_] := (-Subscript[P, x])*(Cos[\[CurlyPhi]] /  
(4*Pi*\[Sigma]*R^2))*Sum[((2*n + 1)/n)*(Subscript[z, 0]/R)^(n -  
1)*LegendreP[n, 1, 1.*Cos[\[CurlyTheta]]], {n, 1, m}];

After this you needn't use arbitrary-precision numbers:

Block[{m, R, \[Sigma], z, P},
   m = 70; R = 0.0912; \[Sigma] = 1.;
   z /: Subscript[z, 0] = 0.075;
   P /: Subscript[P, x] = 50.*10^-9;
   ListPlot[\[CapitalPhi][m],
     PlotRange -> All, PlotStyle -> Hue[0.7],
     AxesLabel -> {"", "\[CapitalPhi]"}, PlotJoined -> True]
]

Another thing worth mentioning is that arbitrary-precision computations  
may fail if the numbers in the input have very low precision:

In[4]:=
Erf[Sqrt[-1 - I*z]] /. z -> 1`2*^-3
Erf[Sqrt[-1 - I z]] /. z -> 1`3*^-3

Out[4]=
0.0006 - 1.128*I

Out[5]=
0.00153 - 1.6504*I

Out[4] has incorrect significant digits.

Maxim Rytin
m.r at inbox.ru


  • Prev by Date: Re: about PATTERNS
  • Next by Date: Re: Mathematica Graphics output in an ASP.NET WebApplication
  • Previous by thread: LegendreP -- up to which order?
  • Next by thread: Re: LegendreP -- up to which order?