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

• To: mathgroup at smc.vnet.net
• Subject: [mg27037] Re: Problem with the GaussianQuadratureWeights[n, a, b]
• From: Jens-Peer Kuska <kuska at informatik.uni-leipzig.de>
• Date: Thu, 1 Feb 2001 03:00:24 -0500 (EST)
• Organization: Universitaet Leipzig
• References: <91pjn3\$5vp@smc.vnet.net> <T1z06.128641\$b16.480415@ralph.vnet.net> <3A77A287.9CBF49CC@kronos.maths.monash.edu.au>
• Sender: owner-wri-mathgroup at wolfram.com

```Hi,

I had several private e-mails with the original poster, because the
mathgroup
is a bit slow.

I have a bug fix in my init.m, it looks like:

If[\$VersionNumber>=4.,

Unprotect[LegendreP]

LegendreP[lm__Integer, x_Real] /; Precision[x] > 16 :=
Module[{z}, LegendreP[lm, z] /. z -> x]

Protect[LegendreP]
]

I have posted this bug fix when Mathematica 4.0 came out
into the mathgroup. It evaluate LegendreP[] with symbolic
x (that is correct in version 4) and finaly replace the
symbolic argument with the high precision numerical argument.
It is a bit slower than buildin but better than the trash
produced by Mathematica 4.x

In my first reply to the poster I had simply forgotten that
my systems have the patch.

package works without any modification.

The other solution is to compute the weights with

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}]]]}]
]

none of the two funcions above are affected by the LegengreP[]
bug.

Hope that helps
Jens

Paul Cally wrote:

>  Aah, but only if N is even!! There is a terrible bug in Mathematica 4.0 and 4.1,
> where LegendreP gives a wildly wrong answer of the precision of the argument
> is higher than \$MachinePrecision. This really stuffs up the Gaussian quadrature
> weights if one of the nodal points is zero (ie, if N is odd). If N is even, it gets
> them
> all wrong by the same amount, but then fortuitously scales all the weights back to
> the correct values. Here are the details:
>
> Mathematica 4.0 for Linux
> Copyright 1988-1999 Wolfram Research, Inc.
>  -- Motif graphics initialized --
>
> In[1]:= LegendreP[5,2,N[1/3]]
>
> Out[1]= -10.3704
>
> In[2]:= LegendreP[5,2,N[1/3,20]]
>
> Out[2]= -0.17283950617283951
>
> In[3]:= %%/%
>
> Out[3]= 60.
>
> ===================================
>
> LegendreP[n,m,x] appears to work properly if x is a machine precision
> number, but if it is higher precision,
> it gives a result which is out by the factor n!/m!
>
> Please note that this also causes the standard package
> GaussianQuadratures to give wrong results for the
> quadrature weights if n is odd and precision is set to greater than
> \$MachinePrecision.
>
>
>
> Out[2]= {{-0.774597, 0.555556}, {0, 0.888889}, {0.774597, 0.555556}}
>
>
> Out[3]= {{-0.77459666924148337703585, 0.98360655737704918033},
>
> >    {0, 0.0327868852459016393443},
>
> >    {0.77459666924148337703585, 0.98360655737704918033}}
>
> ======================================================
>
> fixed in 4.1. As a temporary workaround, I define a function legendreP which
> calls LegendreP, but then multiplies by the required factor to undo the bug.
> You will also need to pull sections out of NumericalMath`GaussianQuadrature`
> and modify them to use legendreP.
>
> Paul Cally
>
> --
>
> +--------------------------------------------------------------------------+
> |Assoc Prof Paul Cally            |    Ph:  +61 3 9905-4471                |
> |Dept of Mathematics & Statistics |    Fax: +61 3 9905-3867                |
> |Monash University                |    paul.cally at sci.monash.edu.au        |
> |PO Box 28M, Victoria 3800        |                                        |
> |AUSTRALIA                        | http://www.maths.monash.edu.au/~cally/ |
> +--------------------------------------------------------------------------+

```

• Prev by Date: Re: Why NestWhile?
• Next by Date: NetCDF Files in Mathematica?