Re: Re: Problem with the GaussianQuadratureWeights[n,a, b]

• To: mathgroup at smc.vnet.net
• Subject: [mg26465] Re: [mg26444] Re: Problem with the GaussianQuadratureWeights[n,a, b]
• From: Jens-Peer Kuska <kuska at informatik.uni-leipzig.de>
• Date: Thu, 21 Dec 2000 22:46:49 -0500 (EST)
• Organization: Universitaet Leipzig
• References: <91pjn3\$5vp@smc.vnet.net> <200012210651.BAA08643@smc.vnet.net> <p05010404b667d7cf36e3@[128.165.240.102]>
• Sender: owner-wri-mathgroup at wolfram.com

```Hi,

GaussianQuadrature[] computes the "weights" by the evaluation of
LegendreP[]. It is a well known fact that the numerical computation of
LegendreP[n,x] for high precision x (>16) is broken in Mathematica 4.0.x

You can check this by

In[]:=LegendreP[3, 0.5]
Out[]=-0.4375

In[]:=LegendreP[3, 0.5`30]
Out[]=-0.07291666666666666666666666667

You can fix this by adding

Unprotect[LegendreP]

LegendreP[n_, x_?NumericQ] /; Precision[x] > 16 :=
Module[{u, p}, p = LegendreP[n, u]; p /. u -> x]

Protect[LegendreP]

and now

gives the correct result.

because I have the feature-fix in my
init.m since version 4.0 came out.

You can still compute the weights by

GWeights1[n_Integer, a_, b_, prec_Integer] :=
Module[{x, xi, q, m},
xi = Sort[N[Solve[LegendreP[n, x] == 0, x], prec]] /.
Complex[r_, i_] :> r;
q = Integrate[LegendreP[n - 1, x]*LegendreP[n - 1, x], {x, -1, 1}];
wfac = {x, q/(LegendreP[n - 1, x]*D[LegendreP[n, x], x])} /. xi;
q = Plus @@ (Last /@ wfac);
{((b - a)*#[[1]] + (a + b))/2, 2*#[[2]]/q} & /@ wfac
]

or

GWeights2[n_Integer, a_, b_, prec_Integer] :=
Module[{x, xi, q, m},
xi = Sort[N[Solve[LegendreP[n, x] == 0, x], prec]] /.
Complex[r_, i_] :> r;

m = Table[LegendreP[i, x], {i, 0, n - 1}] /. xi;
{((b - a)*#[[1]] + (a + b))/2, (b - a)*#[[2]]/2} & /@
Transpose[{ x /. xi,
LinearSolve[Transpose[m],
Prepend[Table[0, {n - 1}],
Integrate[LegendreP[0, x], {x, -1, 1}]]]}]
]

both function evaluate LegendreP[] with symbolic x (that
work correct in version 4.0.x) and
replace it later by the numerical values.

Hope that helps
Jens

Cheng Liu wrote:
>
> Thank you Jens,
>
>         By setting the precision higher than the default (16) solve the
> problem described in my previous post.  But another problem appears,
>
> In[9]:=
> Out[9]=
> \!\({{0.0254460438286207377369051579760743056`24.7499,
>        0.08184679159616065512531181311972918224`22.9004}, \
> {0.1292344072003027800680676133596057328`25.5629,
>        0.176800361948499354889282418503202`22.9376}, \
> {0.2970774243113014165466967939615192127788252`26.1861,
>        0.241352841905111882955195047583155`22.9492}, {1\/2,
>        9.1004562140604214415812457652156224`23.0189*^-9}, \
> {0.70292257568869858345330320603848078722`26.5602,
>        0.241352841905111882955195047583156`22.9064}, \
> {0.8707655927996972199319323866403942672`26.3914,
>        0.176800361948499354889282418500441`22.3581}, \
> {0.9745539561713792622630948420239256944`26.3331,
>        0.081846791596160655125311813120987`22.4951}}\)
>
> Note the point at x=1/2, w is almost zero, which should be around 0.2.
>
> Another observation is that GaussianQuadratureWeights[7, 0, 1, 20] and