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