On 26/05/2008, at 10:38 AM, Hans W. Borchers wrote:
Rolf Turner <r.turner <at> auckland.ac.nz> writes:
[...]
However uniroot(p,c(1.995,2.005)) gives
$root
[1] 1.999993
$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.999999999999982]
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.999999999999982
> 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.999999999999982.
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.
