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

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.

CodePudding user response：

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)

etc. (and I might have made some mistakes above, you'd better check it).

Then use

polyroot(c(a0, a1, a2, a3))

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.

CodePudding user response：

You can create a table with a column for each variable and filter the rows not satisfying your equation:

library(tidyverse)

set.seed(1337)
n <- 1000

tibble(
  a = runif(n),
  b = runif(n)
) |>
  filter(a   b < 1 - a   b) |>
  expand_grid(
    theta = seq(0, 1, by = 1),
    x = seq(0, 1, by = 1)
  ) |>
  filter(
    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
  )
#> # A tibble: 0 × 4
#> # … with 4 variables: a <dbl>, b <dbl>, theta <dbl>, x <dbl>

^{Created on 2022-06-13 by the reprex package (v2.0.0)}

Unfortunately, there is no point in the sampled space satisfying your equation. This is probably due to ==0 instead of <e where e is a very small error. One needs to allow small errors in numerical sampling solutions.

Why just not solve the roots of the equation analytically?