[R] Solving 100th order equation

[Rolf Turner]

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}}