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.