[R] Solving 100th order equation

[Spencer Graves]

Consider the following: 

library(polynomial)
p100 <- polynomial(1:101)
s100 <- solve(p100)
p100r <- 101*poly.calc(s100, tol=.Machine$double.eps)
max(abs(coef(p100r)-(1:101)))
 550767713

      For order 2, it worked perfectly: 

(p2 <- polynomial(1:3))
(s2 <- solve(p2))
(p2r <- 3*poly.calc(s2))
max(abs(coef(p2r)-1:3))
0

      For order 50, it clearly had some problems: 

p50 <- polynomial(1:51)
s50 <- solve(p50)
p50r <- 3*poly.calc(s50)
max(abs(coef(p50r)-1:3))
max(abs(coef(p50r)-1:3))
[1] 2.823529

      To find the theoretical method, the ultimate reference is the 
code.  Reading "polynom:::solve.polynomial" reveals that this passes an 
appropriate matrix to 'eigen'. 

      Hope this helps. 
      Spencer

Shubha Vishwanath Karanth wrote:
> Was also wondering which theoretical method is used to solve this problem?
>
> Thanks,
> Shubha Karanth | Amba Research
> Ph +91 80 3980 8031 | Mob +91 94 4886 4510 
> Bangalore * Colombo * London * New York * San José * Singapore * www.ambaresearch.com
>
> -----Original Message-----
> From: Gabor Grothendieck [mailto:ggrothendieck at gmail.com] 
> Sent: Saturday, May 24, 2008 6:13 PM
> To: Peter Dalgaard
> Cc: Shubha Vishwanath Karanth; r-help at stat.math.ethz.ch; Duncan Murdoch
> Subject: Re: [R] Solving 100th order equation
>
> On Sat, May 24, 2008 at 8:31 AM, Peter Dalgaard
> <p.dalgaard at biostat.ku.dk> wrote:
>   
>> Shubha Vishwanath Karanth wrote:
>>     
>>> To apply uniroot I don't even know the interval values... Does numerical
>>> methods help me? Or any other method?
>>>
>>> Thanks and Regards,
>>> Shubha
>>>
>>> -----Original Message-----
>>> From: Duncan Murdoch [mailto:murdoch at stats.uwo.ca] Sent: Saturday, May 24,
>>> 2008 5:08 PM
>>> To: Shubha Vishwanath Karanth
>>> Subject: Re: [R] Solving 100th order equation
>>>
>>> Shubha Vishwanath Karanth wrote:
>>>
>>>       
>>>> Hi R,
>>>>
>>>>
>>>> I have a 100th order equation for which I need to solve the value for x.
>>>> Is there a package to do this?
>>>>
>>>>
>>>> For example my equation is:
>>>>
>>>>
>>>> (x^100 )- (2*x^99) +(10*x^50)+.............. +(6*x ) = 4000
>>>>
>>>>
>>>> I have only one unknown value and that is x. How do I solve for this?
>>>>
>>>>
>>>>         
>>> uniroot() will find one root.  If you want all of them, I don't know what
>>> is available.
>>>
>>> Duncan Murdoch
>>>
>>>       
>> polyroot() is built for this, but it stops at 48th degree polynomials, at
>> least as currently implemented. Not sure that it (or anything else) would be
>> stable beyond that limit. YACAS perhaps?
>>
>>     
>
> Unfortunately yacas does not seem to be able to handle it:
>
>   
>> library(Ryacas)
>> x <- Sym("x")
>> Solve((x^100 )- (2*x^99) +(10*x^50)+(6*x ) - 4000 == 0, x)
>>     
> [1] "Starting Yacas!"
> expression(list())
>
> Simpler one works ok:
>
>   
>> Solve(x^2 - 1, x)
>>     
> expression(list(x == 1, x == -1))
> This e-mail may contain confidential and/or privileged i...{{dropped:10}}
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>