[R] How to get solution of following polynomial?

[Ravi Varadhan]

Hi Hans,

I actaually meant the "polynom" package (not "polynomial", which was a
typo).  I am curious as to the main differences between "polynom" and
"PolynomF".  

Ron - by vectorizing, I mean that the function fn() can take a vector as an
input and return the function values at all the points in the vector.
"sapply" is an easy way to do this.  

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 Hans W. Borchers
Sent: Monday, January 12, 2009 4:30 AM
To: r-help at stat.math.ethz.ch
Subject: Re: [R] How to get solution of following polynomial?

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 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.