source("http://macosa.dima.unige.it/r.R")    # If I have not already loaded the library
---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
S 08 Integrals and Derivatives (and Limits). Taylor polynomial
 
# deriv(f,VariableName), deriv2, …, deriv6
 
#  Note:  deriv(f,a)  equals  D(body(f),a)
 
G = function(u) sin(u); deriv(G,"u")
# cos(u)
f = function(x) 3*x^3-5*x^2+x+3
graphF(f, -2,4, "blue")
deriv(f,"x")
#  3 * (3 * x^2) - 5 * (2 * x) + 1
# To use the expression produced by "deriv" you need to evaluate it with eval:
df = function(x) eval(deriv(f,"x"))
graph(df,-3,5, "red")
# I could also define  df  by writing  3*(3*x^2)-5*(2*x)+1   [ or 9*x^2-10*x+1 ]
deriv2(f,"x")
#  3 * (3 * (2 * x)) - 5 * 2
d2f = function(x) eval(deriv2(f,"x")); graph(d2f,-3,5, "violet")
deriv3(f,"x")
#  3 * (3 * 2)
d3f = function(x) eval(deriv3(f,"x"))+x-x      # I put the "x" with a trick
graph(d3f,-3,5, "green")
                        
# To evaluate the derivative of log10(x) and log2(x) you must use log(x)/log(10) and
# log(x)/log(2), where log(x) is the inverse of exp(x)
f=function(x) log(x); g=function(x) log(x)/10; h=function(x) log(x)/2
deriv(f,"x"); deriv(g,"x"); deriv(h,"x")
# 1/x  1/x/10  1/x/2
 
# To evaluate the derivative of a function where abs(u) appears, it needs to be
# replaced with sqrt(u^2)
#
# I can calculate the integral of a piecewise linear function (with coordinates in x
# and y), using the areaPol command, if I define:
INTEGRAL = function(x,y) areaPol( c(x[length(x)],x[1],x), c(0,0,y) )
# Example:
                    
x = c(10,13, 17, 19); y = c(7,1,9,8)
Plane(min(x),max(x), 0,max(y)); polyl(x,y,"blue")
INTEGRAL = function(x,y) areaPol( c(x[length(x)],x[1],x), c(0,0,y) )
INTEGRAL(x,y)
# 49
# But I can use the integral(f, a,b) command, if I define such function explicitly:
 
f = function(x) ifelse(x<13, 7-6/3*(x-10), ifelse(x<17, 1+8/4*(x-13), 9-1/2*(x-17)))
Plane(10,19, 0,9); graph(f, 10,19, "blue")
                    
integral(f, 10,19)
# 49
 
# A function that can be described more easily:
 
g = function(x) sin(x)
integral(sin, 0,pi/2); integral(sin, 0,pi); integral(sin, 0,2*pi)
#           1                2                2.032977e-16  ("pratically" 0)
 
# See here for a more a more complex example.
 
# Definite integrals of many funcions can be calculated with the areaF command: see.
# Only when we have focused the theorem fundamental theorem of calculus (see here or
# here) we can interpret areas as defined integrals.
 
#
# Gintegral(f,A,B,col)  traces the graph of x -> [A,x]f from A to B (A and B must be
# finite);  Gintegra, Gintegr plot slimmer graphics.
 
Plane(-pi,2*pi, -1,1); graph1(cos, -pi,2*pi, "blue")
Gintegr(cos, 0,2*pi, "red"); Gintegr(cos, 0,-pi, "violet")
POINT(0,0, "red")
    
 
# The graph on the right:
 
Plane(-pi,2*pi, -1,1); graph1(cos, -pi,2*pi, "blue")
POINT(2.5,0.5, "red")
GintegrK(cos, 2.5,2*pi, 0.5, "brown"); GintegrK(cos, 2.5,-pi, 0.5, "seagreen")
 
# With GintegralK(f,A,B, k, col) (and GintegraK, GintegrK) I can trace the graph
# of x -> k + [A,x]f
 
# In very special cases Gintegral (Gintegra, Gintegr) may not work; in such cases you
# can try using GintegralA or, if it is not working, GintegralB or … GintegralD  (or
# GintegraA, GintegrA, … for slimmer graphics):
f = function(x) ifelse(x<13, 7-6/3*(x-10), ifelse(x<17, 1+8/4*(x-13), 9-1/2*(x-17)))
Plane(10,19, 0,50); Gintegra(f, 10,19, "blue")
# Error in integrate
Plane(10,19, 0,50); GintegraA(f, 10,19, "brown")
 
          
# For the same reasons if integral does not work you can use integralA, …, intergralD
#
# With R we can compute definite integrals. We can also calculate the indefinite
# integrals of polynomial functions.
Q = c(-2,-3,2.4,-1/5)                       # the polynomial -2-3*x+2.4*x^2-1/5*x^3
q = function(x) Q[1]+x*Q[2]+x^2*Q[3]+x^3*Q[4]
# the integral:
I = function(p) {I=NULL; for (i in 1:length(p)) I = c(I,p[i]/i); I}
I(Q); fraction(I(Q))
#  -2.00 -1.50  0.80 -0.05
#     -2  -3/2   4/5 -1/20                   -2*x-3/2*x^2+4/5*x^3-1/20*x^4 + constant
Iq = function(x) x*I(Q)[1]+x^2*I(Q)[2]+x^3*I(Q)[3]+x^4*I(Q)[4]
# ∫[-1,2] q  using I(Q)
Iq(2)-Iq(-1)
#  -4.05
# I can calculate the defined integer in each case, with:
integral(q, -1,2)
#  -4.05
 
# For other functions I can use WolframAlpha. Let's see how to use it in this case:
#  integrate -2-3*x+2.4*x^2-1/5*x^3 dx
#  -0.05 x^4 + 0.8 x^3 - 1.5 x^2 - 2 x + constant
#  integrate -2-3*x+2.4*x^2-1/5*x^3 dx from x=-1 to 2
#  -0.45
#
# See here for the "study of a function"
 
           
 
# It is easy to integrate numerically and graphically functions that are not easy to
# study "by hand". See here for an example.
# An example of integral that you can not study by hand:
F = function(x) sqrt(3*x-x^3); integral(F, 0,3/2)
# 1.720722
 

# Note 1. If the integral between -Inf and Inf is 0 check whether the integer between 0
# and Inf is finite too: if it does not happen it means that the integral diverges (0
# is output only for a formal simplification).
 
# Note 2. To integrate a function between a and b if it is not defined in a point c
# between a and b, make the sum of the integral between a and c and that between c and b.
#
# See here  how to calculate areas using integrals
 
         
 
# If F is a function of at least 2 variables, if two of these are u and v,
# derivxy(F,u,v) computes the second partial derivative with respect to u and v,
# derivyx(F,u,v) that with respect to v and u (xy and yx indicate that the derivatives
# are "crossed", even if the variables are not x and y). deriv2(F,u) computes the 2nd
# derivative respect to u.
# If F is a function of x and y, hessian(F,a,b) (see) computes the hessian of F in (a,b).
# (we remember that if in (x,y) all partial derivatives are zero and the hessian is
#  negative/null/positive, then  it is a saddle point / we cannot conlude anything /
#  it is a min/max point if deriv2 respect to x is > / < 0 )
 
F = function(x,y,z) 5*x^2*y^3+x^3*z
deriv(F,"x"); deriv(F,"z")
# 5*(2*x)*y^3+3*x^2*z    x^3
deriv2(F,"y"); deriv2(F,"z")
# 5*x^2*(3*(2*y))        0
derivxy(F,"x","y"); derivxy(F,"x","z")
# 5*(2*x)*(3*y^2)        3*x^2
#
Q = function(x,y) x^4+6*x^2*y^2+y^4-32*(x^2+y^2)
                                  
# deriv(Q,"x") and deriv(Q,"y") are 0 in (0,0).
# To calculate hessian the variables must be x and y
hessian(Q,0,0); hessian(Q,0,4); hessian(Q,2,-2)
# 4096   16384  -8192
deriv2(Q,"x")
# 4 * (3 * x^2) + 6 * 2 * y^2 - 32 * 2
# In (0,0) the hessian in > 0 and deriv2 respect to x is > 0: max point
# And indeed:
MAXF2(Q, -1,1, -1,1)
#  0 0  0    In (0,0) I have a max, and Q(0,0) = 0
#
# For Taylor polynomials see here.
                                
#
# For integrals of the functions of two variables we can use the command INTEGRAL. See here.
   
# For max and min of functions of two variables see here (commands MINF2, MAXF2).

# The mean value (not of a set of values but) of a variable quantity depends on
# which variable one wants to make it dependent on.
# If f is a continuous positive function on [a,b] we assume as mean of f on [a,b] the
# value m such that the area of the drawn rectangle equals the result of integrating f
# on [a,b]:  m = integral(f, a,b) / (b-a)
                         
# For example, the problem:
#     Calculate the average value of the chord of a semi-circle
#          (of radius 1) that is parallel to the diameter
# depends on whether, for example, we identify the position of the chord using
# L, D or A. I can obtain 1, π/2 or 4/π.
                                  
FL = function(x) x; I = integral(FL, 0,2) / 2; I
# 1
FD = function(x) 2*sqrt(1-x^2); I = integral(FD, 0,1)/1; I
# 1.570796
pi/2
# 1.570796
FA = function(t) 2*sin(t); I = integral(FA, 0,pi/2) / (pi/2); I
# 1.27324
I*pi
# 4

Other examples of use