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