6

I want to find the solution of: 

-x^3+6*x^2+51*x+44=0

but with R. Is it possible?

I found the package Ryacas, but nobody seems to be able to make it work.

May sound trivial, but I'm not able to find an easy way to do this...

Do you have an alternative?

Thanks guys!

M. Beausoleil
  • 3,141
  • 6
  • 29
  • 61
  • What do you mean by "nobody seems to be able to make it work"? –  Sep 16 '15 at 03:26
  • 1
    Maybe it's possible, but R is designed for stats, not math. Try http://www.wolframalpha.com/input/?i=solve+-x%5E3%2B6*x%5E2%2B51*x%2B44%3D0 – Frank Sep 16 '15 at 03:30
  • yes, you use `uniroot`, `f <- function(x) -x^3+6*x^2+51*x+44; uniroot(f, interval=c(-100, 100))` – Rorschach Sep 16 '15 at 03:31
  • 1
    @bunk You only get 1 solution. There are 3. See `curve(f, -10, 15); abline(h=0)`. –  Sep 16 '15 at 03:37
  • @Pascal true, but it's not too hard to adapt the method to the rest (just move the interval). There are other root finding functions as well, I can't remember them though. edit: apparently there is a uniroot.all, whcih sounds promising – Rorschach Sep 16 '15 at 03:40
  • @Pascal I mean, that I did a little research on the internet and it's always the same story with: sh: yacas: command not found Error in socketConnection(host = "127.0.0.1", port = 9734, server = FALSE, : cannot open the connection In addition: Warning message: In socketConnection(host = "127.0.0.1", port = 9734, server = FALSE, : 127.0.0.1:9734 cannot be opened – M. Beausoleil Sep 16 '15 at 03:51
  • 1
    Yes, "yacas" is an open-source software you need to install on your machine. –  Sep 16 '15 at 03:53
  • Do you think it worth the installation? – M. Beausoleil Sep 16 '15 at 03:59

3 Answers3

11

You can use polynom package:

library(polynom)
p <- polynomial(c(44,51,6,-1))
# 44 + 51*x + 6*x^2 - x^3 
solve(p)
# [1] -4 -1 11

But you simply can use the function polyroot from base package:

polyroot(c(44,51,6,-1))
# [1] -1+0i -4+0i 11+0i

If you keep the real part with Re:

Re(polyroot(c(44,51,6,-1)))
# [1] -1 -4 11
6

Here we solve for the roots using the relationship between a matrix and its characteristic polynomial.

Given the polynomial a0 + a1*x^1 + a2*x^2 + x^3, define the matrix:

0   0  -a0
1   0  -a1
0   1  -a2

The eigenvalues of this matrix are the roots of the polynomial.

Substituting y = -x in your polynomial equation gives this 

y^3 + 6*y^2 - 51*y + 44=0

And gives this example

> z <- matrix(c(0,1,0,0,0,1,-44,51,-6),3,3)
> z
     [,1] [,2] [,3]
[1,]    0    0  -44
[2,]    1    0   51
[3,]    0    1   -6
> eigen(z)
$values
[1] -11   4   1

$vectors
           [,1]        [,2]        [,3]
[1,]  0.6172134  0.73827166  0.98733164
[2,] -0.7715167 -0.67115606 -0.15707549
[3,]  0.1543033 -0.06711561 -0.02243936

Or, since we've substituted -y for x:

> eigen(-z)$values
[1] 11 -4 -1

See: http://www-math.mit.edu/~edelman/publications/polynomial_roots.pdf

Matthew Lundberg
  • 42,009
  • 6
  • 90
  • 112
4

I just stumbled upon this question and I am not sure if anything inherently changed around the Ryacas package, but it seems to work great in 2020, here is a helpful vignette to get started: https://cran.r-project.org/web/packages/Ryacas/vignettes/getting-started.html

Following the vignette, things seem to work as expected when I run the code:

library(Ryacas)
# initialize equation:
eq <- "-x^3+6*x^2+51*x+44"
# simplify the equation:
library(glue)
yac_str(glue("Simplify({eq})"))

[1] "6*x^2-x^3+51*x+44"

# factor:
yac_str(glue("Factor({eq})"))

[1] "(-1)*(x-11)*(x+4)*(x+1)"

You can evaluate the expression like this, plugging in whatever values for x:

# evaluate
evaluate(eq,list(x=c(0,1,10,100,-100)))

[[1]]
$src
[1] "-x^3+6*x^2+51*x+44"

attr(,"class")
[1] "source"

[[2]]
[1] "[1]      44     100     154 -934856 1054944\n"

Here you can see the results where x=0 produced an answer of 44, x=1 produced an answer of 100, etc...

If you evaluated the new simplified or factored versions and evaluated those, you would of course end up with the same exact results:

evaluate(yac_str(glue("Simplify({eq})")),list(x=c(0,1,10,100,-100)))

[[1]]
$src
[1] "6*x^2-x^3+51*x+44"

attr(,"class")
[1] "source"

[[2]]
[1] "[1]      44     100     154 -934856 1054944\n"

Notice the formula changed in the $src output, but we get the same results.

Here's the factored one too: evaluate(yac_str(glue("Factor({eq})")),list(x=c(0,1,10,100,-100)))

[[1]]
$src
[1] "(-1)*(x-11)*(x+4)*(x+1)"

attr(,"class")
[1] "source"

[[2]]
[1] "[1]      44     100     154 -934856 1054944\n"

The only real difference between what I outlined here and what's outlined in the vignette is the actual formula, and the fact that I used library(glue) instead of paste0, which is also a fair option.

Ricky
  • 1,005
  • 1
  • 12
  • 15
  • 1
    You can also just use `yac("Solve(-x^3+6*x^2+51*x+44==0, x)")` and it gives the response `"{x==(-1),x==(-4),x==11}"` – Vlad Aug 15 '20 at 04:30