Home > Software engineering >  How to write a function for polynomial functions in R to make them plottable
How to write a function for polynomial functions in R to make them plottable

Time:01-01

I want to plot polynomial functions of high degrees with given coefficiants. I created this function 'f_erzeuger' in order to write a polynomial function 'f' to be able to plot it with ggplot2 and stat_function(fun = f). 'coef' is a vector of coefficiants for the ploynomial function.

f_erzeuger <- function(coef) {
  f <- function(x) {
    x_vec <- c(1, x, x^2, x^3, x^4, x^5, x^6, x^7, x^8, x^9, x^10)
    d <- length(coef) - 1 # degree of the polynomial
    sum(coef * x_vec[1:(d 1)])
  }
  f
} 

However ggplot2 is not able to plot this curve, propably because the function doesn't actually calculate a function term.

f <- f_erzeuger(c(4, 5, 2))
ggplot(data.frame(x = c(-50, 50)), aes(x = x))  
  stat_function(fun = f, color = "blue")  
  stat_function(fun = function(x){3   5*x   2*x^2}, color = "red")

The created function f is beeing plottet as a line even though f works as a function.

ggplot

Does anyone have an idea how to fix this? Thank you in advance!

CodePudding user response:

The problem is that your function is not vectorized. Running ?stat_function and looking at the documentation for fun:

fun: Function to use. Either 1) an anonymous function in the base or rlang formula syntax (see rlang::as_function()) or 2) a quoted or character name referencing a function; see examples. Must be vectorised.

To make the function vectorized, we need to make sure that, for example, f(c(0, 1)) will return c(f(0), f(1)). Note that one issue in your function is the line where you define x_vec = c(1, x, ...). This wouldn't work if x were a vector with more than one element.

There are many ways to vectorize your function. I will do it using the tidyverse (mainly purrr::map()).

f_erzeuger = function(coef) {
  function(xvals) {
    d = length(coef)
    map_dbl(xvals, function(x) {
      x_vec = x ^ (0:(d - 1))
      sum(coef * x_vec)
    })
  }
}

The changes made to this function:

  • Most importantly, the function is now vectorized.
  • Instead of defining x_vec explicitly to degree 10, we can leverage the fact that ^ in R is vectorized (so x^(0:2) is the same as c(1, x, x^2)).
  • We can return the function directly inside f_erzeuger() instead of defining f and then returning it.

Now things will work as expected: enter image description here

CodePudding user response:

Try curve.

f <- f_erzeuger(c(4, 5, 2))
curve(Vectorize(f)(x), -50, 50, col='#c55d0d')

enter image description here

  • Related