iCAL Modern Calculus with
R and Python


Version 2.0 developed from January 2022 for Cal III by Dr. Samuel Shen, Distinguished Professor
San Diego State University, California, USA
https://shen.sdsu.edu
Email:


Basics of R Programming

R As a Smart Calculator

#The text behind the symbol # is a comment for R 
#R can perform basic arithmetic calculations
1+4
## [1] 5
#Plotting a single data point
plot(2, col = 'red')

#pi is circumference ratio approximately 3.1415926 
2 + pi/4-0.8 
## [1] 1.985398
# <- is the assign symbol: assign 1 to x.
x<-1 
y<-2
z<-4
#Equivalent to 2*1ˆ2-4 =-2
t<-2*x^y-z 
t
## [1] -2
# Symbols "=" and "<-" are equivalent in most cases v=3
u=2 
v <- 3
u+v
## [1] 5
sin(u*v) # u*v = 6 is considered radian
## [1] -0.2794155

Write a Function in R

#Defining a function named 'square' that returns the square of the input
square <- function(x) x*x 
square(4)
## [1] 16
#Define function named 'sam'
sam <- function(x) {x*x + sin(x)}

#Calculate 'sam' at x=2 and plot
sam(2)
## [1] 4.909297
plot(sam, col = 'red', lwd = 5)

#Multivariable functoin
fctn <- function(x,y,z) {x+y-z/2} 
fctn(1,2,3)
## [1] 1.5

Plot with R

#Plot the curve of y=sin(x) from -pi to 2 pi 
plot(sin, -2*pi, 2*pi, col = 'blue', lwd =2)

#Plot the curve of y=sin(2x) from -9 to 9 given by 'x'
x = seq(-9, 9, length = 200)
y = sin(2*x)
plot(x, y,  col = 'red', lwd = 2, 
     pch = 16, cex = 0.5, type = 'l') 

#Different ways to plot single point
plot(1,2)

plot(c(1,2))

plot(c(1), c(2)) #the best

x = c(1)
y = c(2)
plot(x,y)

plot(c(1,2), c(3, 4))

x = c(1, 2)
y = c(3, 4)
plot(x,y, 
     xlim = c(-4, 4),
     ylim = c(-4, 4))

x = c(1, 2, 3, 4)
y = c(3, 4, -2, -4)
plot(x, y, 
     xlim = c(-4, 4),
     ylim = c(-4, 4))

#A sequence of points from 1 to 10
x = seq(1, 10)
x
##  [1]  1  2  3  4  5  6  7  8  9 10
y = x^2
y
##  [1]   1   4   9  16  25  36  49  64  81 100
plot(x, y)

n = 10
x = seq(1, n)
x
##  [1]  1  2  3  4  5  6  7  8  9 10
#Square each element of x
y = x^2
y
##  [1]   1   4   9  16  25  36  49  64  81 100
plot(x, y, type = "o")

#Define a function
square <- function(x) x*x 

#Plot the defined function
plot(square, -3, 2)

#Plot a 3D surface using the following code
x = y = seq(-1, 1, length=100)
Z = outer(x, y, function(x, y){1+x^2+y^2}) #outer (x,y, function) is outer product
# yields a 3D surface with perspective angle 310 deg
persp(x=x, y=y, z=Z, theta=310)

Symbolic Calculation for Calculus by R

#Take derivative of x^2 and the answer is 2x
D(expression(x^2,'x'), 'x')
## 2 * x
#Define a symbolic function with variable x D(fx,’x’)
fx = expression(x^2,'x') 
#Differentiate the function and yield a result below
D(fx, 'x')
## 2 * x
#Define a symbolic function with variable x D(fx,’x’)
xt = expression(2*t*cos(t),'t') 
#Differentiate the function and yield a result below
D(xt, 't')
## 2 * cos(t) - 2 * t * sin(t)
D(D(xt, 't'), 't') #second derivative
## -(2 * sin(t) + (2 * sin(t) + 2 * t * cos(t)))
sam = expression((1 + x^2)*sin(x^3),'x') 
sam
## expression((1 + x^2) * sin(x^3), "x")
#First derivative
samD1 = D(sam, 'x')
samD1
## 2 * x * sin(x^3) + (1 + x^2) * (cos(x^3) * (3 * x^2))
#Second derivative
samD2 = D(samD1, 'x')
samD2
## 2 * sin(x^3) + 2 * x * (cos(x^3) * (3 * x^2)) + (2 * x * (cos(x^3) * 
##     (3 * x^2)) + (1 + x^2) * (cos(x^3) * (3 * (2 * x)) - sin(x^3) * 
##     (3 * x^2) * (3 * x^2)))
#Third derivative
samD3 = D(samD2, 'x')
samD3
## 2 * (cos(x^3) * (3 * x^2)) + (2 * (cos(x^3) * (3 * x^2)) + 2 * 
##     x * (cos(x^3) * (3 * (2 * x)) - sin(x^3) * (3 * x^2) * (3 * 
##     x^2))) + (2 * (cos(x^3) * (3 * x^2)) + 2 * x * (cos(x^3) * 
##     (3 * (2 * x)) - sin(x^3) * (3 * x^2) * (3 * x^2)) + (2 * 
##     x * (cos(x^3) * (3 * (2 * x)) - sin(x^3) * (3 * x^2) * (3 * 
##     x^2)) + (1 + x^2) * (cos(x^3) * (3 * 2) - sin(x^3) * (3 * 
##     x^2) * (3 * (2 * x)) - ((cos(x^3) * (3 * x^2) * (3 * x^2) + 
##     sin(x^3) * (3 * (2 * x))) * (3 * x^2) + sin(x^3) * (3 * x^2) * 
##     (3 * (2 * x))))))
#Fourth derivative
D(samD3, 'x')
## 2 * (cos(x^3) * (3 * (2 * x)) - sin(x^3) * (3 * x^2) * (3 * x^2)) + 
##     (2 * (cos(x^3) * (3 * (2 * x)) - sin(x^3) * (3 * x^2) * (3 * 
##         x^2)) + (2 * (cos(x^3) * (3 * (2 * x)) - sin(x^3) * (3 * 
##         x^2) * (3 * x^2)) + 2 * x * (cos(x^3) * (3 * 2) - sin(x^3) * 
##         (3 * x^2) * (3 * (2 * x)) - ((cos(x^3) * (3 * x^2) * 
##         (3 * x^2) + sin(x^3) * (3 * (2 * x))) * (3 * x^2) + sin(x^3) * 
##         (3 * x^2) * (3 * (2 * x)))))) + (2 * (cos(x^3) * (3 * 
##     (2 * x)) - sin(x^3) * (3 * x^2) * (3 * x^2)) + (2 * (cos(x^3) * 
##     (3 * (2 * x)) - sin(x^3) * (3 * x^2) * (3 * x^2)) + 2 * x * 
##     (cos(x^3) * (3 * 2) - sin(x^3) * (3 * x^2) * (3 * (2 * x)) - 
##         ((cos(x^3) * (3 * x^2) * (3 * x^2) + sin(x^3) * (3 * 
##             (2 * x))) * (3 * x^2) + sin(x^3) * (3 * x^2) * (3 * 
##             (2 * x))))) + (2 * (cos(x^3) * (3 * (2 * x)) - sin(x^3) * 
##     (3 * x^2) * (3 * x^2)) + 2 * x * (cos(x^3) * (3 * 2) - sin(x^3) * 
##     (3 * x^2) * (3 * (2 * x)) - ((cos(x^3) * (3 * x^2) * (3 * 
##     x^2) + sin(x^3) * (3 * (2 * x))) * (3 * x^2) + sin(x^3) * 
##     (3 * x^2) * (3 * (2 * x)))) + (2 * x * (cos(x^3) * (3 * 2) - 
##     sin(x^3) * (3 * x^2) * (3 * (2 * x)) - ((cos(x^3) * (3 * 
##     x^2) * (3 * x^2) + sin(x^3) * (3 * (2 * x))) * (3 * x^2) + 
##     sin(x^3) * (3 * x^2) * (3 * (2 * x)))) - (1 + x^2) * (sin(x^3) * 
##     (3 * x^2) * (3 * 2) + ((cos(x^3) * (3 * x^2) * (3 * x^2) + 
##     sin(x^3) * (3 * (2 * x))) * (3 * (2 * x)) + sin(x^3) * (3 * 
##     x^2) * (3 * 2)) + (((cos(x^3) * (3 * (2 * x)) - sin(x^3) * 
##     (3 * x^2) * (3 * x^2)) * (3 * x^2) + cos(x^3) * (3 * x^2) * 
##     (3 * (2 * x)) + (cos(x^3) * (3 * x^2) * (3 * (2 * x)) + sin(x^3) * 
##     (3 * 2))) * (3 * x^2) + (cos(x^3) * (3 * x^2) * (3 * x^2) + 
##     sin(x^3) * (3 * (2 * x))) * (3 * (2 * x)) + ((cos(x^3) * 
##     (3 * x^2) * (3 * x^2) + sin(x^3) * (3 * (2 * x))) * (3 * 
##     (2 * x)) + sin(x^3) * (3 * x^2) * (3 * 2)))))))
#Derivative of e^x
patient = expression(exp(x), 'x')
patient
## expression(exp(x), "x")
D(patient, 'x')
## exp(x)
#Change the expression and use the same derivative command
fx = expression(x^2*sin(x),'x')
D(fx,'x')
## 2 * x * sin(x) + x^2 * cos(x)

Symbolic Calculation Method 2

#by the library mosaicCalc
#install.packages('mosaicCalc')
library(mosaicCalc)
D(x^2 ~ x)
## function (x) 
## 2 * x
antiD(2*x ~ x)
## function (x, C = 0) 
## x^2 + C
D(2*t*cos(2*t) ~ t)
## function (t) 
## {
##     .e1 <- 2 * t
##     2 * cos(.e1) - 4 * (t * sin(.e1))
## }
F = antiD(x ~x)
F(1)-F(0)
## [1] 0.5
fxy = expression(3 - x^2 - 2*y^2, 'x','y')
fxy
## expression(3 - x^2 - 2 * y^2, "x", "y")
D(fxy, 'x')
## -(2 * x)
D(fxy, 'y')
## -(2 * (2 * y))

More variables

#define a function of two or more variables
fxyz = expression(x^2+y^2+z^2, 'x','y','z')
#This gives the expression of the function in terms of x, y and z 
fxyz 
## expression(x^2 + y^2 + z^2, "x", "y", "z")
#This gives the partial derivative with respect to x
D(fxyz,'x') 
## 2 * x
#This gives the partial derivative with respect to y
D(fxyz,'y')
## 2 * y

Exercise 7.12

#define a function of two or more variables
u = expression(sin(x - c*t), 'x','t')
#This gives the expression of the function in terms of x, y and z
u
## expression(sin(x - c * t), "x", "t")
#This gives the partial derivative with respect to x
ux = D(u,'x') 
ux
## cos(x - c * t)
#This gives the partial derivative with respect to t
ut = D(u,'t')
ut
## -(cos(x - c * t) * c)
# functions with constants 
fxa = expression(a*x^2, 'x')
D(fxa, 'x')
## a * (2 * x)
fxn = expression(x^n, 'x')
D(fxn, 'x')
## x^(n - 1) * n
sxan = expression(a*sin(n*x), 'x')
sxan
## expression(a * sin(n * x), "x")
D(sxan, 'x')
## a * (cos(n * x) * n)

Partial derivatives in Midterm 2 Problem 1

#define a function of two or more variables
fxyt = expression(3*cos(3*x + y)*exp(-0.2*t), 'x','y','t')
#This gives the expression of the function in terms of x, y and t
fxyt 
## expression(3 * cos(3 * x + y) * exp(-0.2 * t), "x", "y", "t")
#This gives the partial derivative with respect to x
D(fxyt,'x')
## -(3 * (sin(3 * x + y) * 3) * exp(-0.2 * t))
#This gives the partial derivative with respect to y 
D(fxyt,'y')
## -(3 * sin(3 * x + y) * exp(-0.2 * t))
#This gives the partial derivative with respect to t
D(fxyt,'t')
## -(3 * cos(3 * x + y) * (exp(-0.2 * t) * 0.2))
square = function(x) x^2
integrate (square, 0,1) 
## 0.3333333 with absolute error < 3.7e-15
sam2 = function(x){(1 + x^2) * sin(x^3)}
sam2
## function (x) 
## {
##     (1 + x^2) * sin(x^3)
## }
#Integrate xˆ2 from 0 to 1 equals to 1/3 with details below
integrate(sam2, 0,1)
## 0.3870778 with absolute error < 4.3e-15
fc = function(x){cos(x)}
integrate(fc, 0, pi/2)
## 1 with absolute error < 1.1e-14
#Integrate cos(x) from 0 to pi/2 equals to 1 with details below
integrate(cos,0,pi/2) 
## 1 with absolute error < 1.1e-14
#Symbolic expression() cannot be used for numerics
a =2; n = 3
integrate(sxan,0, pi/2) 
## Error in get(as.character(FUN), mode = "function", envir = envir): object 'sxan' of mode 'function' was not found
#Use function for numerics
sxan2 = function(x){a*sin(n*x)}
integrate(sxan2, 0, pi/2) 
## 0.6666667 with absolute error < 2.2e-14

Antiderivative

#install.packages('mosaicCalc')
library(mosaicCalc)
antiD(x^2 ~ x)
## function (x, C = 0) 
## x^3/3 + C
antiD(sin(2*x) ~ x)
## function (x, C = 0) 
## C - cos(-2 * x)/2
sam = expression(x^2, 'x')
antiD(sam ~ x) #not working
## function (x, C = 0, sam) 
## sam * x + C
sam = function(x) x^2
antiD(sam ~ x) # not working
## function (x, C = 0, sam) 
## sam * x + C
antiD((sin(2*x))^2 ~ x) #not working
## function (x, C = 0) 
## (4 * x + sin(-4 * x))/8 + C

Vectors and Matrices

#Enter data inside c() for a 4X1 column vector 
c(1,6,3,pi,-3) 
## [1]  1.000000  6.000000  3.000000  3.141593 -3.000000
#Generate a sequence from 2 to 6
seq(2,6) 
## [1] 2 3 4 5 6
# Generate a sequence from 1 to 10 with 2 increment
seq(1,10,2) 
## [1] 1 3 5 7 9
x=c(1,-1,1,-1)
#1 is added to each element of x
x+1 
## [1] 2 0 2 0
#2 multiplies each element of x
2*x 
## [1]  2 -2  2 -2
#Each element of x is divided by 2
x/2
## [1]  0.5 -0.5  0.5 -0.5
y=seq(1,4)
# This * multiples each pair of elements
x*y
## [1]  1 -2  3 -4
# Transpose of a matrix
t(x) 
##      [,1] [,2] [,3] [,4]
## [1,]    1   -1    1   -1
#Matrix multiplication: 1X4 matrix times a 4X1 matrix 
#This is equivalent to a dot product
t(x)%*%y 
##      [,1]
## [1,]   -2
u = c(3, 4)
v = c(-2.5, 1.5)
nu = norm(u, type = '2') #type = '2' specifies Euclidean norm calculation
nv = norm(v, type = '2') 
dotuv = sum(u*v)
c = dotuv/(nu*nv)
c
## [1] -0.1028992
acos(c)*180/pi
## [1] 95.90614
u = c(1,1,0)
v = c(1,-1,0)
nu = norm(u, type = '2')
nv = norm(v, type = '2') 
dotuv = sum(u*v)
c = dotuv/(nu*nv)
c
## [1] 0
acos(c)*180/pi
## [1] 90
acos(sum(u*v)/(norm(u, type = '2')*norm(v, type = '2')))*180/pi
## [1] 90
#Cross product, dot product
#Import the required library 
#install.packages('pracma')
library(pracma) 
## 
## Attaching package: 'pracma'
## The following objects are masked from 'package:Matrix':
## 
##     expm, lu, tril, triu
# Taking two vectors 
a = c(3, 5, 4) 
b = c(2, 7, 5) 
# Calculating cross product using cross() 
cross(a, b)
## [1] -3 -7 11
dot(a,b)
## [1] 61
#4X1 matrix times a 1X4 matrix yields a 4X4 matrix
x%*%t(y) 
##      [,1] [,2] [,3] [,4]
## [1,]    1    2    3    4
## [2,]   -1   -2   -3   -4
## [3,]    1    2    3    4
## [4,]   -1   -2   -3   -4
#Convert a vector y into a matrix of 2 rows. 
mx=matrix(x,2)
my=matrix(y,2)
#The matrix elements go by column, first column, second, etc
my
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4
#find dimensions of a matrix

dim(my) 
## [1] 2 2
#Convert a matrix to a vector, also via columns
as.vector(my) 
## [1] 1 2 3 4
#multiplication between each pair of elements
mx*my 
##      [,1] [,2]
## [1,]    1    3
## [2,]   -2   -4
#division between each pair of elements
mx/my  
##      [,1]       [,2]
## [1,]  1.0  0.3333333
## [2,] -0.5 -0.2500000
#Subtraction between each pair of elements and scalar multiplication of each element of my
mx-2*my
##      [,1] [,2]
## [1,]   -1   -5
## [2,]   -5   -9
#matrix multiplication 
mx%*%my 
##      [,1] [,2]
## [1,]    3    7
## [2,]   -3   -7
#determinant of a square matrix
det(my) 
## [1] -2
#Inverse of a matrix
myinv = solve(my) 
myinv
##      [,1] [,2]
## [1,]   -2  1.5
## [2,]    1 -0.5
#verifies the inverse of a matrix 
#Should result in the identity matrix
myinv%*%my 
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
#output the diagonal vector of a matrix
diag(my) 
## [1] 1 4
#yields eigenvalues and unit eigenvectors
myeig=eigen(my)

#$values gives the two eigenvalues
#$vectors gives the two eigenvectors
myeig
## eigen() decomposition
## $values
## [1]  5.3722813 -0.3722813
## 
## $vectors
##            [,1]       [,2]
## [1,] -0.5657675 -0.9093767
## [2,] -0.8245648  0.4159736
#output only eigenvalues 
myeig$values 
## [1]  5.3722813 -0.3722813
#output only eigenvectors
myeig$vectors
##            [,1]       [,2]
## [1,] -0.5657675 -0.9093767
## [2,] -0.8245648  0.4159736
#SVD decomposition of a matrix M=UDV’ 
#SVD can be done for any m-by-n rectangular matrix 
mysvd = svd(my) 

#output d, U, and V
mysvd 
## $d
## [1] 5.4649857 0.3659662
## 
## $u
##            [,1]       [,2]
## [1,] -0.5760484 -0.8174156
## [2,] -0.8174156  0.5760484
## 
## $v
##            [,1]       [,2]
## [1,] -0.4045536  0.9145143
## [2,] -0.9145143 -0.4045536
#output d only, as a vector 
mysvd$d 
## [1] 5.4649857 0.3659662
#Output the U matrix 
U=mysvd$u 
#Generate the D matrix
D=diag(mysvd$d) 

#Output the V matrix
V=mysvd$v 
#Recover the original matrix 
U%*%D%*%t(V) 
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4
my 
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4
#solve linear equations matrix %*% x = b 
ysol=solve(my,c(1,3)) 
#the resulting solution of solve(matrix, b)
ysol 
## [1]  2.5 -0.5
#verifies the solution
my%*%ysol 
##      [,1]
## [1,]    1
## [2,]    3

Statistics

#generate 10 normally distributed numbers 
x=rnorm(10) 
x
##  [1] -0.33269498  1.18688534 -0.15674479 -0.54297390  0.59975444  0.89181936
##  [7]  0.98874502  0.09467425 -1.20516875  0.20222203
#Calculate the mean
mean(x)
## [1] 0.1726518
#Calculate the variance
var(x)
## [1] 0.5750709
#Calculate the standard deviation
sd(x)
## [1] 0.7583343
#Calculate the median
median(x)
## [1] 0.1484481
#Returns the values for the minimum, 1st quartile, 2nd quartile, 3rd quartile, and maximum
quantile(x)
##         0%        25%        50%        75%       100% 
## -1.2051688 -0.2887074  0.1484481  0.8188031  1.1868853
#Yields the min and max of x
range(x) 
## [1] -1.205169  1.186885
#Maximum value
max(x)
## [1] 1.186885
#Yields the box plot of x
boxplot(x) 

#generate 1000 normally distributed random numbers N(0,1) 
w=rnorm(1000) 

#generate 100 random numbers following N(10,5ˆ2)
z=rnorm(10000, mean=10, sd=5) 

#statistical summary of the data sequence
summary(rnorm(12)) 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.3812 -0.2443  0.2995  0.2268  0.6379  2.3073
#plot the histogram of 1000 random numbers
hist(w)

R Tutorials

#Commands to input data files

# mydata <- read.csv("mydata.csv") 
# read csv file named "mydata.csv"

#mydata <- read.table("mydata.txt") 
# read text file named "my data.txt"

#library(gdata) 
# load the gdata package 

#mydata = read.xls("mydata.xls") 
# read an excel file

#library(foreign) 
# load the foreign package 

#mydata = read.mtp("mydata.mtp") 
# read from .mtp file

#library(foreign) 
# load the foreign package 

#mydata = read.spss("myfile", to.data.frame=TRUE)
#read a spss file

#ff <- tempfile()
#cat(file = ff, "123456", "987654", sep = "\n") 
#read.fortran(ff, c("F2.1","F2.0","I2"))
#read a fotran file

# library(ncdf4)
# load the ncdf4 package 

#ncin <- ncdf4::nc_open("ncfname") 
# open a NetCDF file

#lon <- ncvar_get(ncin, "lon") 
#read data "lon" from a netCDF file into R

#library("rjson")
# load the rjson package 

#jd<- fromJSON(file = "dat.json") 
# read data from a JSON file dat.json


#library(ncdf4) 
#Suppose the following error message pops up
#When you load ncdf4 package
#Error in library(ncdf4) : there is no package called ncdf4 

#You can install the R package by
#install.packages("ncdf4")

#R packages: animation, chron, e1071, fields, ggplot2, lattice,
#latticeExtra, maps, mapdata, mapproj, matrixStats, ncdf4,
#NLRoot, RColorBrewer, rgdal, rasterVis, raster, rjson, sp, TTR

#The zipped data for this book can be downloaded from: 
#climatemathematics.sdsu.edu/data.zip

#To load a single package, such as "animation", you can do
#library(animation)

#You can also load all these packages in one shot using
# pacman
#install.packages("pacman")
#library(pacman)
#pacman::p_load(animation, chron, e1071, fields, ggplot2, lattice,
               #latticeExtra, maps, mapdata, mapproj, matrixStats, ncdf4,
               #NLRoot, RColorBrewer, rgdal, rasterVis, raster, rjson, sp, TTR)