Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Roots of a Jacobi polynomial

• To: mathgroup at smc.vnet.net
• Subject: [mg120513] Re: Roots of a Jacobi polynomial
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Wed, 27 Jul 2011 06:11:36 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com
• References: <201107261107.HAA09296@smc.vnet.net>

On 07/26/2011 06:07 AM, Jacopo Bertolotti wrote:
> Dear MathGroup,
> Mathematica implements Jacobi polynomials as JacobiP[n,a,b,x] where n is
> the order of the polynomial. As it can be checked plotting it a Jacobi
> polynomial has n real roots in the interval [-1,1] and it goes rapidly
> to infinity outside this interval (at least when both a and b are>-1).
> The problem arise when you try to find the roots of such a polynomial
> for a relatively high value of n. As an example the command
> NSolve[JacobiP[20, -0.5, -0.5, x] == 0, x] correctly returns
>
> {{x ->  -0.99702}, {x ->  -0.972111}, {x ->  -0.92418}, {x ->  -0.852398},
> {x ->  -0.760555}, {x ->  -0.649375}, {x ->  -0.522529}, {x ->  -0.382672},
> {x ->  -0.233449}, {x ->  -0.0784582}, {x ->  0.0784591}, {x ->  0.233445},
> {x ->  0.382684}, {x ->  0.522504}, {x ->  0.649423}, {x ->  0.760466}, {x
> ->  0.852539}, {x ->  0.924002}, {x ->  0.972267}, {x ->  0.996958}}
>
> while NSolve[JacobiP[25, -0.5, -0.5, x] == 0, x] gives
>
> {{x ->  -1.01869}, {x ->  -0.979859 - 0.0479527 I}, {x ->  -0.979859 +
> 0.0479527 I}, {x ->  -0.870962 - 0.070991 I}, {x ->  -0.870962 + 0.070991
> I}, {x ->  -0.71378 - 0.0505783 I}, {x ->  -0.71378 + 0.0505783 I}, {x ->
> -0.571283}, {x ->  -0.486396}, {x ->  -0.367829}, {x ->  -0.248377}, {x ->
> -0.125513}, {x ->  -0.0000434329}, {x ->  0.125442}, {x ->  0.2489}, {x ->
> 0.365644}, {x ->  0.496977}, {x ->  0.555743}, {x ->  0.717741- 0.0573399
> I}, {x ->  0.717741+ 0.0573399 I}, {x ->  0.87423- 0.0652273 I}, {x ->
> 0.87423+ 0.0652273 I}, {x ->  0.977876- 0.0422422 I}, {x ->  0.977876+
> 0.0422422 I}, {x ->  1.01494}}
>
> i.e. both complex roots and roots outside the [-1,1] interval.
> Substituting any of these values back into the polynomial easily show
> that these values are not roots at all. Also notice that using the
> command Root (e.g. Sort@Table[Root[JacobiP[25, -0.5, -0.5, x], i], {i,
> 1, 25}]) gives different but still wrong results.
> On a related note: NIntegrate sometimes gives wrong results when
> integrating function of the form (1-x)^a (1+x)^b f[x] and I feel this
> might related to the Gauss-Jacobi quadrature failing to retrieve the
> correct roots of a Jacobi polynomial.
>
> Do anyone have a solution for that?
>
> Thank you
>
> Jacopo

This is not a well conditioned problem. To do better you'll need to
allow NSolve to use higher precision by not giving machine numbers in 
the input.

In[102]:= solns = Sort[NSolve[JacobiP[25, -1/2, -1/2, x] == 0, x]]

Out[102]= {{x -> -0.998027}, {x -> -0.982287}, {x -> -0.951057}, {x \
-> -0.904827}, {x -> -0.844328}, {x -> -0.770513}, {x -> -0.684547}, \
{x -> -0.587785}, {x -> -0.481754}, {x -> -0.368125}, {x -> \
-0.24869}, {x -> -0.125333}, {x -> 0.}, {x -> 0.125333}, {x ->
    0.24869}, {x -> 0.368125}, {x -> 0.481754}, {x -> 0.587785}, {x ->
    0.684547}, {x -> 0.770513}, {x -> 0.844328}, {x -> 0.904827}, {x ->
     0.951057}, {x -> 0.982287}, {x -> 0.998027}}

These might seem imperfect when you check residuals (see below). In fact 
they are not at all bad. We first show that by comparing to higher 
precision versions that do validate fairly well as giving smallish 
residuals.

In[104]:= solns30 =
  Sort[NSolve[JacobiP[25, -1/2, -1/2, x] == 0, x,
    WorkingPrecision -> 30]]

Out[104]= {{x -> -0.998026728428271561952336806863}, {x -> \
-0.982287250728688681085641742865}, {x -> \
-0.951056516295153572116439333379}, {x -> \
-0.904827052466019527713668647933}, {x -> \
-0.844327925502015078548558063967}, {x -> \
-0.770513242775789230803009636396}, {x -> \
-0.684547105928688673732283357621}, {x -> \
-0.587785252292473129168705954639}, {x -> \
-0.481753674101715274987191502872}, {x -> \
-0.368124552684677959156947147493}, {x -> \
-0.248689887164854788242283746006}, {x -> \
-0.125333233564304245373118759817}, {x -> 0}, {x ->
    0.125333233564304245373118759817}, {x ->
    0.248689887164854788242283746006}, {x ->
    0.368124552684677959156947147493}, {x ->
    0.481753674101715274987191502872}, {x ->
    0.587785252292473129168705954639}, {x ->
    0.684547105928688673732283357621}, {x ->
    0.770513242775789230803009636396}, {x ->
    0.844327925502015078548558063967}, {x ->
    0.904827052466019527713668647933}, {x ->
    0.951056516295153572116439333379}, {x ->
    0.982287250728688681085641742865}, {x ->
    0.998026728428271561952336806863}}

This shows the two results are quite close.

In[106]:= Max[Abs[(x /. solns) - (x /. solns30)]]
Out[106]= 4.47615*10^-10

This shows the high precision solutions give reasonable residuals.

In[108]:= Max[Abs[JacobiP[25, -1/2, -1/2, x] /. solns30]]
Out[108]= 0.*10^-12

The machine precision results will not give good residuals. There are 
two possible (and related) reasons. One is that the conditioning is such 
that we are bound to lose many digits in assessing residuals (as seems 
to happen above). The other is that the form of JacobiP polynomials 
might just be inducing cancellation error.

Certainly this second item is part of the problem:

In[109]:= Max[Abs[Expand[JacobiP[25, -1/2, -1/2, x]] /. solns30]]
Out[109]= 0.*10^-21

Here are the corresponding machine precision evaluations of residuals.

In[112]:= Max[Abs[JacobiP[25, -1/2, -1/2, x] /. solns]]
Out[112]= 17.

In[113]:= Max[Abs[Expand[JacobiP[25, -1/2, -1/2, x]] /. solns]]
Out[113]= 9.31323*10^-9

Daniel Lichtblau
Wolfram Research

• References:
  • Roots of a Jacobi polynomial
    • From: Jacopo Bertolotti <J.Bertolotti@utwente.nl>