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. Volumes.
# 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)
#
# If f is a piecewise-defined function we must derive every sub-function (see).
#
# Note. df calculated by the software can have a wider domain than f. Example:
f = function(x) log(1+tan(x/2))
BF=5; HF=3; a=-2*pi; b=2*pi; Plane(a,b, -5,5)
graph2(f, a,b, "brown"); df = function(x) eval( deriv(f,"x") ); deriv(f,"x")
# 1/2/cos(x/2)^2/(1 + tan(x/2))
graph2(df, a,b, "seagreen")
FDOM=f; graph(DOM, a,b, "red") # DOM is defined 0 in the domain of FDOM
pointO(pi*c(-1,-1/2,1,2-1/2),rep(0,4),"blue")
text(-5,-1/2,"domain",font=2,col="red"); text(-1/2,-1.5,"f",font=2,col="brown")
text(-1/2,2.5,"df",font=2,col="seagreen"); text(-2.5,-1.5,"?",font=2,col="seagreen")
#
# 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
# We can use areaPol also for spline functions (see):
x = c(-0.2, 1.3, 1.7, 3.2, 4.1, 4.6)
y = c(9.1, 17.5, 20.7, 25.4, 35.7, 51.2)
Plane(-1,6,-10,70); POINT(x,y,1); L=spline(x,y,1000); polyline(L$x,L$y,"red")
INTEGRAL(L$x,L$y) # INTEGRAL defined above
# 108.6348
# 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)
# Here, how to calculate the effective force on the mast of a racing sailbot
# (from "Numerical Methods", Chapra & Canale, McGraw-Hill)
# See here for integrals of F from A to B when A or B are infinite or F is undefined
# at A or B .
# 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
# Gintegral0 (or GintegraA, GintegrA, … Gintegra0, Gintegr0, 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
# or intergral0.
#
# In some case, some manipulations are necessary to achieve a result: a simple example.
#
# With R we can compute definite integrals. We can also calculate the indefinite
# integrals (ie antiderivatives) 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
#
# The elementary functions (functions of one variable which are the composition of
# arithmetic operations, powers, exponentials, logarithms, trigonometric functions,
# inverse trigonometric functions) are closed under differentiation but not under
# integration. Some examples of elementary functions that have not an elementary
# antiderivative: sqrt(1+x^3) sqrt(1-x^4) 1/sqrt(1+x^4) 1/log(x) log(log(x))
# exp(x^2) exp(-x^2) exp(x)/x exp(exp(x)) x^2*exp(-x^2) sin(x^2) cos(exp(x))
# sin(x)/x. Hence the importance of the numerical integration we have seen.
#
# We can also trace the direction field of the solutions of an indefinite integral (see):
#
# See here for the compute of limits of a function
#
# See here for the "study of a function"
# See here for the correct use of x -> x^x and similar function that use ^.
| x -> x^x | |  |
# 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
#
# The gradient of Q, grad(Q), is x,y -> ( d(Q(x,y)/dx, d(Q(x,y)/dy ):
c( deriv(Q,"x"), deriv(Q,"y") )
# 4*x^3+12*x*y^2-64*x 12*x^2*y+4*y^3-64*y
# For the use of this concept see here.
#
# For Taylor polynomials see here.
#
# For integrals of the functions of two variables we can use the command INTEGRAL. See here.
# For volumes of solids of revolution 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