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.

Reply via email to