Solving 100th order equation
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}}