Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Roots of a Jacobi polynomial

• To: mathgroup at smc.vnet.net
• Subject: [mg120544] Re: Roots of a Jacobi polynomial
• From: Jacopo Bertolotti <J.Bertolotti at utwente.nl>
• Date: Thu, 28 Jul 2011 07:54:42 -0400 (EDT)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com
• References: <201107261107.HAA09296@smc.vnet.net> <4E2EFAAD.60403@wolfram.com>

On 07/26/2011 07:34 PM, Daniel Lichtblau wrote:
> 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
>
Thank you Daniel. And thank you also to everyone who wrote me 
suggestions (both here and in private). The combination of not using 
machine number in the input and increasing the WorkingPrecision seems to 
work fine. Apparently just using one of the two tricks do not yield any 
sensible result.

Jacopo

• Follow-Ups:
  • Re: Roots of a Jacobi polynomial
    • From: Richard Kendall-Smith <richpks@gmail.com>

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