Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2001
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

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


  • Prev by Date: Re: Why NestWhile?
  • Next by Date: NetCDF Files in Mathematica?
  • Previous by thread: Re: Problem with the GaussianQuadratureWeights[n, a, b]
  • Next by thread: Re: Why NestWhile?