RON70 <ron_michael70 <at> yahoo.com> writes:
>
>
> 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 means 'PolynomF' which is an improved version of the old polynomial
package.
You do not need to recreate the polynomial from points. Instead, calculate
the exact polynomial:
library(PolynomF)
z <- polynom()
p11 <- 1 - A1[1,1]*z - A2[1,1]*z^2 - A3[1,1]*z^3 - A4[1,1]*z^4
# ...
p <- p11*p22 - p12*p21
There is probably a shorter way to generate these four polynoms p11, ..., p22.
Anyway, the result is
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
solve(p)
# [1] -1.365976+0.000000i -0.737852-1.639581i -0.737852+1.639581i
# [4] -0.012071-1.287727i -0.012071+1.287727i 1.000000+0.000000i
# [7] 1.388794-0.281841i 1.388794+0.281841i
and the real solutions are 1.0 and -1.365976 !
Regards, Hans Werner
> 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: rvaradhan <at> jhmi.edu
> >
______________________________________________
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.
