NUMERICAL ANALYSIS
SALVADOR CASTRO
Numerical analysis is a field of mathematics and computer science which studies algorithms for obtaining numerical approximations to problems involving continuous variables which do not have an easy algebraic solution.
1 Locating Roots of Equations
A root-finding algorithm is a numerical method, or algorithm, for finding a value \(x\) such that \(f(x) = 0\), for a given function \(f\). Such an \(x\) is called a root of the function \(f\).
1.1 Plotting the Roots of Functions
INPUT: Function \(f\), endpoint values \(from, to,\) root(s) of the function f
OUTPUT: Plot of f between the end points \(from, to,\) and the location of its roots
plot.root <- function(f, from, to, root = NULL, title = NULL) {
library(shape, quietly = TRUE, warn.conflicts = FALSE)
x <- seq(from, to, length.out = 100)
y <- f(x)
plot(x, y,
axes = FALSE,
xlim = c(from, to),
ylim = c(floor(min(f(x))), ceiling(max(f(x)))),
xlab = "x",
ylab = "y",
main = title,
col.main = "tomato",
type = "l",
lwd = 5,
col = "tomato")
# x-axis
axis(1, pos = 0)
# y-axis
axis(2, pos = 0)
# coordinate point
points(root, rep(0, length(root)),
col = "white",
bg = "tomato",
pch = 21,
cex = 1)
# root
axis(1, at = round(root, 2),
pos = (ceiling(max(f(x))) - floor(min(f(x))))/10,
tick = FALSE,
col.axis = "darkred")
} # end plot.root
# save to the working directory getwd(), ls()
dump("plot.root", file = "plot.root.R")
1.1.1 Plot the roots of \(x^2-1\) on \([-2, 2]\)
Solution:
\[ \begin{aligned} x^2-1 &= 0 \\ x^2 &= 1 \\ \sqrt{x^2} &= \sqrt{1} \\ x &= \pm 1 \\ \end{aligned} \]
plot.root(f = function(x){x^2 - 1},
from = -2, to = 2, root = c(-1, 1),
title = expression(paste("y = " * x^2)))
1.2 The Bisection Method
The bisection method in mathematics is a root-finding method that repeatedly bisects an interval and then selects a subinterval in which a root must lie for further processing. It is a very simple and robust method, but it is also relatively slow.
As described in (Kincaid and Cheney 2008, 76):
INPUT: Function \(f\), endpoint values \(a, b,\) tolerance epsilon, maximum iterations \(NMAX\)
CONDITIONS: \(a < b\), \(f(a)\) and \(f(b)\) have opposite signs i.e. \(f(a) \cdot f(b) < 0\)
OUTPUT: value which differs from a root of \(f(x) = 0\) by less than epsilon
bisection <-
function(f, from, to, epsilon = 1e-6, NMAX = 101, plot = FALSE, title = NULL) {
source("plot.root.R")
a <- from
b <- to
error <- b - a
fa <- f(a)
fb <- f(b)
for (n in 1:NMAX) { # limit iterations to prevent infinite loop
error <- error/2
c <- a + error
fc <- f(c)
if (abs(error) < epsilon) {
if (plot)
plot.root(f, from, to, c, title)
return(c)
} # end if
if (sign(fa) != sign(fc)) {
b <- c
fb <- fc
}
else {
a <- c
fa <- fc
} # end if
} # end for n
} # end bisection
# save to the working directory getwd(), ls()
dump("bisection", file = "bisection.R")
1.2.1 The root of \(\mathrm{cosine}\) function on \(\left[\frac{-\pi}{3}, \frac{4\pi}{3}\right]\)
bisection(f = function(x){return(cos(x))}, from = -pi/3, to = 4*pi/3, plot = TRUE, title = expression(paste("y = " * cos(x))))
## [1] 1.5707971.2.2 The roots of a quadratic function
Bisection algorithm ensures to locate only one root of the equation at a time given it exists between the initial interval provided by the user. When an equation has more than one roots, the choice of initial interval determines which root is located.
f <- function(x) x^2 - 1
(r1 <- bisection(f, from = -2, to = 0))## [1] -1.000001(r2 <- bisection(f, from = 0, to = 2))## [1] 0.999999plot.root(f, from = -2, to = 2, root = c(r1, r2), title = expression(paste("y = " * x^2)))
1.2.3 The root of \(\displaystyle{\ln \left(\frac{1+x}{1-x^2} \right)}\) on [-0.5, 0.5]
bisection(function(x){log((1 + x )/(1 - x^2))}, from = -.5, to = .5, plot = TRUE,
title = expression(paste("y = ln", bgroup("(",frac(1+x, 1-x^2),")")), line = 1.6))
## [1] -9.536743e-07
1.2.4 Board in hall problem
Taken from (Kincaid and Cheney 2008, 89):
In a building, two intersecting halls with widths \(\omega_1 = 9\) feet and \(\omega_2 = 7\) feet meet at an angle \(\alpha = 125^{\circ}\), as shown:
Board in hall problem
Assuming a two-dimensional situation, what is the longest board that can negotiate the turn? Ignore the thickness of the board. The relationship between the angles \(\Theta\) and the length of the board \((\mathscr{l} = \mathscr{l}_1 + \mathscr{l}_2)\) is
\[ \mathscr{l}_1 = \omega_1 \cdot \mathrm{csc}(\beta) \]
and
\[ \mathscr{l}_2 = \omega_2 \cdot \mathrm{csc}(\gamma) \]
where
\[ \beta = \pi - \alpha - \gamma \]
thus
\[ \mathscr{l} = \omega_1 \cdot \mathrm{csc}(\pi - \alpha - \gamma) + \omega_2 \cdot \mathrm{csc}(\gamma). \]
The maximum length of the board that can make the turn is found by minimizing \(\mathscr{l}\) as a function of \(\gamma\) . Taking the derivative and setting \(\frac{d\mathscr{l}}{d\gamma} = 0\), we obtain
\[ \omega_1 \cdot \mathrm{cot}(\pi - \alpha - \gamma) \cdot \mathrm{csc}(\pi - \alpha - \gamma) - \omega_2 \cdot \mathrm{cot}(\gamma) \cdot \mathrm{csc}(\gamma) = 0 \]
Substitute in the known values and numerically solve the nonlinear equation.
Solution:
\[ 9 \cdot \mathrm{cot}(\pi - \frac{125^{\circ}\cdot \pi}{180^{\circ}} - \gamma) \cdot \mathrm{csc}(\pi - \frac{125^{\circ}\cdot \pi}{180^{\circ}} - \gamma) - 7 \cdot \mathrm{cot}(\gamma) \cdot \mathrm{csc}(\gamma) = 0 \]
Approximate the value of angle \(\gamma\)
gamma <- bisection(function(gamma){9*1/tan(pi - (125*pi/180) - gamma)*1/sin(pi - (125*pi/180) - gamma) - 7*1/tan(gamma)*1/sin(gamma)}, from = 0, to = pi/4)
gamma## [1] 0.4511963then substitute it in the length equations, \(\mathscr{l}_1\) and \(\mathscr{l}_2\)
# l2
(l2 <- 7 * 1/sin(gamma))## [1] 16.05348# l1
(l1 = 9 * 1/sin(pi - (125*pi/180) - gamma))## [1] 18.47774# total length l
(l = l1 + l2)## [1] 34.53122
1.3 The Method of False Position (regula falsi)
The false position method, like the bisection method, trappes the root of an equation in a sequence of narrowing intervals. Unlike the bisection method, the method of false position does not selec the midpoint of each interval. Instead, it uses the point where the secant lines intersect the x-axis.
As described in (Burden and Faires 2011, 73):
INPUT: Initial approximations \(a, b\); tolerance epsilon; maximum number of iterations \(NMAX\) and optional plot which is set to FALSE by default
CONDITIONS: \(a < b\), \(f(a)\) and \(f(b)\) have opposite signs i.e. \(f(a) \cdot f(b) < 0\)
OUTPUT: Approximate solution \(p\) or message of failure.
regula.falsi <-
function(f, from, to, epsilon = 1e-6, NMAX = 101, plot = FALSE, title = NULL) {
source("plot.root.R")
i <- 2
a <- from
b <- to
fa <- f(a)
fb <- f(b)
while (i < NMAX) {
p <- b - fb*(b - a)/(fb - fa)
if (abs(p - b) < epsilon) {
if (plot)
plot.root(f, from, to, p, title)
return(p)
} # end if
i <- i + 1
q <- f(p)
if (q*fb < 0) {
a <- b
fa <- fb
} # end if
b <- p
fb <- q
} # end while
} # end regula.falsi
# save to the working directory getwd(), ls()
dump("regula.falsi", file = "regula.falsi.R")
1.3.1 The root of \(\mathrm{sine}\) function on \(\left[\frac{\pi}{16}, \; \frac{31\pi}{16}\right]\)
regula.falsi(f = function(x){sin(x)}, from = pi/16, to = 31*pi/16, plot = TRUE,
title = expression(paste("y = " * sin(x))))
## [1] 3.141593
1.4 Newton-Raphson Method
Newton-Raphson method in an iterative algorithm that starts with an initial guess which is reasonably close to the true root, then the approximation is refined by computing the tangent line of the function at the point that corresponds to the initial guess. The x-intercept of this tangent line will typically be a better approximation to the function’s root than the original guess, and the procedure is repeated until the approximation approaches the root sufficiently close. Convergence is assumed when the ratio of the function and its derivative is less than a predetermined margin of error.
As described in (Kincaid and Cheney 2008, 92):
INPUT: Function \(f\), its derivative \(f'\), endpoint values \(a, b,\) tolerance epsilon, delta, maximum iterations \(NMAX\) and optional plot which is set to FALSE by default
CONDITIONS: \(a < b\), \(f(a)\) and \(f(b)\) have opposite signs i.e. \(f(a) \cdot f(b) < 0\)
OUTPUT: Value which differs from a root of \(f(x) = 0\) by less than epsilon
Newton.Raphson <-
function(f, f_prime, x = 0, nmax = 101, epsilon = 1e-6, delta = .01) {
fx <- f(x)
for (n in 1:nmax) {
fp <- f_prime(x)
if (abs(fp) < delta)
return(cat("small derivative"))
d <- fx/fp
x <- x - d
fx <- f(x)
if (abs(d) < epsilon)
return(x)
} # end for
} # end Newton.Raphson
# save to the working directory getwd(), ls()
dump("Newton.Raphson", file = "Newton.Raphson.R")
1.4.1 The root of \(x^3 + 2x^2 + x - 3\) on \([1, 3]\)
\[ f(x) = x^3 + 2x^2 + x - 3 \]
and
\[ \frac{d}{dx} \; f(x) = 3x^2 + 4x + 1 \]
f <- function(x) ((x - 2)*x + 1)*x - 3
f_prime <- function(x) (3*x - 4)*x + 1
Newton.Raphson(f, f_prime)## [1] 2.174559
1.4.2 Computing square roots
Compute the square root of \(\frac{1}{2}\) by using Newton’s method. Note that any positive real number has two square roots, one positive and one negative.
Solution:
A square root of \(x\) is a number \(r\) such that \(r^2=x\). For \(x = \frac{1}{2}\) we have
\[ \begin{aligned} r &= \sqrt{\frac{1}{2}} \\ \\ r^2 &= \frac{1}{2} \\ \\ r^2-\frac{1}{2} &= 0. \end{aligned} \]
so wee need to find the roots of the function
\[ f(r) = r^2-\frac{1}{2}. \]
Its derivative is
\[ \frac{d}{dr} \; f(r) = 2r. \]
f <- function(r) r^2 - (1/2)
f_prime <- function(r) 2*r
(root1 <- Newton.Raphson(f, f_prime, x = 1))## [1] 0.7071068(root2 <- Newton.Raphson(f, f_prime, x = -1))## [1] -0.7071068plot.root(function(r){r^2 - (1/2)}, from = -1.5, to = 1.5, root = c(root1, root2),
title = expression(paste("y = " * r^2 - frac(1,2))))
1.5 Secant Method
Secant method is a variation of the Newton-Raphson method that avoids differentiation by using an approximation of the functions derivative:
\[ f'(x) \approx \frac{f(x+h) - f(x)}{h}, \] for a small \(h\)
As described in (Kincaid and Cheney 2008, 122):
INPUT: Function \(f\), endpoint values \(a, b,\) tolerance epsilon, maximum iterations \(NMAX\)
CONDITIONS: \(a < b\), \(f(a)\) and \(f(b)\) have opposite signs i.e. \(f(a) f(b) < 0\)
OUTPUT: Value which differs from a root of \(f(x) = 0\) by less than epsilon
secant <-
function(f, from, to, nmax = 101, epsilon = 1e-6, plot = FALSE, title = NULL) {
source("plot.root.R")
a <- from
b <- to
fa <- f(a)
fb <- f(b)
if (abs(fa) > abs(b)) {
temp <- a
a <- b
b <- temp
temp <- fa
fa <- fb
fb <- temp
} #end if
for (n in 2:nmax) {
if (abs(fa) > abs(fb)) {
temp <- a
a <- b
b <- temp
temp <- fa
fa <- fb
fb <- temp
} # end if
d <- (b - a)/(fb - fa)
b <- a
fb <- fa
d <- d*fa
if (abs(d) < epsilon) {
if (plot)
plot.root(f, from, to, root = a, title = title)
return(a)
}
a <- a - d
fa <- f(a)
} # end for
} # end procedure secant
# save to the working directory getwd(), ls()
dump("secant", file = "secant.R")
Solve \(x^3 - x - 1 = 0\)
secant(f = function(x){x^(3) - x - 1}, -1.5, 1.5, epsilon = 1e-6, plot = T, title = expression(paste("y = " * x^3 - x - 1)))
## [1] 1.324718
1.5.1 The arc-length of catenary by root finding
Taken from (Kincaid and Cheney 2008, 76):
An electric power cable is suspended (at points of equal height) from two towers that are 100 meters apart. The cable is allowed to dip 10 meters in the middle. How long is the cable?
Solution:
Catenary problem
Numerically solve for lambda: It is known that the curve assumed by a suspended cable is a catenary. When the y-axis passes through the lowest point, we can assume an equation of the form \[ y(x) = \lambda \cdot \cosh(x/\lambda). \]
The conditions of the problem are that \(y(50) = y(0) + 10\). Hence, we obtain
\[ \begin{aligned} \lambda \cdot \cosh(\frac{50}{\lambda}) &= \lambda + 10 \\ \\ \lambda \cdot \cosh(\frac{50}{\lambda}) - \lambda - 10 &= 0 \end{aligned} \]
(lambda <- secant(function(lambda){lambda*cosh(50/lambda) - lambda - 10}, from = 100, to = 200))## [1] 126.6324Substitute \(\lambda\) into the arc length formula of the catenary
\[ y(x) = \lambda \cdot \cosh(x/\lambda) \]
arc.length.catenary <- function(x) {return(lambda*sinh(x/lambda))}
2 * arc.length.catenary(50)## [1] 102.61871.6 Points of intersection of two graphs of functions
Taken from (Kincaid and Cheney 2008, 85):
Approximate where the graphs of \(y = 3^x\) and \(y = e^x\) intersect.
Soution:
Graphs intersect when \(3x = e^x\) or \(3x - e^x = 0\)
f <- function(x){3*x - exp(x)}
(root1 <- regula.falsi(f, from = 0, to = 1))## [1] 0.6190616(root2 <- regula.falsi(f, from = 1, to = 2))## [1] 1.512134plot.root(f, from = 0, to = 2, root = c(root1, root2), title = expression(paste("y = ", e^x - 3*x)))
Plots of \(f(x) = e^x\) and \(f(x) = 3x\)
# x continuum
x <- seq(0.3, 1.8, by = .1)
# plot# e^x
plot(x, exp(x),
type = "l",
lwd = 5,
ylab = "f(x)",
col = "tomato",
xlim = c(0, 2),
ylim = c(0, 7))
# plot 3x
lines(x, 3*x,
lwd = 5,
col = "dodgerblue")
# tracing lines
segments(x0 = c(root1, root2), y0 = c(3*root1, exp(root2)),
x1 = c(root1, root2), y1 = -.8,
lty = 3,
lwd = .5)
# arrows
Arrowhead(x0 = c(root1, root2), y0 = -.25,
angle = 270,
npoint = 25,
arr.lwd = .2,
arr.length = .5,
arr.type = "curved",
arr.adj = 1,
lty = 1)
# coordinate point
points(x = c(root1, root2), y = c(3*root1, exp(root2)),
col = "white",
bg = "black",
pch = 21,
cex = .7)
# plot title
title(expression(e ^ x*phantom(" vs. 3*x")),
cex.main = 2,
col.main = "tomato")
title(expression(phantom(e ^ x)*" vs. "*phantom(3*x)),
cex.main = 1.5,
col.main = "black")
title(expression(phantom("e^x vs. ")*3*x),
cex.main = 2,
col.main = "dodgerblue")
1.7 Fixed-point iteration
A fixed point (also known as an invariant point) of a function is an element of the function’s domain that is mapped to itself by the function. That is to say, is a fixed point of the function \(f(x)\) if \(f(c) = c\). In graphical terms, a fixed point means the point \((x, f(x))\) is on the line \(y=x\). In other words, the graph of has a point in common with that line.
Ex: Solve \(f(x) = x^3 - x - 1 = 0\) on \((1, 2)\)
Notice \(f(1) = -1\) and \(f(2) = 5\), as a result the graph of \(f(x)\) crosses the \(x-axis\). In other words, a root of \(f(x)\) exists on \((1, 2).\)
- Rearrange the equation \(f(x)=0\) in the form \(x=g(x)\). \[ \begin{aligned} x^3 - x - 1 &= 0 \\ \\ x^3 &= x + 1 \text{ by adding x+1 to both sides.} \\ \\ x &= (x+1)^{\frac{1}{3}} \text{ by taking the cubic root of both sides.} \\ \\ g(x) &= (x+1)^{\frac{1}{3}}. \text{ substitute g(x) for x.} \\ \end{aligned} \]
- then an iterative method may be written as \(\quad x_{n+1} = g(x_n), n = 0, 1, 2, \dots\)
Plot of \(f(x) = \sqrt[3](x + 1)\)
# x continuum
x <- seq(-1, 1.5, by = .01)
root <- secant(f = function(x){x^(3) - x - 1}, -1.5, 1.5, epsilon = 1e-6)
# plot y = x
plot(x, x,
lwd = 1,
type = "l",
col = "black",
xlab = "",
ylab = "",
asp = 1,
axes = F,
xlim = c(-1.2, 1.5),
ylim = c(0, 1.5))
# plot function
lines(x, (x + 1)^(1/3),
type = "l",
lwd = 5,
ylab = "f(x)",
col = "tomato",
asp = 1)
# tracing lines
segments(x0 = root, y0 = root,
x1 = root, y1 = 0,
lty = 3,
lwd = .5)
# arrows
Arrowhead(x0 = root, y0 = 0,
angle = 270,
npoint = 25,
arr.lwd = .2,
arr.length = .5,
arr.type = "curved",
arr.adj = 1,
lty = 1)
# coordinate point
points(x = root, y = root,
col = "white",
bg = "tomato",
pch = 21,
cex = .7)
# plot title
title(main = expression(paste("g(x) = " * sqrt(x+1, 3))),
cex.main = 1.5,
col.main = "tomato")
# axes
axis(1, pos = 0)
axis(2, pos = 0)
axis(1, at = round(root, 2),
pos = 0,
col = "dodgerblue",
lwd.tick = 2,
col.axis = "dodgerblue",
font = 2)
Point of intersection of \(y =\sqrt[3](x+1)\) and \(y=x\).
plot.root(f = function(x){(x + 1)^(1/3) - x}, from = -1, to = 1.5, root = root, title = expression(paste("y = ", sqrt(x+1, 3), " - x")))
2 Numerical Integration
2.1 Plotting the Areas under Curves
INPUT: Function f; initial approximations \(a, b\); tolerance epsilon; number of shapes \(n\)
CONDITIONS: \(a < b\), \(f(a)\) and \(f(b)\) have opposite signs i.e. \(f(a) f(b) < 0\)
OUTPUT: Approximate solution of the definite integral of f.
plot.integral <-
function(f, a, b, from = a, to = b, title = NULL) {
x <- seq(from, to, length.out = 100)
y <- f(x)
# function plot
plot(x, y,
xlim = c(from, to),
ylim = c(ifelse(min(y) < 0, min(y), 0), max(y)),
xlab = "x",
ylab = "y",
main = title,
col.main = "tomato",
type = "l",
lwd = 3,
col = "tomato")
# area shading
x <- seq(a, b, length.out = 100)
y <- f(x)
polygon(c(x, b, a, a), c(y, 0, 0, f(a)),
border = adjustcolor("dodgerblue", alpha.f = 0.3),
col = adjustcolor("dodgerblue", alpha.f = 0.3))
} # end plot.integral
# save to the working directory getwd(), ls()
dump("plot.integral", file = "plot.integral.R")
2.1.1 Plot the area under \(e^{cos(x)}\) between \(x=0\) and \(x=\pi\)
plot.integral(f = function(x) exp(cos(x)),
a = 0, b = pi, from = -3*pi/2, to = 3*pi/2,
title = expression(paste("y = " * e^{cos(x)})))
2.2 Riemann Sums
Riemann summs is used for approximating the area of functions by dividing the region up into simpler shapes (rectangles) area of which could be easily computed. These shapes together form a region that is similar to the region being measured.
As described in (Kincaid and Cheney 2008, 185):
INPUT: Function \(f\), endpoint values \(a, b,\) tolerance epsilon, number of shapes \(n\)
CONDITIONS: \(a < b\), \(f(a)\) and \(f(b)\) have opposite signs i.e. \(f(a) f(b) < 0\)
OUTPUT: Approximate solution of the definite integral of f.
Riemann <-
function(f, a, b, epsilon = 1e-6, n = 1e5) {
h <- (b - a)/n
sum <- 0
for (i in n:1) {
x <- a + i*h
sum <- sum + f(x)
} # end for
sum.lower <- h*sum
sum.upper <- sum.lower + h*(f(a) - f(b))
return(list(sum.lower, sum.upper))
} # end Riemann
# save to the working directory getwd(), ls()
dump("Riemann", file = "Riemann.R")
2.2.1 Approximate the area under \(\displaystyle{\int\limits_{0}^{\pi} \mathrm{e}^{\mathrm{cos}(x)} \; dx}\) between \(0\) and \(\pi\).
Riemann(f = function(x) exp(cos(x)), a = 0, b = pi)## [[1]]
## [1] 3.977426
##
## [[2]]
## [1] 3.9775
2.2.2 Approximate the area under \(\displaystyle{\int\limits_{0}^{1}\mathrm{e}^{-x^2} \; dx}\) between \(0\) and \(1\).
Riemann(f = function(x) 1/exp(x^2), a = 0, b = 1)## [[1]]
## [1] 0.746821
##
## [[2]]
## [1] 0.7468273plot.integral(f = function(x) 1/exp(x^2),
a = 0, b = 1, from = -2, to = 2,
title = expression(paste("y = " * frac(1, e^{x^2}))))
2.2.3 The error function
\(\mathrm{erf}(x)\) is the error function encountered in integrating the normal distribution (which is a normalized form of the Gaussian function)
Approximate the area between \(0\) and \(1\) under the error function
\[ \mathrm{erf}(x) = \displaystyle{\frac{2}{\sqrt{\pi}}\int\limits_{0}^{x} \mathrm{e}^{-t^2} \; dt} \]
# erf(1) = 0.8427
Riemann(function(x) 2/(sqrt(pi)*exp(x^2)), a = 0, b = 1)## [[1]]
## [1] 0.8426972
##
## [[2]]
## [1] 0.8427044plot.integral(f = function(x) 2/(sqrt(pi)*exp(x^2)),
a = 0, b = 1, from = -2, to = 2,
title = expression(paste("y = " * frac(2, sqrt(pi*e^{x^2})))))
# compare
erf <- function(x) 2 * pnorm(x * sqrt(2)) - 1
erf(1)## [1] 0.8427008
2.3 Trapezoid Rule
Unlike the Riemann sums, the trapezoid method is based on an estimation of the area beneath a curve using trapezoids and considered an improvement in terms of reducing the error by fitting the function closer.
As described in (Kincaid and Cheney 2008, 191):
INPUT: Function \(f\), endpoint values \(a, b,\) tolerance epsilon, number of shapes \(n\)
CONDITIONS: \(a < b\), \(f(a)\) and \(f(b)\) have opposite signs i.e. \(f(a) f(b) < 0\)
OUTPUT: Approximate solution of the definite integral of f.
trapezoid <-
function(f, a, b, from = a, to = b, n = 1e4, plot = FALSE, title = NULL) {
source("plot.integral.R")
h <- (b - a)/n
sum <- .5*(f(a) + f(b))
for (i in 1:n) {
x <- a + i*h
sum <- sum + f(x)
} # end for
if (plot)
plot.integral(f, a, b, from, to, title)
return(sum*h)
} # end trapezoid
# save to the working directory getwd(), ls()
dump("trapezoid", file = "trapezoid.R")
2.3.1 Approximate the area under \(\displaystyle{\int\limits_{0}^{1} \mathrm{e}^{-x^2} \; dx}\) between \(0\) and \(1\).
trapezoid(function(x) exp(-x^2), a = 0, b = 1)## [1] 0.7468609plot.integral(function(x) exp(-x^2),
a = 0, b = 1, from = -2.5, to = 2.5,
title = expression(paste("y = ", (e^{-x^2}))))
2.3.2 Approximate the area under \(\displaystyle{\int\limits_{0}^{1} \frac{\sin(x)}{x} \; dx}\) between \(0\) and \(1\).
trapezoid(function(x) sin(x)/x, a = 0.000000001, b = 1)## [1] 0.9461672plot.integral(function(x) sin(x)/x,
a = 0.000000001, b = 1,
title = expression(paste("y = " * frac(sin(x), x))))
2.3.3 Arc-length of catenary by numerical integration
Taken from (Kincaid and Cheney 2008, 189):
Verify the numerical approximation given for the arc length in the catenary example above.
Solution:
From calculus, the length of a curve is
\[ \displaystyle{\int\limits_{a}^{b} \sqrt{1+ \left[ f'(x)\right]^2} \; dx} \]
where \(f\) is a function whose graph is the curve on the interval \(a \leq x \leq b\).
The arc length of a catenary is
\[ f(\theta) = \mathrm{cosh} \left(\frac{\theta}{\lambda} \right) \]
and its derivative is
\[ \frac{d}{d\theta} \; f(\theta) = \mathrm{sinh} \left(\frac{\theta}{\lambda} \right) \]
Substituting this into the arc-length formula we get
\[ \displaystyle{\int\limits_{a}^{b} \sqrt{1+ \left[ \mathrm{sinh} \left(\frac{\theta}{\lambda} \right)\right]^2} \; dx} \]
trapezoid(function(theta){sqrt(1 + (sinh(theta/lambda))^2)}, a = -50, b = 50)## [1] 102.6295
2.4 Romberg Algorithm
Romberg algorithm applies Richardson extrapolation to the trapezoid rule. Richardson extrapolation is a sequence acceleration method used to improve the rate of convergence of a sequence. The Romberg algorithm produces a triangular array of numbers, all of which are numerical estimates of a definite integral. The first column of this table contains estimates of the integral obtained by the recursive trapezoid formula with decreasing values of the step size. The successive columns in the Romberg array are generated by Richardson’s extrapolation formula:
\[ R(n, m) = R(n, m - 1) + \frac{1}{4^m - 1} \cdot \left[R(n, m - 1) - R(n - 1, m - 1)\right] \]
As described in (Kincaid and Cheney 2008, 206):
INPUT: Function \(f\), endpoint values \(a, b,\) tolerance epsilon, number of shapes \(n\)
CONDITIONS: \(a < b\), \(f(a)\) and \(f(b)\) have opposite signs i.e. \(f(a) f(b) < 0\)
OUTPUT: Approximate solution of the definite integral of f.
Romberg <-
function(f, a, b, n = 5) {
r <- matrix(NA, n, n)
h <- b - a
r[1, 1] <- (h/2)*(f(a) + f(b))
for (i in 2:n) {
h <- h/2
sum <- 0
for (k in seq(1, (2^(i-1)), by = 2)) {
sum <- sum + f(a + k*h)
} # end for k
r[i, 1] <- .5*r[i-1, 1] + sum*h
for (j in 2:i) {
r[i, j] <- r[i, j-1] + (r[i, j-1] - r[i-1, j-1])/(4^j - 1)
} # end for j
} # end for i
return(r)
} # end Romberg
# save to the working directory getwd(), ls()
dump("Romberg", file = "Romberg.R")
2.4.1 Approximate the area under \(\displaystyle{\int\limits_{0}^{1} \frac{4}{1+x^2} \; dx}\) between \(0\) and \(1\).
options(digits = 3)
Romberg(function(x) 4/(1 + x^2), a = 0, b = 1)## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.00 NA NA NA NA
## [2,] 3.10 3.11 NA NA NA
## [3,] 3.13 3.13 3.13 NA NA
## [4,] 3.14 3.14 3.14 3.14 NA
## [5,] 3.14 3.14 3.14 3.14 3.14
2.4.2 Number of prime integers by the logarithmic integral
Taken from (Kincaid and Cheney 2008, 189):
The logarithmic integral is a special mathematical function defined by the equation
\[ \mathrm{li}(x) = \displaystyle{\int\limits_{2}^{x} \frac{1}{\ln{t}} \; dt} \]
For large \(x\), the number of prime integers less than or equal to \(x\) is closely approximated by \(\mathrm{li}(x)\). For example, \(\mathrm{li}(200)\) is around 50 and there are 46 primes less than 200.
Verify that \(\mathrm{li}(200) \approx 50\).
Solution:
Romberg(f = function(x) 1/log(x), a = 2, b = 200, n = 7)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 161.5 NA NA NA NA NA NA
## [2,] 102.2 98.3 NA NA NA NA NA
## [3,] 73.5 71.6 71.2 NA NA NA NA
## [4,] 59.9 59.0 58.8 58.8 NA NA NA
## [5,] 53.6 53.2 53.1 53.1 53.1 NA NA
## [6,] 50.9 50.7 50.6 50.6 50.6 50.6 NA
## [7,] 49.7 49.7 49.6 49.6 49.6 49.6 49.6
2.5 Simpson’s Rule
Unlike the straight line segments used in the trapezoidal rule, Simpson’s rule approximates the integral of a function using quadratic polynomials (i.e. parabolic arcs).
As described in (Kincaid and Cheney 2008, 224):
INPUT: Function \(f\), endpoint values \(a, b,\) tolerance epsilon, number of shapes \(n\)
CONDITIONS: \(a < b\), \(f(a)\) and \(f(b)\) have opposite signs i.e. \(f(a) f(b) < 0\)
OUTPUT: Approximate solution of the definite integral of f.
Simpson <-
function(f, a, b, from = a, to = b, epsilon = 1e-6, level = 1, level.max = 101,
plot = FALSE, title = NULL) {
source("plot.integral.R")
level <- level + 1
h <- b - a
c <- (a + b)/2
one_simpson <- h*(f(a) + 4*f(c) + f(b))/6
d <- (a + c)/2
e <- (c + b)/2
two_simpson <- h*(f(a) + 4*f(d) + 2*f(c) + 4*f(e) + f(b))/12
if (level >= level.max) {
simpson_result <- two_simpson
cat("maximum level reached")
}
else if (abs(two_simpson - one_simpson) < 15*epsilon) {
simpson_result <- two_simpson + (two_simpson - one_simpson)/15
}
else {
left_simpson <- Simpson(f, a, c, epsilon/2, level, level.max)
right_simpson <- Simpson( f, c, b, epsilon/2, level, level.max)
simpson_result <- left_simpson + right_simpson
} # end if
if (plot)
plot.integral(f, a, b, from, to, title)
return(simpson_result)
} # end function Simpson
# save to the working directory getwd(), ls()
dump("Simpson", file = "Simpson.R")
2.5.1 Approximate the area under \(\displaystyle{\int\limits_{0}^{\frac{5}{4}\pi} \frac{\cos(2x)}{\mathrm{e}^x} \; dx}\) between \(0\) and \(\frac{5}{4}\pi\).
Simpson(function(x){cos(2*x)/exp(x)}, a = 0, b = 5*pi/4, from = -pi/4, to = 3*pi/2,
plot = TRUE, epsilon = 1e-10, level.max = 1e3,
title = expression(paste("y = ", frac(cos(2*x), e^{x}))))
## [1] 0.208
2.5.2 Normal distribution
Approximate the area under the standanrd normal curve
\[ \displaystyle{\int\limits_{-\infty}^{\infty} \frac{1}{\sigma \sqrt{2\pi}} \mathrm e^{-\frac{(x-\mu)^2}{2\sigma^2}} \; dx} \]
About 68% of values drawn from a normal distribution are within one standard deviation, \(\sigma\) away from the mean, \(\mu\).
options(digits = 10)
Simpson(function(x){exp(-.5*x^2)/sqrt(2*pi)}, a = -1, b = 1, from = -4, to = 4,
plot = TRUE, epsilon = 1e-10, level.max = 1e5,
title = expression(paste("y = ", frac(e^{-.5*x^2}, sqrt(2*pi))))) # 0.682689492137
## [1] 0.6826878338# about 0.954499736104 of the values lie within two standard deviations
Simpson(function(x){exp(-.5*x^2)/sqrt(2*pi)}, a = -2, b = 2, from = -4, to = 4,
plot = TRUE, title = expression(paste(paste("y = ", frac(1, sigma * sqrt(2 * pi)),
" ", e^{-frac((x - mu)^2, 2 * sigma^2)}))))
## [1] 0.9548815529# and about 0.997300203937 are within three standard deviations
Simpson(function(x){exp(-.5*x^2)/sqrt(2*pi)}, a = -3, b = 3, from = -4, to = 4,
plot = TRUE, title = expression(paste("y = ", frac(1, sigma * sqrt(2 * pi)),
" ", e^{-frac((x - mu)^2, 2 * sigma^2)})))
## [1] 1.001993116# 0.999936657516
Simpson(function(x){exp(-.5*x^2)/sqrt(2*pi)}, a = -4, b = 4, from = -4, to = 4,
plot = TRUE, title = expression(paste("y = ", frac(1, sigma * sqrt(2 * pi)),
" ", e^{-frac((x - mu)^2, 2 * sigma^2)})))
## [1] 1.006782818# 0.5 of the values lie below the mean
Simpson(function(x){exp(-.5*x^2)/sqrt(2*pi)}, a = -6, b = 0, from = -4, to = 4,
plot = TRUE, title = expression(paste("y = ", frac(1, sigma * sqrt(2 * pi)),
" ", e^{-frac((x - mu)^2, 2 * sigma^2)})))
## [1] 0.50241359323 Systems of Linear Equations
3.1 Naive Gaussian Elimination
Gaussian elimination is the standard method for solving an algebraic linear system Ax = b for the unknown vector x when the coefficient matrix A and right-hand side vector b are known.
As described in (Kincaid and Cheney 2008, 252):
INPUT: Coefficient matrix A and right-hand side vector b
OUTPUT: Unknown vector x
Naive.Gauss <-
function(A, b) {
a <- as.matrix(A)
n <- nrow(a)
x <- matrix(NA, n, dimnames = list(c(paste("x", 1:n, sep = "")), NULL))
for (k in 1:(n - 1)) {
for (i in (k + 1):n) {
xmult <- a[i,k]/a[k,k]
a[i,k] <- xmult
for (j in (k + 1):n) {
a[i,j] <- a[i,j] - (xmult)*a[k,j]
} #end for
b[i] <- b[i] - (xmult)*b[k]
} # end for
} # end for
x[n] <- b[n]/a[n,n]
for (i in (n - 1):1) {
sum <- b[i]
for (j in (i + 1):n) {
sum <- sum - a[i,j] * x[j]
} # end for
x[i] <- sum/a[i,i]
} #end for
return(x)
} # end Naive Gauss
# save to the working directory getwd(), ls()
dump("Naive.Gauss", file = "Naive.Gauss.R")
Solve the following system
\[ \begin{cases} 6x_1 - 2x_2 + 2x_3 + 4x_4 = 16 \\ 12x_1 - 8x_2 + 6x_3 + 10x_4 = 10 \\ 3x_1 - 13x_2 + 9x_3 + 3x_4 = -19 \\ -6x_1 + 4x_2 + x_3 - 18x_4 = -34 \\ \end{cases} \]
The coefficient and matrix and right-hand side vector:
# the coefficient matrix
A <- matrix(c(6, -2, 2, 4,
12, -8, 6, 10,
3, -13, 9, 3,
-6, 4, 1, -18),
4, 4, byrow = TRUE)
# right-hand side vector
b <- matrix(c(16, 26, -19, -34))Solution:
x <- Naive.Gauss(A, b); x## [,1]
## x1 3
## x2 1
## x3 -2
## x4 1# check
A %*% x # Ax = b## [,1]
## [1,] 16
## [2,] 26
## [3,] -19
## [4,] -34b## [,1]
## [1,] 16
## [2,] 26
## [3,] -19
## [4,] -34
3.1.1 Electrical Network Problem
Taken from (Kincaid and Cheney 2008, 245):
A simple electrical network contains a number of resistances and a single source of electromotive force (a battery) as shown below.
Electrical network
Using Kirchhoff’s laws and Ohm’s law, we can write a system of linear equations that govern this circuit. If x1, x2, x3, and x4 are the loop currents as shown, then the equations are
\[ \begin{matrix} 15x_1 & - & 2x_2 & - & 6x_3 & & & = & 300 \\ \\ -2x_1 & + & 12x_2 & - & 4x_3 & - & x_4 & = & 0 \\ \\ -6x_1 & - & 4x_2 & + & 19x_3 & - & 9x_4 & = & 0 \\ \\ & - & x_2 & - & 9x_3 & + & 21x_4 & = & 0 \\ \end{matrix} \]
Approximate the unknown vector x.
Solution:
# The coefficient matrix
A <- matrix(c(15, -2, -6, 0,
-2, 12, -4, -1,
-6, -4, 19, -9,
0, -1, -9, 21),
4, 4, byrow = TRUE)
# right-hand side vector
b <- matrix(c(300, 0, 0, 0))options(digits = 2)
Naive.Gauss(A, b)## [,1]
## x1 26.5
## x2 9.4
## x3 13.3
## x4 6.1
3.2 Gaussian Elimination with Scaled Partial Pivoting
As described in (Kincaid and Cheney 2008, 267):
INPUT: Coefficient matrix A and right-hand side vector b
OUTPUT: Unknown vector x
Gauss <-
function(A, b) {
a <- as.matrix(A)
b <- as.matrix(b)
# forward elimination with scaled partial pivoting
n <- nrow(a)
s <- matrix(NA, n)
l <- matrix(NA, n)
for (i in 1:n) {
l[i] <- i
smax <- 0
for (j in 1:n) {
smax <- max(smax, abs(a[i, j]))
} # end for
s[i] <- smax
} # end for
for (k in 1:(n - 1)) {
rmax <- 0
for (i in k:n) {
r <- abs(a[l[i], k] / s[l[i]])
if (r > rmax) {
rmax <- r
j <- i
} # end if
} # end for
temp <- l[j]
l[j] <- l[k]
l[k] <- temp
for (i in (k + 1):n) {
xmult <- a[l[i], k] / a[l[k], k]
a[l[i], k] <- xmult
for (j in (k + 1):n) {
a[l[i], j] <- a[l[i], j] - xmult*a[l[k], j]
} # end for
} # end for
} # end for
# The Back Substitution Phase
x <- matrix(NA, n, dimnames = list(c(paste("x", 1:n, sep = "")), NULL))
for (k in 1:(n - 1)) {
for (i in (k + 1):n) {
b[l[i]] <- b[l[i]] - a[l[i],k]*b[l[k]]
} # end for i
} # end for k
x[n] <- b[l[n]]/a[l[n],n]
for (i in (n - 1):1) {
sum <- b[l[i]]
for (j in (i + 1):n) {
sum <- sum - a[l[i],j]*x[j]
} # end for j
x[i] <- sum/a[l[i],i]
} #end for i
return(x)
} # end Gauss
# save to the working directory getwd(), ls()
dump("Gauss", file = "Gauss.R")
Solve the following system
\[ \begin{cases} 2x_1 + 4x_2 - 2x_3 &= \quad 6 \\ \:\:x_1 + 3x_2 + 4x_3 &= \:\:-1 \\ 5x_1 + 2x_2 &= \quad 2 \\ \end{cases} \]
Solution:
# The coefficient matrix
A <- matrix(c(2, 4, -2,
1, 3, 4,
5, 2, 0),
3, 3, byrow = TRUE)
# right-hand side vector
b <- matrix(c(6, -1, 2))Gauss(A, b)## [,1]
## x1 0
## x2 1
## x3 -1References
Burden, R. L., and J. D. Faires. 2011. Numerical Analysis. Brooks/Cole.
Kincaid, D. R., and E. W. Cheney. 2008. Numerical Mathematics and Computing. Cengage Learning.