Re: Roots of a Jacobi polynomial
- To: mathgroup at smc.vnet.net
- Subject: [mg120543] Re: Roots of a Jacobi polynomial
- From: Emu <samuel.thomas.blake at gmail.com>
- Date: Thu, 28 Jul 2011 07:54:31 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <j0ool2$kr9$1@smc.vnet.net>
On Jul 27, 8:18 pm, Bob Hanlon <hanl... at cox.net> wrote:
> Don't use machine precision
>
> eqn1 = JacobiP[25, -1/2, -1/2, x] == 0;
>
> soln1 = NSolve[eqn1, x, WorkingPrecision -> 15];
>
> Length[soln1]
>
> 25
>
> And @@ (eqn1 /. soln1)
>
> True
>
> And @@ (-1 < # < 1 & /@ (x /. soln1))
>
> True
>
> eqn2 = JacobiP[25, -0.5`20, -0.5`20, x] == 0;
>
> soln2 = NSolve[eqn2, x, WorkingPrecision -> 15];
>
> soln1 == soln2
>
> True
>
> Bob Hanlon
>
> ---- Jacopo Bertolotti <J.Bertolo... at utwente.nl> 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
Here's one way to get them in Alpha.
http://www.wolframalpha.com/input/?i=numerical+solutions+to+JacobiP%5B25%2C+-1%2F2%2C+-1%2F2%2C+x%5D+%3D%3D+0
Sam