[R] complex conjugates roots from polyroot?

[Ravi Varadhan]

Hi Spencer,

Here is a simple approach to detect conjugate pairs:

is.conj <- function(z1, z2, tol=1.e-10) {
# determine if two complex numbers are conjugates
cond1 <- abs(Re(z1) - Re(z2)) < tol
cond2 <- abs(Im(z1) + Im(z2)) < tol
cond1 & cond2
}

set.seed(123)
z <- polyroot(sample(1:5, size=8, rep=T))
zmat <- which(outer(z, z, FUN="is.conj"), arr.ind=T)
zmat[zmat[,1] < zmat[,2], ]

# result
     row col
[1,]   1   3
[2,]   5   6
[3,]   4   7
>>

We see that (1,3), (4,7), and (5,6) are the conjugate pairs.

This doesn't address the issue of numerical round-off (there is no argument
in polyroot that governs the accuracy of the roots).

Best,
Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 

----------------------------------------------------------------------------
--------

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Spencer Graves
Sent: Friday, November 23, 2007 12:08 PM
To: r-help at r-project.org
Subject: [R] complex conjugates roots from polyroot?

Hi, All: 

      Is there a simple way to detect complex conjugates in the roots 
returned by 'polyroot'?  The obvious comparison of each root with the 
complex conjugate of the next sometimes produces roundoff error, and I 
don't know how to bound its magnitude: 

(tst <- polyroot(c(1, -.6, .4)))
tst[-1]-Conj(tst[-2])
[1] 3.108624e-15+2.22045e-16i
abs(tst[-1]-Conj(tst[-2]))/abs(tst[-1])
1.971076e-15
.Machine$double.neg.eps
1.110223e-16

      Testing (abs(tst[-1]-Conj(tst[-2]))/abs(tst[-1]) < 
(20*.Machine$double.neg.eps)) would catch this example, but it might not 
catch others. 
   
      The 'polyroot' help page says, "See Also ... the 'zero' example in 
the demos directory."  Unfortunately, I've so far been unable to find 
that. 

      Any suggestions? 
      Thanks in advance. 
      Spencer Graves

______________________________________________
R-help at 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.