MathGroup Archive 2001

[Date Index] [Thread Index] [Author Index]

Search the Archive

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/ |
+--------------------------------------------------------------------------+





  • Prev by Date: Re: Re: Q: Factor with Polynominals?
  • Next by Date: Re: Why NestWhile?
  • Previous by thread: RE: Appending to Lists
  • Next by thread: Re: Problem with the GaussianQuadratureWeights[n, a, b]