Hi,
You can use the "polynomial" package to solve your problem.
The key step is to find the exact polynomial representation of fn(). Noting
that it is a 8-th degree polynomial, we can get its exact form using the
poly.calc() function. Once we have that, it is a simple matter of finding the
roots using the solve() function.
require(polynomial)
a <- c(-0.07, 0.17)
b <- c(1, -4)
cc <- matrix(c(0.24, 0.00, -0.08, -0.31), 2)
d <- matrix(c(0, 0, -0.13, -0.37), 2)
e <- matrix(c(0.2, 0, -0.06, -0.34), 2)
A1 <- diag(2) + a %*% t(b) + cc
A2 <- -cc + d
A3 <- -d + e
A4 <- -e
# I am vectorizing your function
fn <- function(z)
{
sapply(z, function(z) {
y <- diag(2) - A1*z - A2*z^2 - A3*z^3 - A4*z^4
det(y)
})
}
x <- seq(-5, 5, length=9) # note we need 9 points to exactly determine a 8-th
degree polynomial
y <- fn(x)
p <- poly.calc(x, y) # uses Lagrange interpolation to determine polynomial form
p
> 1 - 1.18*x + 0.2777*x^2 - 0.2941*x^3 - 0.1004*x^4 + 0.3664*x^5 - 0.0636*x^6 +
> 0.062*x^7 - 0.068*x^8
# plot showing that p is the exact polynomial representation of fn(z)
pfunc <- as.function(p)
x1 <-seq(-5, 5, length=100)
plot(x1, fn(x1),type="l")
lines(x1, pfunc(x1), col=2, lty=2)
solve(p) # gives you the roots (some are, of course, complex)
Hope this helps,
Ravi.
____________________________________________________________________
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
----- Original Message -----
From: RON70 <ron_michae...@yahoo.com>
Date: Sunday, January 11, 2009 1:05 am
Subject: [R] How to get solution of following polynomial?
To: r-help@r-project.org
> Hi, I want find all roots for the following polynomial :
>
> a <- c(-0.07, 0.17); b <- c(1, -4); cc <- matrix(c(0.24, 0.00, -0.08,
> -0.31), 2); d <- matrix(c(0, 0, -0.13, -0.37), 2); e <- matrix(c(0.2,
> 0,
> -0.06, -0.34), 2)
> A1 <- diag(2) + a %*% t(b) + cc; A2 <- -cc + d; A3 <- -d + e; A4 <- -e
> fn <- function(z)
> {
> y <- diag(2) - A1*z - A2*z^2 - A3*z^3 - A4*z^4
> return(det(y))
> }; uniroot(fn, c(-10, 1))
>
> Using uniroot function, I got only one solution of that. Is there any
> function to get all four solutions? I looked at polyroot() function,
> but I
> do not think it will work for my problem, because, my coef. are
> matrix, nor
> number
>
> Thanks
>
>
>
> --
> View this message in context:
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help@r-project.org mailing list
>
> PLEASE do read the posting guide
> and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help@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.
