Num Fun!
Numerical Analysis in R - Sorting, Searching, Root Finding and Integration
CSJP
Sorting Begins …
Selection Sort
set.seed(1234)
a = sample(10)
a
## [1] 2 6 5 8 9 4 1 7 10 3
selsort = function(x) {
for (i in 1:(length(x) - 1)) {
min = i
for (j in (i + 1):length(x)) if (x[j] < x[min])
min = j # find the ith smallest number
t = x[min] # make exchanges btw ith smallest one and ith position element
x[min] = x[i]
x[i] = t
}
x
}
selsort(a)
## [1] 1 2 3 4 5 6 7 8 9 10
Searching Steps in …
Binary Search
set.seed(1234)
v = sample(1:26, 1)
v
## [1] 3
binsearch = function(x) {
p = 1
r = length(LETTERS[1:26])
if (x <= 0 | x > r)
return(print("No Such LETTER!"))
while (r >= p) {
h = floor((p + r)/2) # divide and conquer
if (x < h)
r = h - 1 # fall in lower half
else p = h + 1 # fall in upper half
if (x == h)
return(LETTERS[h]) # find it
}
}
binsearch(v)
## [1] "C"
v = 30
binsearch(v)
## [1] "No Such LETTER!"
v = 26
binsearch(v)
## [1] "Z"
v = 14
binsearch(v)
## [1] "N"
Finding Roots …
Plot Function f(x)
# Plot Function f(x)
f = function(x) x^3 - 9 * x
curve(f, -5, 5, ylab = expression(f(x) == x^3 - 9 %*% x), main = expression(paste("Plot of ",
f(x) == x^3 - 9 %*% x)))
abline(h = 0, v = 0)
grid()
points(c(-3, 0, 3), c(0, 0, 0), pch = 20)
text(c(-3, -0.35, 3), c(8, -8, -8), c("(-3, 0)", "(0, 0)", "(3, 0)"))
Root Finding (1) …
Bisect Method
# Bisect method
bisect = function(f, x0, x1, tol, itmax) {
if (sign(f(x0)) == sign(f(x1)))
return(print("Wrong Initial Conditions!"))
for (i in 1:itmax) {
x = (x0 + x1)/2
if (abs(f(x)) < tol)
return(x)
if (sign(f(x)) == sign(f(x0)))
x0 = x else x1 = x
}
return(print("Root Not Found!"))
}
bisect(f, -4, -2, 1e-05, 100)
## [1] -3
bisect(f, 2, 4, 1e-05, 100)
## [1] 3
Root Finding (2) …
Secant Method
# Secant method
f = function(x) x^3 - 9 * x
secant = function(f, x0, x1, tol, itmax) {
for (i in 1:itmax) {
x = x1 - f(x1) * (x1 - x0)/(f(x1) - f(x0))
if (abs(f(x)) < tol)
return(x)
x0 = x1
x1 = x
}
return(print("Root Not Found!"))
}
secant(f, -4, -2, 1e-05, 100)
## [1] -3
secant(f, 2, 4, 1e-05, 100)
## [1] 3
Root Finding (3) …
Newton Method
# Newton method
f = expression(x^3 - 9 * x)
fprime = D(f, "x")
newton = function(f, fprime, tol, itmax) {
for (i in 1:itmax) {
x = x - eval(f)/eval(fprime)
if (abs(eval(f)) < tol)
return(x)
}
return(print("Root Not Found!"))
}
x = -2
newton(f, fprime, 1e-05, 100)
## [1] -3
x = 2
newton(f, fprime, 1e-05, 100)
## [1] 3
Finally, Integration Comes …
Using R function
# Using R function
f = function(x) 1/x
fint = integrate(f, 1, 2)
eval(fint)
## 0.6931 with absolute error < 7.7e-15
Simple Quadrature Method
# simple quadrature method
f = function(x) 1/x
intrect = function(f, a, b, N) {
r = 0
w = (b - a)/N
for (i in 1:N) r = r + w * f(a - w/2 + i * w)
return(r)
}
intrect(f, 1, 2, 10)
## [1] 0.6928
intrect(f, 1, 2, 1000)
## [1] 0.6931
my sessionInfo …
print(sessionInfo(), locale = FALSE)
## R version 3.0.2 (2013-09-25)
## Platform: i386-w64-mingw32/i386 (32-bit)
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] cluster_1.14.4 foreach_1.4.1 plyr_1.8
## [4] reshape2_1.2.2 caret_6.0-21 lattice_0.20-24
## [7] randomForest_4.6-7 e1071_1.6-2 class_7.3-9
## [10] nnet_7.3-7 ada_2.0-3 rpart_4.1-4
## [13] gtools_3.2.1 gdata_2.13.2 caTools_1.16
## [16] KernSmooth_2.23-10 MASS_7.3-29 ROCR_1.0-5
## [19] gplots_2.12.1 ggplot2_0.9.3.1 scales_0.2.3
## [22] RColorBrewer_1.0-5 colorspace_1.2-4 vcd_1.3-1
## [25] boot_1.3-9 knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-6 car_2.0-19 codetools_0.2-8 dichromat_2.0-0
## [5] digest_0.6.4 evaluate_0.5.1 formatR_0.10 gtable_0.1.2
## [9] iterators_1.0.6 labeling_0.2 munsell_0.4.2 proto_0.3-10
## [13] stringr_0.6.2 tools_3.0.2