[R] Solving 100th order equation

[Rolf Turner]

On 25/05/2008, at 1:10 PM, <Bill.Venables at csiro.au>  
<Bill.Venables at csiro.au> 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

	<snip>

>> library(PolynomF)
>> x <- polynom()
>> p <- x^100 - 2*x^99 + 10*x^50 + 6*x - 4000
>> z <- solve(p)
>> z
>   [1] -1.0741267+0.0000000i -1.0739999-0.0680356i -1.0739999 
> +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.0000000+0.0000000i
>>
>
> 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.999993

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