Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Re: Problem with the GaussianQuadratureWeights[n,a, b]

• To: mathgroup at smc.vnet.net
• Subject: [mg26465] Re: [mg26444] Re: Problem with the GaussianQuadratureWeights[n,a, b]
• From: Jens-Peer Kuska <kuska at informatik.uni-leipzig.de>
• Date: Thu, 21 Dec 2000 22:46:49 -0500 (EST)
• Organization: Universitaet Leipzig
• References: <91pjn3$5vp@smc.vnet.net> <200012210651.BAA08643@smc.vnet.net> <p05010404b667d7cf36e3@[128.165.240.102]>
• Sender: owner-wri-mathgroup at wolfram.com

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

Cheng Liu wrote:
> 
> 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

• References:
  • Re: Problem with the GaussianQuadratureWeights[n, a, b]
    • From: Jens-Peer Kuska <kuska@informatik.uni-leipzig.de>