Hi Ravi, Thanks for this reply. However I could not understand meaning of "vectorizing the function". Can you please be little bit elaborate on that? Secondly the package "polynomial" is not available in CRAN it seems. What is the alternate package?
Thanks, Ravi Varadhan wrote: > > 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. > > -- View this message in context: http://www.nabble.com/How-to-get-solution-of-following-polynomial--tp21396441p21406350.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ 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.