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