Function for polynomials of arbitrary order (symbo

Question

I've found polynomial coefficients from my data:

R <- c(0.256,0.512,0.768,1.024,1.28,1.437,1.594,1.72,1.846,1.972,2.098,2.4029)
Ic <- c(1.78,1.71,1.57,1.44,1.25,1.02,0.87,0.68,0.54,0.38,0.26,0.17)
NN <- 3
ft <- lm(Ic ~ poly(R, NN, raw = TRUE))
pc <- coef(ft)

So I can create a polynomial function:

f1 <- function(x) pc[1] + pc[2] * x + pc[3] * x ^ 2 + pc[4] * x ^ 3

And for example, take a derivative:

g1 <- Deriv(f1)

How to create a universal function so that it doesn't have to be rewritten for every new polynomial degree NN?

Answer 1

Now after quite much effort in demonstrating how we can work out this question ourselves, consider using R package polynom. As a small package, it aims at implementing construction, derivatives, integration, arithmetic and roots-finding of univariate polynomials. This package is written completely with R language, without any compiled code.

## install.packages("polynom")
library(polynom)

We still consider the cubic polynomial example used before.

pc <- 1:4 / 10

## step 1: making a "polynomial" object as preparation
pcpoly <- polynomial(pc)
#0.1 + 0.2*x + 0.3*x^2 + 0.4*x^3 

## step 2: compute derivative
expr <- deriv(pcpoly)

## step 3: convert to function
g1 <- as.function(expr)

#function (x) 
#{
#    w <- 0
#    w <- 1.2 + x * w
#    w <- 0.6 + x * w
#    w <- 0.2 + x * w
#    w
#}
#<environment: 0x9f4867c>

Note, by step-by-step construction, the resulting function has all parameters inside. It only requires a single argument for x value. In contrast, functions in the other two answers will take coefficients and derivative order as mandatory arguments, too. We can call this function

g1(seq(0, 1, 0.2))
# [1] 0.200 0.368 0.632 0.992 1.448 2.000

To produce the same graph we see in other two answers, we get other derivatives as well:

g0 <- as.function(pcpoly)  ## original polynomial

## second derivative
expr <- deriv(expr)
g2 <- as.function(expr)
#function (x) 
#{
#    w <- 0
#    w <- 2.4 + x * w
#    w <- 0.6 + x * w
#    w
#}
#<environment: 0x9f07c68>

## third derivative
expr <- deriv(expr)
g3 <- as.function(expr)
#function (x) 
#{
#    w <- 0
#    w <- 2.4 + x * w
#    w
#}
#<environment: 0x9efd740>

Perhaps you have already noticed that I did not specify nderiv, but recursively take 1 derivative at a time. This may be a disadvantage of this package. It does not facilitate higher order derivatives.

Now we can make a plot

## As mentioned, `g0` to `g3` are parameter-free
curve(g0(x), from = 0, to = 5)
curve(g1(x), add = TRUE, col = 2)
curve(g2(x), add = TRUE, col = 3)
curve(g3(x), add = TRUE, col = 4)

Answer 2

My original answer may not be what you really want, as it was numerical rather symbolic. Here is the symbolic solution.

## use `"x"` as variable name
## taking polynomial coefficient vector `pc`
## can return a string, or an expression by further parsing (mandatory for `D`)
f <- function (pc, expr = TRUE) {
  stringexpr <- paste("x", seq_along(pc) - 1, sep = " ^ ")
  stringexpr <- paste(stringexpr, pc, sep = " * ")
  stringexpr <- paste(stringexpr, collapse = " + ")
  if (expr) return(parse(text = stringexpr))
  else return(stringexpr)
  }

## an example cubic polynomial with coefficients 0.1, 0.2, 0.3, 0.4
cubic <- f(pc = 1:4 / 10, TRUE)

## using R base's `D` (requiring expression)
dcubic <- D(cubic, name = "x")
# 0.2 + 2 * x * 0.3 + 3 * x^2 * 0.4

## using `Deriv::Deriv`
library(Deriv)

dcubic <- Deriv(cubic, x = "x", nderiv = 1L)
# expression(0.2 + x * (0.6 + 1.2 * x))

Deriv(f(1:4 / 10, FALSE), x = "x", nderiv = 1L)  ## use string, get string
# [1] "0.2 + x * (0.6 + 1.2 * x)"

Of course, Deriv makes higher order derivatives easier to get. We can simply set nderiv. For D however, we have to use recursion (see examples of ?D).

Deriv(cubic, x = "x", nderiv = 2L)
# expression(0.6 + 2.4 * x)

Deriv(cubic, x = "x", nderiv = 3L)
# expression(2.4)

Deriv(cubic, x = "x", nderiv = 4L)
# expression(0)

If we use expression, we will be able to evaluate the result later. For example,

eval(cubic, envir = list(x = 1:4))  ## cubic polynomial
# [1]  1.0  4.9 14.2 31.3

eval(dcubic, envir = list(x = 1:4))  ## its first derivative
# [1]  2.0  6.2 12.8 21.8

The above implies that we can wrap up an expression for a function. Using a function has several advantages, one being that we are able to plot it using curve or plot.function.

fun <- function(x, expr) eval.parent(expr, n = 0L)

Note, the success of fun requires expr to be an expression in terms of symbol x. If expr was defined in terms of y for example, we need to define fun with function (y, expr). Now let's use curve to plot cubic and dcubic, on a range 0 < x < 5:

curve(fun(x, cubic), from = 0, to = 5)  ## colour "black"
curve(fun(x, dcubic), add = TRUE, col = 2)  ## colour "red"

The most convenient way, is of course to define a single function FUN rather than doing f + fun combination. In this way, we also don't need to worry about the consistency on the variable name used by f and fun.

FUN <- function (x, pc, nderiv = 0L) {
  ## check missing arguments
  if (missing(x) || missing(pc)) stop ("arguments missing with no default!")
  ## expression of polynomial
  stringexpr <- paste("x", seq_along(pc) - 1, sep = " ^ ")
  stringexpr <- paste(stringexpr, pc, sep = " * ")
  stringexpr <- paste(stringexpr, collapse = " + ")
  expr <- parse(text = stringexpr)
  ## taking derivatives
  dexpr <- Deriv::Deriv(expr, x = "x", nderiv = nderiv)
  ## evaluation
  val <- eval.parent(dexpr, n = 0L)
  ## note, if we take to many derivatives so that `dexpr` becomes constant
  ## `val` is free of `x` so it will only be of length 1
  ## we need to repeat this constant to match `length(x)`
  if (length(val) == 1L) val <- rep.int(val, length(x))
  ## now we return
  val
  }

Suppose we want to evaluate a cubic polynomial with coefficients pc <- c(0.1, 0.2, 0.3, 0.4) and its derivatives on x <- seq(0, 1, 0.2), we can simply do:

FUN(x, pc)
# [1] 0.1000 0.1552 0.2536 0.4144 0.6568 1.0000

FUN(x, pc, nderiv = 1L)
# [1] 0.200 0.368 0.632 0.992 1.448 2.000

FUN(x, pc, nderiv = 2L)
# [1] 0.60 1.08 1.56 2.04 2.52 3.00

FUN(x, pc, nderiv = 3L)
# [1] 2.4 2.4 2.4 2.4 2.4 2.4

FUN(x, pc, nderiv = 4L)
# [1] 0 0 0 0 0 0

Now plotting is also easy:

curve(FUN(x, pc), from = 0, to = 5)
curve(FUN(x, pc, 1), from = 0, to = 5, add = TRUE, col = 2)
curve(FUN(x, pc, 2), from = 0, to = 5, add = TRUE, col = 3)
curve(FUN(x, pc, 3), from = 0, to = 5, add = TRUE, col = 4)

Answer 3

Since my final solution with symbolic derivatives eventually goes too long, I use a separate session for numerical calculations. We can do this as for polynomials, derivatives are explicitly known so we can code them. Note, there will be no use of R expression here; everything is done directly by using functions.

So we first generate polynomial basis from degree 0 to degree p - n, then multiply coefficient and factorial multiplier. It is more convenient to use outer than poly here.

## use `outer`
g <- function (x, pc, nderiv = 0L) {
  ## check missing aruments
  if (missing(x) || missing(pc)) stop ("arguments missing with no default!")
  ## polynomial order p
  p <- length(pc) - 1L
  ## number of derivatives
  n <- nderiv
  ## earlier return?
  if (n > p) return(rep.int(0, length(x)))
  ## polynomial basis from degree 0 to degree `(p - n)`
  X <- outer(x, 0:(p - n), FUN = "^")
  ## initial coefficients
  ## the additional `+ 1L` is because R vector starts from index 1 not 0
  beta <- pc[n:p + 1L]
  ## factorial multiplier
  beta <- beta * factorial(n:p) / factorial(0:(p - n))
  ## matrix vector multiplication
  drop(X %*% beta)
  }

We still use the example x and pc defined in the symbolic solution:

x <- seq(0, 1, by = 0.2)
pc <- 1:4 / 10

g(x, pc, 0)
# [1] 0.1000 0.1552 0.2536 0.4144 0.6568 1.0000

g(x, pc, 1)
# [1] 0.200 0.368 0.632 0.992 1.448 2.000

g(x, pc, 2)
# [1] 0.60 1.08 1.56 2.04 2.52 3.00

g(x, pc, 3)
# [1] 2.4 2.4 2.4 2.4 2.4 2.4

g(x, pc, 4)
# [1] 0 0 0 0 0 0

The result is consistent with what we have with FUN in the the symbolic solution.

Similarly, we can plot g using curve:

curve(g(x, pc), from = 0, to = 5)
curve(g(x, pc, 1), from = 0, to = 5, col = 2, add = TRUE)
curve(g(x, pc, 2), from = 0, to = 5, col = 3, add = TRUE)
curve(g(x, pc, 3), from = 0, to = 5, col = 4, add = TRUE)