Function for polynomials of arbitrary order (preferred symbolic method)

Question

I 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 the derivative:

 g1 <- Deriv(f1) 

How to create a universal function so that it does not correspond for each new degree of the polynomial NN ?

Answer 1

My initial answer may not be what you really want, as it was pretty conditional. Here is a 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 `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 it easier to obtain higher order derivatives. We can just install nderiv . However, for D we must use recursion (see Examples ?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 an expression, we can evaluate the result later. For instance,

 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 

It follows from the above that we can complete the expression for the function. Using a function has several advantages, one of which is that we can build it using curve or plot.function .

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

Please note: fun requires an expr expression in terms of x . If expr was defined in terms of y , for example, we need to define fun with function (y, expr) . Now use curve to plot cubic and dcubic on the 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 one fun function instead of the f + fun combination. Thus, we also do not need to worry about the consistency of 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 the plot is also simple:

 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 2

Since my final solution with symbolic derivatives is ultimately too large, I use a separate session for numerical calculations. We can do this the same way as for polynomials, the derivatives are clearly known, so we can encode them. Please note: there will be no use of the expression R; everything is done directly using functions.

So, first we create a polynomial basis from degree 0 to degree p - n , then we multiply the coefficient and factor factor. It’s 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 x and pc example 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 a symbolic solution.

Similarly, we can build 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) 

Answer 3

Now, having spent quite a bit of effort to demonstrate how we can resolve this issue ourselves, consider using the polynom R package. As a small package, it is designed to implement construction, derivatives, integration, arithmetic, and the search for the roots of one-dimensional polynomials. This package is completely written with R-language without any compiled code.

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

We are still considering the example of the cubic polynomial used earlier.

 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 that the step-by-step construction of the resulting function has all the parameters inside. For x , only one argument is required. In contrast, functions in the other two answers will also take coefficients and derived order as required arguments. 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 get the same graph that we see in the other two answers, we get other derivatives:

 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> 

You may have already noticed that I did not specify nderiv , but recursively took 1 derivative at a time. This may be a flaw in this package. This does not favor higher order derivatives.

Now we can make the 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)