Differentiation and integration in R

N.B. Venkateswarlu, GVPCEW
venkat_ritch@yahoo.com

R base system supports both differentiation and integration. We can use function D() to calculate differentiation of a function. Argument for this function is expression object having the mathematical expression for which we want to find differentiation.

First, we take the equation as an expression like the following.

f=expression(x^2+3*x)

To calculate first derivative of f, we use D() function and ‘x’ to specify that derivation has to be carried out with respect to x (Of course, the above expression, f, contains only x as the independent variable).

D(f,'x')
## 2 * x + 3

We get 2 * x + 3 as output.

The following R command can be used to find second derivative of the above f.

D(D(f,'x'),'x')
## [1] 2

We get 2 as output.

Assuming that, the returned value of function D() can be assigned to an object and send the same as argument to D() to compute second derivative of f.

ff<-D(f,'x')

D(ff,'x')
## [1] 2

Here, also we get 2 as second derivative.

We can also evaluate expression values using eval() function.

x<-1:5

eval(f)

eval(ff)

The following figure 1, illustrates the working of the above concepts.

Figure 1.

Figure 1.

Exercise: Tabulate function f (f=x3-2x2+3) values and its derivative values in the range -3,3 with the increment value of 0.1 using R statements.

Answer: See the following figure 2 which contains our workout.

Figure 2.

Figure 2.

If the expression is having more than one independent variable, we can calculate differentiation with respect to each of them.

f<-expression(x^2+y^2+2*x*y-3*x+4*y+4)
D(f,'x')
## 2 * x + 2 * y - 3
D(f,'y')
## 2 * y + 2 * x + 4

We can also calculate derivatives of a function involving trigonometric functions also. For example, we have defined an expression, exp, with sine and cosine terms.

exp<- expression(sin(cos(x + y^2)))

The following command’s give derivative with respect to x and displays the same.

dx<-D(exp,“x”); dx

The following command’s compute derivative with respect to y and displays the same.

dy<-D(exp,“y”); dy

We can also calculate all the (partial) derivatives of an equation using deriv() function. It returns the result as an expression.

dxy<-deriv(exp,c(“x”,“y”));dxy;typeof(dxy)

In the following, we are defining some set of values for x and y and then computing the partial derivatives with respect to x and y using their representation in expression form in dxy.

x<-seq(-pi,pi,pi/4)

y<-pi

eval(dxy)

We can also use an option func=T with derive function such that it returns all the derivatives of the expression as a function which we can invoke at any time. pp<-deriv(exp,c(“x”,“y”),func=T) pp(x,y)

The following figure 90, illustrates our workout of the above R commands at the R console.

Figure 3

Figure 3

Let us explore how to compute higher order derivatives. We use a function DD() for this. We advise readers to come to this concept once they are exposed to function writing in R. The function, DD, takes expression for which we want to compute derivative as first argument, with which variable we want derivative as a second argument and derivative number (first or second or third, etc) as the third argument. It is implemented as a recursive function as shown below.

DD<-function(expr,name,order=1){
if(order<1) stop("Order must be >=1")
if(order==1)D(expr,name)
else DD(D(expr,name),name,order-1)
}

We have defined an expression.

exp<-expression(x^3)

To compute first derivative of exp, run the following.

DD(exp,"x",1)
## 3 * x^2

To compute second derivative of exp, run the following.

DD(exp,"x",2)
## 3 * (2 * x)

To compute third derivative of exp, run the following.

DD(exp,"x",3)
## 3 * 2

The following is our workout(Figure 4) at the R console. We can confirm that the first derivative of exp (x3) is 3x2, 6x and 6 which are coming from the DD function.

Figure 2

Figure 2

To calculate, integration value of a function, we first define a function (with name f or some other name0 for the function as shown below.

f<-function(x) x^2+3*x

Note: As we have not introduced about R functions, we request readers to comeback to this sections once after going through the chapter on functions.

We use integrate() function and supply lower and upper limits along the with function f.

As we want integration between 0 and 1, we have given the following command.

integrate(f,0,1)
## 1.833333 with absolute error < 2e-14

We get the following output: 1.833333 with absolute error < 2e-14

We want to integrate between 0 to infinity, we have used the following.

integrate(f,0,Inf)

We get the following output: Error in integrate(f, 0, Inf) : the integral is probably divergent

We can also use the following type integrate function call.

integrate(f, lower=0, upper=Inf)

The following is the snapshot (Figure 5) of our workout with integrate() function.

Figure 5

Figure 5

If we want to compute double, triple and higher order integrals, we can use adaptIntegrate() function of cubature package.

For instance, assume that we want to calculate volume of a cube as shown in the following figure 6.

Figure 6

Figure 6

We know from calculus fundamentals that the volume of the above cube can be written as a triple integral.

Figure 7

Figure 7

The following is our workout(Figure 8) at R console. We know the area of the cube is 43=64. The function adaptIntegrate() is also giving the same as its output. First, we have defined function f, with the equation that has to be used for triple integration. In this case, equation is simple 1 as shown in the above integral equation. Then, we call adaptIntegrate() function with lower and upper limits as 0,4 in all the three integrations.

Figure 8

Figure 8

Example: How to calculate the following triple integral?

Figure 9

Figure 9

Answer:

library(cubature) # load the package "cubature" 
## Warning: package 'cubature' was built under R version 3.3.3
f <- function(x) { 2/3 * (x[1] + x[2] + x[3]) } 
# "x" is vector x[1], x[2], x[3] are referring to x1, x2 and x3 respectively.
adaptIntegrate(f, lowerLimit = c(0, 0, 0), upperLimit = c(0.5, 0.5, 0.5))
## $integral
## [1] 0.0625
## 
## $error
## [1] 1.387779e-17
## 
## $functionEvaluations
## [1] 33
## 
## $returnCode
## [1] 0
If we want differentiation and integration in discrete data(not for equations), we can do use diff() and cumsum() function like the following. 
> x<-c(10,20,30,20,12,20,20)
> diff(x)
[1] 10 10 -10 -8 8 0
> cumsum(x)
[1] 10 30 60 80 92 112 132