Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Problem with the GaussianQuadratureWeights[n, a, b]

87 views
Skip to first unread message

Cheng Liu

unread,
Dec 20, 2000, 1:30:59 AM12/20/00
to
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.
ƒ::indet :
Indeterminate expression - 1672.96 + ComplexInfinity +
ComplexInfinity + \[LeftSkeleton]6\[RightSkeleton] + ComplexInfinity +
ComplexInfinity encountered.
Out[43]=
{{0.000800191, Indeterminate}, {0.00421136, Indeterminate}, {0.0103287,
Indeterminate}, {0.0191203, Indeterminate}, {0.0305382,
Indeterminate}, {0.0445201, Indeterminate}, {0.0609897,
Indeterminate}, {0.079857, Indeterminate}, {0.101019,
Indeterminate}, {0.12436, Indeterminate}, {0.149753,
Indeterminate}, {0.177058, Indeterminate}, {0.206128,
Indeterminate}, {0.236802, Indeterminate}, {0.268914,
Indeterminate}, {0.302288, Indeterminate}, {0.336742,
Indeterminate}, {0.372087, Indeterminate}, {0.408132,
Indeterminate}, {0.444677, Indeterminate}, {0.481526,
Indeterminate}, {0.518474, Indeterminate}, {0.555323,
Indeterminate}, {0.591868, Indeterminate}, {0.627913,
Indeterminate}, {0.663258, Indeterminate}, {0.697712,
Indeterminate}, {0.731086, Indeterminate}, {0.763198,
Indeterminate}, {0.793872, Indeterminate}, {0.822942,
Indeterminate}, {0.850247, Indeterminate}, {0.87564,
Indeterminate}, {0.898981, Indeterminate}, {0.920143,
Indeterminate}, {0.93901, Indeterminate}, {0.95548,
Indeterminate}, {0.969462, Indeterminate}, {0.98088,
Indeterminate}, {0.989671, Indeterminate}, {0.995789,
Indeterminate}, {0.9992, Indeterminate}}

When n > 50, the above behavior disappears. Of course, I didn't check
cases when n is really big (it's also not necessary in view of ordinary
numerical computation). The version of Mathematica I am using is
"4.0 for Power Macintosh (July 26, 1999)". Any explanations?

Thanks.

Cheng
--
====================================================
Cheng Liu, Ph.D.
MST-8, Structure/Property Relations
Materials Science and Technology Division
Los Alamos National Laboratory
Los Alamos, New Mexico 87545

Phone: 505-665-6892 (office), 505-667-9950 (lab)
Fax: 505-667-8021
email: cl...@lanl.gov
====================================================

Jens-Peer Kuska

unread,
Dec 21, 2000, 9:31:15 PM12/21/00
to
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
> --


Cheng Liu

unread,
Dec 21, 2000, 10:48:49 PM12/21/00
to
Thank you Jens,

By setting the precision higher than the default (16) solve the
problem described in my previous post. But another problem appears,

In[9]:=
GaussianQuadratureWeights[7, 0, 1, 20]
Out[9]=
\!\({{0.0254460438286207377369051579760743056`24.7499,
0.08184679159616065512531181311972918224`22.9004}, \
{0.1292344072003027800680676133596057328`25.5629,
0.176800361948499354889282418503202`22.9376}, \
{0.2970774243113014165466967939615192127788252`26.1861,
0.241352841905111882955195047583155`22.9492}, {1\/2,
9.1004562140604214415812457652156224`23.0189*^-9}, \
{0.70292257568869858345330320603848078722`26.5602,
0.241352841905111882955195047583156`22.9064}, \
{0.8707655927996972199319323866403942672`26.3914,
0.176800361948499354889282418500441`22.3581}, \
{0.9745539561713792622630948420239256944`26.3331,
0.081846791596160655125311813120987`22.4951}}\)

Note the point at x=1/2, w is almost zero, which should be around 0.2.

Another observation is that GaussianQuadratureWeights[7, 0, 1, 20] and
GaussianQuadratureWeights[7, 0, 1] give different answers to the weight
w. Any help would be appreciated. Thanks.

Cheng

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

--

Jens-Peer Kuska

unread,
Dec 21, 2000, 10:53:55 PM12/21/00
to
Hi,

GaussianQuadrature[] computes the "weights" by the evaluation of
LegendreP[]. It is a well known fact that the numerical computation of
LegendreP[n,x] for high precision x (>16) is broken in Mathematica 4.0.x

You can check this by

In[]:=LegendreP[3, 0.5]
Out[]=-0.4375


In[]:=LegendreP[3, 0.5`30]
Out[]=-0.07291666666666666666666666667

You can fix this by adding

Unprotect[LegendreP]

LegendreP[n_, x_?NumericQ] /; Precision[x] > 16 :=
Module[{u, p}, p = LegendreP[n, u]; p /. u -> x]

Protect[LegendreP]

to your init.m file.

and now

GaussianQuadratureWeights[7, -1, 1, 20]

gives the correct result.

In my first answer I don't think about it,
because I have the feature-fix in my
init.m since version 4.0 came out.

You can still compute the weights by

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

both function evaluate LegendreP[] with symbolic x (that
work correct in version 4.0.x) and
replace it later by the numerical values.

Hope that helps
Jens

Paul Cally

unread,
Feb 1, 2001, 3:00:26 AM2/1/01
to
Jens-Peer Kuska 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....@sci.monash.edu.au |
|PO Box 28M, Victoria 3800 | |
|AUSTRALIA | http://www.maths.monash.edu.au/~cally/ |
+--------------------------------------------------------------------------+


Jens-Peer Kuska

unread,
Feb 1, 2001, 3:10:26 AM2/1/01
to
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

0 new messages