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