Search code examples
rpolynomialsrgl

Creating a 3D Plot of a Polynomial Function with Uniform Distributed Values

I have an equation which goes like this,

2* (1-x-a-b)^2 * x * *theta* + 2 * (1-a-b-x) * x^2 * *theta*  - 2 * b * x^2 + 2 * a * (1-a-b-x)^2 = 0

I want to create a function in R, that selects a and b with restriction (a + b < 1 - a + b) from an uniform distribution. After selecting, I want it to find the solutions for x (both negative and positive).

I want to repeat this process t amount of time in a for loop where I will give the theta value as an input.

After that I want it to create a 3D density plot where solutions are shown with respect to values of a,b on two axes and x on one axis.

So far I have tried to use polynom package and solve function. But I am having hard time with R when it comes to mathematics.

Solution

  • You need to rewrite the polynomial in standard form a0 + a1*x + a2*x^2 + a3*x^3, then you can use the base function polyroot() to find the roots. For example,

    a0 <- 2 * a * (1 - a - b)^2
a1 <- 2 * (1 - a - b)^2 * theta - 4 * a * (1 - a - b)
a2 <- -4 * (1 - a - b) * theta + 2 * (1 - a - b) * theta - 2  * b + 2 * a
a3 <- 0

    So this is a quadratic equation, not a cubic as it appears at first glance.

    Then use

    polyroot(c(a0, a1, a2))

    to find the roots. Select the real roots, and put them together into a matrix roots with columns a, b, root, then use rgl::plot3d(roots) to display them.

    I think you have a typo in your restriction, so I'll ignore it, and this is the plot I get for theta == 1:

    theta <- 1
a <- runif(1000)
b <- runif(1000)

a0 <- 2*a*(1-a-b)^2
a1 <- 2*(1-a-b)^2*theta -4*a*(1-a-b)
a2 <- -4*(1-a-b)*theta + 2*(1-a-b)*theta-2*b+2*a

result <- matrix(numeric(), ncol = 3, dimnames = list(NULL, c("a", "b", "root")))
for (i in seq_along(a)) {
  root <- polyroot(c(a0[i], a1[i], a2[i]))
  if (max(Im(root)) < 1.e8)
    result <- rbind(result, cbind(a[i], b[i], Re(root)))
}
library(rgl)
plot3d(result)

    Created on 2022-06-14 by the reprex package (v2.0.1)

    Most of the roots are really small, but for some of them a2 is nearly zero, and then they can be very large.