[R] Solving 100th order equation
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?
BR, Shubha
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
This e-mail may contain confidential and/or privileged i...{{dropped:13}}
__
R-help@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.
Re: [R] Solving 100th order equation
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:[EMAIL PROTECTED]
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
>
>
>
>
>
>
> BR, Shubha
>
> 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
>
>
>
> This e-mail may contain confidential and/or privileged i...{{dropped:13}}
>
>
>
>
> __
> R-help@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.
>
This e-mail may contain confidential and/or privileged i...{{dropped:10}}
__
R-help@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.
Re: [R] Solving 100th order equation
Shubha Vishwanath Karanth wrote:
To apply uniroot I don't even know the interval values... Does numerical
methods help me? Or any other method?
How about plotting the function? If you write it as f(x) = x^100 - ...
- 6*x - 4000, then clearly f(0) is negative, and large enough x will be
positive. Just find
out what "large enough" means (it depends on the ...) and use that to
start uniroot().
Duncan Murdoch
Thanks and Regards,
Shubha
-Original Message-
From: Duncan Murdoch [mailto:[EMAIL PROTECTED]
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
BR, Shubha
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
This e-mail may contain confidential and/or privileged i...{{dropped:13}}
__
R-help@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.
This e-mail may contain confidential and/or privileged i...{{dropped:10}}
__
R-help@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.
__
R-help@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.
Re: [R] Solving 100th order equation
Shubha Vishwanath Karanth wrote:
To apply uniroot I don't even know the interval values... Does numerical
methods help me? Or any other method?
I forgot: we also have polyroot().
Duncan Murdoch
Thanks and Regards,
Shubha
-Original Message-
From: Duncan Murdoch [mailto:[EMAIL PROTECTED]
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
BR, Shubha
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
This e-mail may contain confidential and/or privileged i...{{dropped:13}}
__
R-help@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.
This e-mail may contain confidential and/or privileged...{{dropped:8}}
__
R-help@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.
Re: [R] Solving 100th order equation
Shubba,
It is hard to imagine why you would want to solve a 100th order
equation. You might want to make absolutely certain that you are
thinking about the correct problem. What are you trying to do? To do
what you want you will need a very large dataset. Perhaps you need to
fit some time of smoothing function (e.g. smoothing spline) to your
data.
John
John David Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
>>> "Shubha Vishwanath Karanth" <[EMAIL PROTECTED]> 5/24/2008
7:52 AM >>>
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:[EMAIL PROTECTED]
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
>
>
>
>
>
>
> BR, Shubha
>
> 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
>
>
>
> This e-mail may contain confidential and/or privileged
i...{{dropped:13}}
>
>
>
>
> __
> R-help@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.
>
This e-mail may contain confidential and/or privileged
i...{{dropped:10}}
__
R-help@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.
Confidentiality Statement:
This email message, including any attachments, is for th...{{dropped:6}}
__
R-help@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.
Re: [R] Solving 100th order equation
See ?polyroot. However, the numerical stability for such a high-order problem is almost certainly a serious issue. On Sat, 24 May 2008, 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:[EMAIL PROTECTED] 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 BR, Shubha 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 -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595__ R-help@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.
Re: [R] Solving 100th order equation
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:[EMAIL PROTECTED] 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? -- O__ Peter Dalgaard Øster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ R-help@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.
Re: [R] Solving 100th order equation
On Sat, May 24, 2008 at 8:31 AM, Peter Dalgaard
<[EMAIL PROTECTED]> 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:[EMAIL PROTECTED] 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))
__
R-help@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.
Re: [R] Solving 100th order equation
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:[EMAIL PROTECTED]
Sent: Saturday, May 24, 2008 6:13 PM
To: Peter Dalgaard
Cc: Shubha Vishwanath Karanth; [EMAIL PROTECTED]; Duncan Murdoch
Subject: Re: [R] Solving 100th order equation
On Sat, May 24, 2008 at 8:31 AM, Peter Dalgaard
<[EMAIL PROTECTED]> 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:[EMAIL PROTECTED] 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@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.
Re: [R] Solving 100th order equation
At 07:58 AM 5/24/2008, Duncan Murdoch 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? I forgot: we also have polyroot(). Duncan Murdoch Thanks and Regards, Shubha -Original Message- From: Duncan Murdoch [mailto:[EMAIL PROTECTED] 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? Finding all of the roots of a 100-th order polynomial on a computer is not an easy task. It would take move than 100 digits of precision to do so. A similar problem to finding all 100 eigenvalues of a 100x100 ill-conditioned matrix. The suggestion to use graphics to find small intervals localizing the roots of interest first would make the most sense, at least for the real roots. It is highly unlikely that all 100 roots are needed. Usually the ones of interest are the one with the largest real part, or the ones inside the unit circle, or the real roots, etc. Is this a homework problem? Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail: [EMAIL PROTECTED] Least Cost Formulations, Ltd.URL: http://lcfltd.com/ 824 Timberlake Drive Tel: 757-467-0954 Virginia Beach, VA 23464-3239Fax: 757-467-2947 "Vere scire est per causas scire" __ R-help@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.
Re: [R] Solving 100th order equation
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:[EMAIL PROTECTED]
Sent: Saturday, May 24, 2008 6:13 PM
To: Peter Dalgaard
Cc: Shubha Vishwanath Karanth; [EMAIL PROTECTED]; Duncan Murdoch
Subject: Re: [R] Solving 100th order equation
On Sat, May 24, 2008 at 8:31 AM, Peter Dalgaard
<[EMAIL PROTECTED]> 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:[EMAIL PROTECTED] 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@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.
__
R-help@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.
Re: [R] Solving 100th order equation
Robert A LaBudde wrote: At 07:58 AM 5/24/2008, Duncan Murdoch 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? I forgot: we also have polyroot(). Duncan Murdoch Thanks and Regards, Shubha -Original Message- From: Duncan Murdoch [mailto:[EMAIL PROTECTED] 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? Finding all of the roots of a 100-th order polynomial on a computer is not an easy task. It would take move than 100 digits of precision to do so. A similar problem to finding all 100 eigenvalues of a 100x100 ill-conditioned matrix. The suggestion to use graphics to find small intervals localizing the roots of interest first would make the most sense, at least for the real roots. It is highly unlikely that all 100 roots are needed. Usually the ones of interest are the one with the largest real part, or the ones inside the unit circle, or the real roots, etc. And uniroot, polyroot or solve.polynom should provide marvelous starting values. Is this a homework problem? Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail: [EMAIL PROTECTED] Least Cost Formulations, Ltd.URL: http://lcfltd.com/ 824 Timberlake Drive Tel: 757-467-0954 Virginia Beach, VA 23464-3239Fax: 757-467-2947 "Vere scire est per causas scire" __ R-help@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. __ R-help@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.
Re: [R] Solving 100th order equation
x_cartesian(0,
-1), x == complex_cartesian(cos(12 * pi/25), -sin(12 *
pi/25)), x == complex_cartesian(cos(23 * pi/50), -sin(23 *
pi/50)), x == complex_cartesian(cos(11 * pi/25), -sin(11 *
pi/25)), x == complex_cartesian(cos(21 * pi/50), -sin(21 *
pi/50)), x == complex_cartesian((root(5, 2) - 1)/4, -sin(2 *
pi/5)), x == complex_cartesian(cos(19 * pi/50), -sin(19 *
pi/50)), x == complex_cartesian(cos(9 * pi/25), -sin(9 *
pi/25)), x == complex_cartesian(cos(17 * pi/50), -sin(17 *
pi/50)), x == complex_cartesian(cos(8 * pi/25), -sin(8 *
pi/25)), x == complex_cartesian(cos(3 * pi/10), -sin(3 *
pi/10)), x == complex_cartesian(cos(7 * pi/25), -sin(7 *
pi/25)), x == complex_cartesian(cos(13 * pi/50), -sin(13 *
pi/50)), x == complex_cartesian(cos(6 * pi/25), -sin(6 *
pi/25)), x == complex_cartesian(cos(11 * pi/50), -sin(11 *
pi/50)), x == complex_cartesian(cos(pi/5), -sin(pi/5)),
x == complex_cartesian(cos(9 * pi/50), -sin(9 * pi/50)),
x == complex_cartesian(cos(4 * pi/25), -sin(4 * pi/25)),
x == complex_cartesian(cos(7 * pi/50), -sin(7 * pi/50)),
x == complex_cartesian(cos(3 * pi/25), -sin(3 * pi/25)),
x == complex_cartesian(cos(pi/10), -((root(5, 2) - 1)/4)),
x == complex_cartesian(cos(2 * pi/25), -sin(2 * pi/25)),
x == complex_cartesian(cos(3 * pi/50), -sin(3 * pi/50)),
x == complex_cartesian(cos(pi/25), -sin(pi/25)), x ==
complex_cartesian(cos(pi/50),
-sin(pi/50)), x == 1))
On Sat, May 24, 2008 at 8:56 AM, Shubha Vishwanath Karanth
<[EMAIL PROTECTED]> 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:[EMAIL PROTECTED]
> Sent: Saturday, May 24, 2008 6:13 PM
> To: Peter Dalgaard
> Cc: Shubha Vishwanath Karanth; [EMAIL PROTECTED]; Duncan Murdoch
> Subject: Re: [R] Solving 100th order equation
>
> On Sat, May 24, 2008 at 8:31 AM, Peter Dalgaard
> <[EMAIL PROTECTED]> 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:[EMAIL PROTECTED] 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...{{dropped:12}}
__
R-help@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.
Re: [R] Solving 100th order equation
> library(PolynomF) > x <- polynom() > p <- x^100 - 2*x^99 + 10*x^50 + 6*x - 4000 > z <- solve(p) > z [1] -1.0741267+0.000i -1.073-0.0680356i -1.073+0.0680356i -1.0655699-0.1354644i [5] -1.0655699+0.1354644i -1.0568677-0.2030274i -1.0568677+0.2030274i -1.0400346-0.2687815i ... [93] 1.0595174+0.2439885i 1.0746575-0.1721335i 1.0746575+0.1721335i 1.0828132-0.1065591i [97] 1.0828132+0.1065591i 1.0879363-0.0330308i 1.0879363+0.0330308i 2.000+0.000i > Now to check how good they are: > > range(Mod(p(z))) [1] 1.062855e-10 1.548112e+15 > Not brilliant, but not too bad. Bill Venables CSIRO Laboratories PO Box 120, Cleveland, 4163 AUSTRALIA Office Phone (email preferred): +61 7 3826 7251 Fax (if absolutely necessary): +61 7 3826 7304 Mobile: +61 4 8819 4402 Home Phone: +61 7 3286 7700 mailto:[EMAIL PROTECTED] http://www.cmis.csiro.au/bill.venables/ -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Gabor Grothendieck Sent: Sunday, 25 May 2008 10:09 AM To: Shubha Vishwanath Karanth Cc: [EMAIL PROTECTED]; Duncan Murdoch; Peter Dalgaard Subject: Re: [R] Solving 100th order equation Actually maybe I was premature. It does not handle the polynomial I tried it on in the example earlier in this thread but it does seem to work with the following very simple polynomials of order 100. At any rate it would not take long to try it on the real problem and see. > Solve(x^100 - 1, x) [1] "Starting Yacas!" expression(list(x == complex_cartesian(cos(pi/50), sin(pi/50)), x == complex_cartesian(cos(pi/25), sin(pi/25)), x == complex_cartesian(cos(3 * pi/50), sin(3 * pi/50)), x == complex_cartesian(cos(2 * pi/25), sin(2 * pi/25)), x == complex_cartesian(cos(pi/10), (root(5, 2) - 1)/4), x == complex_cartesian(cos(3 * pi/25), sin(3 * pi/25)), x == complex_cartesian(cos(7 * pi/50), sin(7 * pi/50)), x == complex_cartesian(cos(4 * pi/25), sin(4 * pi/25)), x == complex_cartesian(cos(9 * pi/50), sin(9 * pi/50)), x == complex_cartesian(cos(pi/5), sin(pi/5)), x == complex_cartesian(cos(11 * pi/50), sin(11 * pi/50)), x == complex_cartesian(cos(6 * pi/25), sin(6 * pi/25)), x == complex_cartesian(cos(13 * pi/50), sin(13 * pi/50)), x == complex_cartesian(cos(7 * pi/25), sin(7 * pi/25)), x == complex_cartesian(cos(3 * pi/10), sin(3 * pi/10)), x == complex_cartesian(cos(8 * pi/25), sin(8 * pi/25)), x == complex_cartesian(cos(17 * pi/50), sin(17 * pi/50)), x == complex_cartesian(cos(9 * pi/25), sin(9 * pi/25)), x == complex_cartesian(cos(19 * pi/50), sin(19 * pi/50)), x == complex_cartesian((root(5, 2) - 1)/4, sin(2 * pi/5)), x == complex_cartesian(cos(21 * pi/50), sin(21 * pi/50)), x == complex_cartesian(cos(11 * pi/25), sin(11 * pi/25)), x == complex_cartesian(cos(23 * pi/50), sin(23 * pi/50)), x == complex_cartesian(cos(12 * pi/25), sin(12 * pi/25)), x == complex_cartesian(0, 1), x == complex_cartesian(-cos(12 * pi/25), sin(12 * pi/25)), x == complex_cartesian(-cos(23 * pi/50), sin(23 * pi/50)), x == complex_cartesian(-cos(11 * pi/25), sin(11 * pi/25)), x == complex_cartesian(-cos(21 * pi/50), sin(21 * pi/50)), x == complex_cartesian(-((root(5, 2) - 1)/4), sin(2 * pi/5)), x == complex_cartesian(-cos(19 * pi/50), sin(19 * pi/50)), x == complex_cartesian(-cos(9 * pi/25), sin(9 * pi/25)), x == complex_cartesian(-cos(17 * pi/50), sin(17 * pi/50)), x == complex_cartesian(-cos(8 * pi/25), sin(8 * pi/25)), x == complex_cartesian(-cos(3 * pi/10), sin(3 * pi/10)), x == complex_cartesian(-cos(7 * pi/25), sin(7 * pi/25)), x == complex_cartesian(-cos(13 * pi/50), sin(13 * pi/50)), x == complex_cartesian(-cos(6 * pi/25), sin(6 * pi/25)), x == complex_cartesian(-cos(11 * pi/50), sin(11 * pi/50)), x == complex_cartesian(-cos(pi/5), sin(pi/5)), x == complex_cartesian(-cos(9 * pi/50), sin(9 * pi/50)), x == complex_cartesian(-cos(4 * pi/25), sin(4 * pi/25)), x == complex_cartesian(-cos(7 * pi/50), sin(7 * pi/50)), x == complex_cartesian(-cos(3 * pi/25), sin(3 * pi/25)), x == complex_cartesian(-cos(pi/10), (root(5, 2) - 1)/4), x == complex_cartesian(-cos(2 * pi/25), sin(2 * pi/25)), x == complex_cartesian(-cos(3 * pi/50), sin(3 * pi/50)), x == complex_cartesian(-cos(pi/25), sin(pi/25)), x == complex_cartesian(-cos(pi/50), sin(pi/50)), x == -1, x == complex_cartesian(-cos(pi/50), -sin(pi/50)), x == complex_cartesian(-cos(pi/25), -sin(pi/25)), x == complex_cartesian(-cos(3 * pi/50), -sin(3 * pi/50)), x == complex_cartesian(-cos(2 * pi/25), -sin(2 * pi/25)), x == complex_cartesian(-cos(pi/10), -((root(5, 2) - 1)/4)), x == complex_carte
Re: [R] Solving 100th order equation
Oops! Actually spectacularly bad. I didn't see the positive exponent! Curiously, there appears to be just one bad apple: > sort(Mod(p(z))) [1] 1.062855e-10 1.062855e-10 1.328999e-10 1.328999e-10 2.579625e-10 [6] 2.579625e-10 3.834721e-10 3.834721e-10 3.875288e-10 3.875288e-10 [11] 5.287459e-10 5.287459e-10 5.306241e-10 5.306241e-10 6.678424e-10 [16] 6.678424e-10 6.876300e-10 6.876300e-10 8.519089e-10 8.519089e-10 [21] 9.345531e-10 9.345531e-10 9.369015e-10 9.369015e-10 9.789510e-10 [26] 9.789510e-10 1.041601e-09 1.041601e-09 1.046996e-09 1.046996e-09 [31] 1.149073e-09 1.149073e-09 1.209875e-09 1.209875e-09 1.232793e-09 [36] 1.232793e-09 1.244167e-09 1.244167e-09 1.388979e-09 1.388979e-09 [41] 1.556354e-09 1.556354e-09 1.596424e-09 1.596424e-09 1.610661e-09 [46] 1.610661e-09 1.722577e-09 1.722577e-09 1.727842e-09 1.727842e-09 [51] 1.728728e-09 1.728728e-09 1.769644e-09 1.769644e-09 1.863983e-09 [56] 1.863983e-09 1.895296e-09 1.895296e-09 1.907509e-09 1.907509e-09 [61] 1.948662e-09 1.948662e-09 2.027021e-09 2.027021e-09 2.061063e-09 [66] 2.061063e-09 2.116735e-09 2.116735e-09 2.185764e-09 2.185764e-09 [71] 2.209158e-09 2.209158e-09 2.398479e-09 2.398479e-09 2.404217e-09 [76] 2.404217e-09 2.503279e-09 2.503279e-09 2.643436e-09 2.643436e-09 [81] 2.654788e-09 2.654788e-09 2.695735e-09 2.695735e-09 2.921933e-09 [86] 2.921933e-09 2.948185e-09 2.948185e-09 2.953596e-09 2.953596e-09 [91] 3.097433e-09 3.097433e-09 3.420593e-09 3.420593e-09 3.735880e-09 [96] 3.735880e-09 4.190042e-09 4.554964e-09 4.554964e-09 1.548112e+15 > Bill Venables CSIRO Laboratories PO Box 120, Cleveland, 4163 AUSTRALIA Office Phone (email preferred): +61 7 3826 7251 Fax (if absolutely necessary): +61 7 3826 7304 Mobile: +61 4 8819 4402 Home Phone: +61 7 3286 7700 mailto:[EMAIL PROTECTED] http://www.cmis.csiro.au/bill.venables/ -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Venables, Bill (CMIS, Cleveland) Sent: Sunday, 25 May 2008 10:58 AM To: [EMAIL PROTECTED]; [EMAIL PROTECTED] Cc: [EMAIL PROTECTED]; [EMAIL PROTECTED]; [EMAIL PROTECTED] Subject: [ExternalEmail] Re: [R] Solving 100th order equation > library(PolynomF) > x <- polynom() > p <- x^100 - 2*x^99 + 10*x^50 + 6*x - 4000 > z <- solve(p) > z [1] -1.0741267+0.000i -1.073-0.0680356i -1.073+0.0680356i -1.0655699-0.1354644i [5] -1.0655699+0.1354644i -1.0568677-0.2030274i -1.0568677+0.2030274i -1.0400346-0.2687815i ... [93] 1.0595174+0.2439885i 1.0746575-0.1721335i 1.0746575+0.1721335i 1.0828132-0.1065591i [97] 1.0828132+0.1065591i 1.0879363-0.0330308i 1.0879363+0.0330308i 2.000+0.000i > Now to check how good they are: > > range(Mod(p(z))) [1] 1.062855e-10 1.548112e+15 > Not brilliant, but not too bad. Bill Venables CSIRO Laboratories PO Box 120, Cleveland, 4163 AUSTRALIA Office Phone (email preferred): +61 7 3826 7251 Fax (if absolutely necessary): +61 7 3826 7304 Mobile: +61 4 8819 4402 Home Phone: +61 7 3286 7700 mailto:[EMAIL PROTECTED] http://www.cmis.csiro.au/bill.venables/ -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Gabor Grothendieck Sent: Sunday, 25 May 2008 10:09 AM To: Shubha Vishwanath Karanth Cc: [EMAIL PROTECTED]; Duncan Murdoch; Peter Dalgaard Subject: Re: [R] Solving 100th order equation Actually maybe I was premature. It does not handle the polynomial I tried it on in the example earlier in this thread but it does seem to work with the following very simple polynomials of order 100. At any rate it would not take long to try it on the real problem and see. > Solve(x^100 - 1, x) [1] "Starting Yacas!" expression(list(x == complex_cartesian(cos(pi/50), sin(pi/50)), x == complex_cartesian(cos(pi/25), sin(pi/25)), x == complex_cartesian(cos(3 * pi/50), sin(3 * pi/50)), x == complex_cartesian(cos(2 * pi/25), sin(2 * pi/25)), x == complex_cartesian(cos(pi/10), (root(5, 2) - 1)/4), x == complex_cartesian(cos(3 * pi/25), sin(3 * pi/25)), x == complex_cartesian(cos(7 * pi/50), sin(7 * pi/50)), x == complex_cartesian(cos(4 * pi/25), sin(4 * pi/25)), x == complex_cartesian(cos(9 * pi/50), sin(9 * pi/50)), x == complex_cartesian(cos(pi/5), sin(pi/5)), x == complex_cartesian(cos(11 * pi/50), sin(11 * pi/50)), x == complex_cartesian(cos(6 * pi/25), sin(6 * pi/25)), x == complex_cartesian(cos(13 * pi/50), sin(13 * pi/50)), x == complex_cartesian(cos(7 * pi/25), sin(7 * pi/25)), x == complex_cartesian(cos(3 * pi/10), sin(3 * pi/10)), x == complex_cartesian(cos(8 * pi/25), sin(8 * pi/25)), x == complex_cartesian(cos(17 * pi/50), sin(17 * pi/50)), x == complex_cartesian(cos(9 * pi/25), sin(9 * pi/25)), x == complex_cartesian(cos(19 * pi/50), sin(19 * p
Re: [R] Solving 100th order equation
On 25/05/2008, at 1:10 PM, <[EMAIL PROTECTED]>
<[EMAIL PROTECTED]> wrote:
Oops! Actually spectacularly bad. I didn't see the positive
exponent!
Curiously, there appears to be just one bad apple:
sort(Mod(p(z)))
[1] 1.062855e-10 1.062855e-10 1.328999e-10 1.328999e-10 2.579625e-10
[6] 2.579625e-10 3.834721e-10 3.834721e-10 3.875288e-10 3.875288e-10
[11] 5.287459e-10 5.287459e-10 5.306241e-10 5.306241e-10 6.678424e-10
[16] 6.678424e-10 6.876300e-10 6.876300e-10 8.519089e-10 8.519089e-10
[21] 9.345531e-10 9.345531e-10 9.369015e-10 9.369015e-10 9.789510e-10
[26] 9.789510e-10 1.041601e-09 1.041601e-09 1.046996e-09 1.046996e-09
[31] 1.149073e-09 1.149073e-09 1.209875e-09 1.209875e-09 1.232793e-09
[36] 1.232793e-09 1.244167e-09 1.244167e-09 1.388979e-09 1.388979e-09
[41] 1.556354e-09 1.556354e-09 1.596424e-09 1.596424e-09 1.610661e-09
[46] 1.610661e-09 1.722577e-09 1.722577e-09 1.727842e-09 1.727842e-09
[51] 1.728728e-09 1.728728e-09 1.769644e-09 1.769644e-09 1.863983e-09
[56] 1.863983e-09 1.895296e-09 1.895296e-09 1.907509e-09 1.907509e-09
[61] 1.948662e-09 1.948662e-09 2.027021e-09 2.027021e-09 2.061063e-09
[66] 2.061063e-09 2.116735e-09 2.116735e-09 2.185764e-09 2.185764e-09
[71] 2.209158e-09 2.209158e-09 2.398479e-09 2.398479e-09 2.404217e-09
[76] 2.404217e-09 2.503279e-09 2.503279e-09 2.643436e-09 2.643436e-09
[81] 2.654788e-09 2.654788e-09 2.695735e-09 2.695735e-09 2.921933e-09
[86] 2.921933e-09 2.948185e-09 2.948185e-09 2.953596e-09 2.953596e-09
[91] 3.097433e-09 3.097433e-09 3.420593e-09 3.420593e-09 3.735880e-09
[96] 3.735880e-09 4.190042e-09 4.554964e-09 4.554964e-09 1.548112e+15
library(PolynomF)
x <- polynom()
p <- x^100 - 2*x^99 + 10*x^50 + 6*x - 4000
z <- solve(p)
z
[1] -1.0741267+0.000i -1.073-0.0680356i -1.073
+0.0680356i -1.0655699-0.1354644i
[5] -1.0655699+0.1354644i -1.0568677-0.2030274i -1.0568677
+0.2030274i -1.0400346-0.2687815i
...
[93] 1.0595174+0.2439885i 1.0746575-0.1721335i 1.0746575
+0.1721335i 1.0828132-0.1065591i
[97] 1.0828132+0.1065591i 1.0879363-0.0330308i 1.0879363
+0.0330308i 2.000+0.000i
Now to check how good they are:
range(Mod(p(z)))
[1] 1.062855e-10 1.548112e+15
It is fascinating to muck about with curve(p,from=a,t=b) for various
values
of a and b.
One eventually discerns that there is a real root just a tad less
than 2.
Doing uniroot(p,c(1.5,2.5)) gives a root of just *over* 2 with a
function
value of about 1.9e25.
However uniroot(p,c(1.995,2.005)) gives
$root
[1] 1.93
$f.root
[1] -4.570875e+24
$iter
[1] 4
$estim.prec
[1] 6.103516e-05
What a difference 7.214144e-06 makes! When you're dealing with
polynomials of degree 100.
Bozhemoi!
cheers,
Rolf Turner
P.S. I suspect that this was a numerical analysis homework question
designed to teach
a salutary lesson to the nonchalant neophytes.
R. T.
##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
__
R-help@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.
Re: [R] Solving 100th order equation
Rolf Turner auckland.ac.nz> writes: > > [...] > > However uniroot(p,c(1.995,2.005)) gives > > $root > [1] 1.93 > > $f.root > [1] -4.570875e+24 > > $iter > [1] 4 > > $estim.prec > [1] 6.103516e-05 > > What a difference 7.214144e-06 makes! When you're dealing with > polynomials of degree 100. > I don't think R is the right tool to solve this kind of questions that belong to the realm of Computer Algebra Systems. Yacas is much to weak for such high-order polybomials, but we can apply more powerful CASystems, for instance the free Maxima system. Applying (%i)nroots(x^100 - 2*x^99 + 10*x^50 + 6*x - 4000, minf, inf) immediately shows that there are only *two* real solutions, and then (%i)a: realroots(x^100 - 2*x^99 + 10*x^50 + 6*x - 4000, 1e-15); (%i)float(a) (%o)[x=-1.074126672042147,x=1.982] will provide real solutions with 15 decimals (does not change when more decimals are used). So the difference that counts is actually much smaller. The complex roots -- if needed -- will require special treatment. I believe it would be fair to point such questions to Computer Algebra mailing lists and not try to give the appearance they could be solved satisfyingly with R. The same as a complicated statistics question in a Matlab or Mathematica mailing list should be pointed to R-help. Hans Werner > Bozhemoi! > > cheers, > > Rolf Turner __ R-help@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.
Re: [R] Solving 100th order equation
On 26/05/2008, at 10:38 AM, Hans W. Borchers wrote:
Rolf Turner auckland.ac.nz> writes:
[...]
However uniroot(p,c(1.995,2.005)) gives
$root
[1] 1.93
$f.root
[1] -4.570875e+24
Whoops! I read that e+24 as e-24, so scrub all that I said.
(You'd think, that having seen Bill Venables make a similar
error --- which he corrected in the follow-up posting to which
I was replying --- I would've been more careful. Well, you'd
think wrong. Actually it *was* e-24 when I posted; then the
Gremlins got in and changed everything. :-) )
$iter
[1] 4
$estim.prec
[1] 6.103516e-05
What a difference 7.214144e-06 makes! When you're dealing with
polynomials of degree 100.
I don't think R is the right tool to solve this kind of questions
that belong to
the realm of Computer Algebra Systems. Yacas is much to weak for
such high-order
polybomials, but we can apply more powerful CASystems, for instance
the free
Maxima system. Applying
(%i)nroots(x^100 - 2*x^99 + 10*x^50 + 6*x - 4000, minf, inf)
immediately shows that there are only *two* real solutions, and then
(%i)a: realroots(x^100 - 2*x^99 + 10*x^50 + 6*x - 4000, 1e-15);
(%i)float(a)
(%o)[x=-1.074126672042147,x=1.982]
will provide real solutions with 15 decimals (does not change when
more decimals
are used). So the difference that counts is actually much smaller.
The complex roots -- if needed -- will require special treatment.
I believe it would be fair to point such questions to Computer
Algebra mailing
lists and not try to give the appearance they could be solved
satisfyingly with
R. The same as a complicated statistics question in a Matlab or
Mathematica
mailing list should be pointed to R-help.
Be that as it were, it doesn't look like Maxima is doing so wonderfully
either. Consider:
> library(PolynomF)
> x <- polynom()
> p <- x^100-2*x^99+10*x^50+6*x-4000
> a <- 1.982
> p(a)
[1] -1.407375e+14
And notice the + sign in the exponent. :-)
I tried doing
> uniroot(p,c(1.9995,2.0005),tol=.Machine$double.eps)
and got
$root
[1] 2
$f.root
[1] 912
$iter
[1] 5
$estim.prec
[1] 8.881784e-16
Well, 912 is a lot closer to 0 than -4.570875e+24, or even -1.407375e
+14.
But it still ain't 0. Also that ``2'' of course isn't really 2.
I set r <- uniroot(p,c(1.9995,2.0005),tol=.Machine$double.eps)$root
I get p(r) = 912.
But print(r,digits=22) gives 1.982.
And print(a,digits=22) gives the identical result.
Although p(r) and p(a) are wildly different.
I guess it's possible that p(a) --- where ``a'' is as it appears,
written out to 14 decimal places --- really is quite close to zero,
and that -1.407375e+14 is due to round-off error in the calculation
process. But it's also possible that p(a) really is pretty close to
-1.407375e+14, i.e. Maxima isn't really finding the correct root
either.
Or something in between.
I guess that to get a really ``correct'' answer, and check it properly,
you need a system that will do ``arbitrary precision arithmetic''.
Which
isn't R. I don't know if Maxima was using arbitrary precision
arithmetic
in the calculations shown.
Anyhow, I think the real point is that solving 100-th degree
polynomials
is a pretty silly thing to do, whatever package or system or
language you use.
That is, R may be the wrong tool for the job, but it's the wrong job.
And if you *must* do it, then you should do some careful checking
and investigation
of the results --- again irrespective of whatever package or system
or language
you use.
cheers,
Rolf Turner
##
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
__
R-help@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.
