Re: Roots of a Jacobi polynomial
- To: mathgroup at smc.vnet.net
- Subject: [mg120569] Re: Roots of a Jacobi polynomial
- From: Richard Kendall-Smith <richpks at gmail.com>
- Date: Fri, 29 Jul 2011 04:43:47 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <201107261107.HAA09296@smc.vnet.net> <4E2EFAAD.60403@wolfram.com> <201107281154.HAA04077@smc.vnet.net>
Hi,
I'm trying to solve a set of non-linear equations by doing the following:
eqlist = {101.74 == a + ((0.055/(479248/c))*(352343 + ((479248/c)*a) -
(158194*b/c))),
47205 + ((158194/c)*a) - ((61842/c)*b) == 68043.98 + (299.73*b),
460584 + (479248*(a^2)/(c^2)) - (158194*b*a/(c^2)) + (61842*(b^2)/(c^2))
== 21986 + (3576*c)}
NSolve[eqlist, {a, b, c}, Reals]
But it interprets c as Musical Notation. I have tried altering the variables
to x y and z but then it comes up with a Visual Form error. I'm new to
Mathematica but I don't see why this is happening because aren't we
specifying a, b, c as variables to be solved?
Thanks,
Richard
On Thu, Jul 28, 2011 at 5:54 AM, Jacopo Bertolotti
<J.Bertolotti at utwente.nl>wrote:
> 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
>
>
- References:
- Roots of a Jacobi polynomial
- From: Jacopo Bertolotti <J.Bertolotti@utwente.nl>
- Re: Roots of a Jacobi polynomial
- From: Jacopo Bertolotti <J.Bertolotti@utwente.nl>
- Roots of a Jacobi polynomial