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

*To*: mathgroup at smc.vnet.net*Subject*: [mg27030] Re: Problem with the GaussianQuadratureWeights[n, a, b]*From*: Paul Cally <cally at kronos.maths.monash.edu.au>*Date*: Thu, 1 Feb 2001 03:00:20 -0500 (EST)*Organization*: Monash University*References*: <91pjn3$5vp@smc.vnet.net> <T1z06.128641$b16.480415@ralph.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Jens-Peer Kuska wrote: > Hi, > > setting the precision will help and > > GaussianQuadratureWeights[42, 0, 1, 32] > > work fine. > > Cheng Liu wrote: > > > > Hi, > > > > I accidently run into this strange behavior. I am using the > > package NumericalMath`GaussianQuadrature` and the function > > GaussianQuadratureWeights[n, a, b] to obtain the abscissas and weights > > for the n-point Gaussian quadrature formula on the interval (a, b). > > The specific interval in my problem is (0, 1), I noticed that > > when 40 <= n <= 50, if n is odd, the function gives correct answer, > > but when n is even, the function gives the following error message: > > > > In[43]:= > > GaussianQuadratureWeights[42, 0, 1] > > \!\(Power::"infy" \(\(:\)\(\ \)\) > > "Infinite expression \!\(1\/0.`\) encountered."\) > > \!\(Power::"infy" \(\(:\)\(\ \)\) > > "Infinite expression \!\(1\/0.`\) encountered."\) > > \!\(Power::"infy" \(\(:\)\(\ \)\) > > "Infinite expression \!\(1\/0.`\) encountered."\) > > General::stop : > > Further output of Power::infy will be suppressed during this calculation. > > f::indet : > > Indeterminate expression - 1672.96 + ComplexInfinity + > > ComplexInfinity + \[LeftSkeleton]6\[RightSkeleton] + ComplexInfinity + > > ComplexInfinity encountered. > SNIPP SNAPP > > Thanks. > > > > Cheng > > -- 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/ | +--------------------------------------------------------------------------+