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. When the code above is loaded, the NumericalMath`GaussianQuadrature` 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. > > In[1]:= << NumericalMath`GaussianQuadrature` > > In[2]:= GaussianQuadratureWeights[3, -1, 1] > > Out[2]= {{-0.774597, 0.555556}, {0, 0.888889}, {0.774597, 0.555556}} > > In[3]:= GaussianQuadratureWeights[3, -1, 1,17] > > Out[3]= {{-0.77459666924148337703585, 0.98360655737704918033}, > > > {0, 0.0327868852459016393443}, > > > {0.77459666924148337703585, 0.98360655737704918033}} > > ====================================================== > > Wolfram were told about this months ago. It is very disappointing that it was not > 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/ | > +--------------------------------------------------------------------------+