This project originated in an Elementary Numerical Methods course at The Ohio State University in 2007, initially implemented in JAVA. The work was later translated to R to enhance its statistical applications. It focuses on fundamental numerical computation methods for solving linear systems, root-finding algorithms, and numerical integration.

The original text, code, and methodology have been refined and enhanced using AI tools to improve readability, efficiency, and accuracy.

1 Numerical Analysis

Numerical analysis is a field of mathematics and computer science that focuses on developing algorithms for obtaining accurate numerical approximations to problems involving continuous variables, particularly those with no tractable algebraic solution. The primary goal is to design and analyze techniques that provide efficient and reliable approximate solutions to complex problems.

1.1 Types of Numerical Methods

1.1.1 Direct Methods

These algorithms compute the solution in a finite number of steps and would yield the exact result if performed with infinite precision arithmetic.

1.1.2 Iterative Methods

Unlike direct methods, iterative approaches do not necessarily terminate in a finite number of steps. Instead:

  • They generate successive approximations from an initial guess.
  • They converge toward the exact solution in the limit.
  • A convergence criterion (often based on residual error) determines when the approximation is sufficiently accurate.

By balancing computational efficiency with precision, numerical analysis enables the practical solution of analytically intractable problems.

1.2 Locating Roots of Equations

A root-finding algorithm is a numerical method for identifying a value \(x\) such that \(f(x) = 0\) for a given function \(f\). Such an \(x\) is called a root (or zero) of \(f\).

Key Properties of Root-Finding Methods

  1. Iterative Approach: Most methods generate a sequence of approximations \({x_n}\) that converge to the root as \(n \to \infty\).

  2. Initial Guess: Requires a starting value \((x_0)\) near the suspected root.

  3. Approximate Solutions: The iteration terminates once the result meets a predefined tolerance, yielding an approximation rather than an exact solution.

1.2.1 Convergence Criteria

The algorithm stops when one (or a combination) of the following is satisfied:

  • Absolute Error: \(\left|x_{n+1} - x_n\right| < \epsilon\)

  • Relative Error: \(\displaystyle{\frac{\left|x_{n+1} - x_n\right|}{\left|x_n\right|} < \epsilon}\)

  • Residual Threshold: \(\left|f(x_n)\right| < \epsilon\)

where \(\epsilon\) is a user-defined tolerance.

Common Challenges

  • Convergence Dependence: Speed and success often depend on the choice of \(x_0\).
  • Multiple Roots: Some roots may be missed if the initial guess is poorly chosen.
  • Non-Convergence: Not all functions guarantee convergence for arbitrary \(x_0\).

1.3 Plotting the Roots of Functions

#' Visualize a Function and Its Roots
#'
#' Plots a function over a specified interval and marks its roots on the x-axis.
#'
#' @param f The function to be plotted (must accept vector inputs).
#' @param a Lower bound of the plotting range.
#' @param b Upper bound of the plotting range.
#' @param root Optional vector of root locations to mark.
#' @param title Plot title (optional).
#'
#' @return None (produces a plot)
#' @export
#'
#' @examples
#' f <- function(x) x^2 - 4
#' plot_root(f, a = -3, b = 3, root = c(-2, 2), title = "Quadratic Function")

plot_root <- function(f, a, b, root = NULL, title = NULL) {
  if (!is.function(f)) stop("f must be a function")
  if (length(a) != 1 || length(b) != 1) stop("Range must be numeric values")
  if (a >= b) stop("'a' must be less than 'b'")
  
  x <- seq(a, b, length.out = 100)
  y <- f(x)
  
  y_range <- range(y, na.rm = TRUE)
  y_pad <- max(0.1 * diff(y_range), 1e-6)
  y_limits <- c(y_range[1] - y_pad, y_range[2] + y_pad)
  
  plot(x, y,
       type = "l",
       lwd = 3,
       col = "tomato",
       xlim = c(a, b),
       ylim = y_limits,
       xlab = "x",
       ylab = "f(x)",
       main = title,
       col.main = "tomato",
       axes = FALSE,
       panel.first = {
         grid(col = "gray90")
         axis(1, pos = 0)
         axis(2, pos = 0)
       })
  
  if (!is.null(root) && length(root) > 0) {
    points(root, rep(0, length(root)),
           col = "white",
           bg = "tomato",
           pch = 21,
           cex = 1)
    
    axis(1, at = round(root, 2), 
         pos = (ceiling(max(f(x))) - floor(min(f(x)))) / 10,
         tick = FALSE,
         col.axis = "darkred")
  }
}


# save to the working directory getwd(), ls()
dump("plot_root", file = "plot_root.R")

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}, 
          a = -2, b = 2, root = c(-1, 1), 
          title = expression(paste("y = " * x^2)))

1.3.1 The Bisection Method

The bisection method (also known as the interval halving method, binary search method, or dichotomy method) is a robust root-finding algorithm that systematically narrows down an interval containing a root of a continuous function.

Key Properties

  • Bracketing Method: Requires two initial points \(a\) and \(b\) such that \(f(a)\) and \(f(b)\) have opposite signs (i.e., \(f(a) \cdot f(b) < 0\)), guaranteeing a root exists in \([a,b]\) by the Intermediate Value Theorem.

  • Linear Convergence: The error decreases by approximately half each iteration, with convergence rate \(O(1/2^n)\).

  • Absolute Error Bound: After \(n\) iterations, the maximum error is \(\frac{b-a}{2^n}\).

Algorithm Steps

  1. Initialization: Choose \(a_0\), \(b_0\) with \(f(a_0)f(b_0) < 0\)

  2. Iteration: For each step \(n\):

  • Compute midpoint \(c_n = \frac{a_n + b_n}{2}\)
  • Evaluate \(f(c_n)\)
  • Update interval:
  • If \(f(a_n)f(c_n) < 0\), set \(b_{n+1} = c_n\)
  • Else, set \(a_{n+1} = c_n\)
  1. Termination: Stop when \(\left|b_n - a_n\right| < \epsilon\) (desired tolerance)

1.3.1.1 Advantages & Limitations

Advantages Disadvantages
✓ Guaranteed convergence for continuous functions ✗ Slow convergence (linear rate)
✓ Simple to implement ✗ Requires bracketing interval
✓ Error bound known exactly ✗ Doesn’t generalize to higher dimensions

1.3.1.2 Practical Considerations

  • Often used to generate initial guesses for faster methods (e.g., Newton-Raphson)

  • Optimal for functions where evaluations are computationally expensive

  • Convergence criteria typically combine absolute and relative error thresholds: \[ \text{Stop if } \left|b_n - a_n\right| < \epsilon_{\text{abs}}+\epsilon_{\text{rel}}\left|c_n\right| \]

As described in (Kincaid & Cheney, 2008, p. 76):

#' Bisection Method for Root Finding (Modified for Plotting)
#'
#' Finds a root of a continuous function within a given interval, with optional plotting.
#'
#' @param f The function to find roots for.
#' @param a Lower bound of the interval.
#' @param b Upper bound of the interval.
#' @param e Tolerance for convergence (default = 1e-8).
#' @param max_iter Maximum number of iterations (default = 1000).
#' @param plot Whether to generate a plot (default = FALSE).
#' @param title Plot title (optional).
#' @return The approximate root.
#' @export
#'
#' @examples
#' bisection(cos, 0, 2, plot = TRUE, title = "cos(x)")

bisection <- function(f, a, b, tol = 1e-8, max_iter = 1000, plot = FALSE, title = NULL, ...) {
  
  if (f(a) * f(b) >= 0) {
    stop("Function must have opposite signs at endpoints (f(a) * f(b) < 0)")
  }
  
  iterations <- 0
  root <- NULL
  converged <- FALSE
  
  while (iterations < max_iter) {
    iterations <- iterations + 1
    c <- (a + b) / 2
    
    fc <- f(c)
    
    if (abs(fc) < tol || (b - a)/2 < tol) {
      root <- c
      converged <- TRUE
      break
    }
    
    if (f(a) * fc < 0) {
      b <- c
    } else {
      a <- c
    }
  }
  
  if (plot && converged) {
    plot_root(f, a, b, root, title, ...)
  }
  
  list(root = root, value = if (!is.null(root)) f(root) else NULL, iterations = iterations, converged = converged)
}

# save to the working directory getwd(), ls()
dump("bisection", file = "bisection.R")

Example Applications:

  1. The root of \(\mathrm{cosine}\) function on \(\left[\frac{-\pi}{3}, \frac{4\pi}{3}\right]\)
# Example usage with cos(x)
bisection(f = function(x){return(cos(x))}, a = -pi/3, b = 4*pi/3, plot = TRUE, title = expression(paste("y = " * cos(x))))

## $root
## [1] 1.570796
## 
## $value
## [1] 6.123234e-17
## 
## $iterations
## [1] 1
## 
## $converged
## [1] TRUE
  1. Roots of Quadratic Functions and the Bisection Method

A general quadratic function: \[ f(x) = ax^2 + bx + c \quad (a \neq 0) \]

has roots given by the quadratic formula: \[ x = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a} \]

Three possible cases:

  1. Two distinct real roots: When discriminant \(D = b^2 - 4ac > 0\)

  2. One real double root: When \(D = 0\) (vertex touches x-axis)

  3. Complex conjugate roots: When \(D < 0\)

Bisection Method Behavior with Quadratics

Key Characteristics:

  • Single-Root Finding:

  • The bisection method can only find one root per execution

  • Which root is found depends entirely on the initial bracketing interval \([a,b]\)

  • Interval Selection Rules:

  • For two real roots at \(r_1 < r_2\):

  • To find \(r_1\): Choose \(a < r_1\) and \(b \in (r_1,r_2)\)

  • To find \(r_2\): Choose \(a \in (r_1,r_2)\) and \(b > r_2\)

  • Cannot find both roots simultaneously

  • Special Cases:

  • Double root: Will fail to find root unless \(a\) or \(b\) equals the root (as \(f(x)\) doesn’t change sign)

  • No real roots: Method will error out (as no bracketing interval exists)

Practical Example: Consider \(f(x) = x^2 - 3x + 2\) with roots at \(x=1\) and \(x=2\):

f <- function(x) {x^2 - 3*x + 2}
# Find left root (x=1):
bisection(f, 0, 1.5)  # Returns ≈1
## $root
## [1] 1
## 
## $value
## [1] -7.450581e-09
## 
## $iterations
## [1] 26
## 
## $converged
## [1] TRUE
# Find right root (x=2):
bisection(f, 1.5, 3)  # Returns ≈2
## $root
## [1] 2
## 
## $value
## [1] -7.450581e-09
## 
## $iterations
## [1] 26
## 
## $converged
## [1] TRUE
# Define the quadratic function
f <- function(x) x^2 - 3*x + 2

# Find both roots using bisection
root1 <- bisection(f, 0, 1.5, e = 1e-8)$root
root2 <- bisection(f, 1.5, 3, e = 1e-8)$root

# Plot the function with both roots marked
plot_root(f, 
          a = -1,  # Adjusted for better visualization
          b = 4,     # Extended to see full behavior
          root = c(root1, root2),  # Using found roots, not hardcoded values
          title = "Quadratic Function f(x) = x² - 3x + 2")

# Verify the roots
cat(sprintf("Found roots at x = %.6f and x = %.6f\n", root1, root2))
## Found roots at x = 1.000000 and x = 2.000000
  1. Root of logarithmic functions: \(\ln\left(\frac{1 + x}{1 - x^2}\right)\)
f <- function(x) log((1 + x)/(1 - x^2))
bisection(f, a = -0.99, b = 0.99, tol = 1e-10, plot = F)
## $root
## [1] 0
## 
## $value
## [1] 0
## 
## $iterations
## [1] 1
## 
## $converged
## [1] TRUE

4. Board in hall problem

Taken from (Kincaid & Cheney, 2008, p. 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
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)}, a = 0, b = pi/4)

gamma
## $root
## [1] 0.4511959
## 
## $value
## [1] 5.061549e-07
## 
## $iterations
## [1] 27
## 
## $converged
## [1] TRUE

then substitute it in the length equations, \(\mathscr{l}_1\) and \(\mathscr{l}_2\)

# l2
(l2 <- 7 * 1/sin(gamma$root))
## [1] 16.0535
# l1
(l1 = 9 * 1/sin(pi - (125*pi/180) - gamma$root))
## [1] 18.47772
# total length l
(l = l1 + l2)
## [1] 34.53122

2 Method of False Position (Regula Falsi)

The False Position Method, also known as Regula Falsi, is a root-finding algorithm that improves upon the bisection method. Instead of halving the interval, it uses linear interpolation (a secant line) between the endpoints to estimate the root. This often leads to faster convergence while still ensuring the root remains bracketed between the endpoints.

Mathematical Formulation

Given a continuous function $ f(x) $ with $ f(a) $ and $ f(b) $ having opposite signs, the next estimate $ c $ is computed as:

\[ c = b - \frac{f(b)(b - a)}{f(b) - f(a)} \] This formula corresponds to the x-intercept of the secant line connecting the points \((a,f(a))\) and \((b,f(b))\).

R Implementation

As described in (Burden & Faires, 2011, p. 73):

#' Regula Falsi (False Position) Method for Root Finding
#'
#' Finds a root of a continuous function within a specified interval using the false position method.
#' Uses plot_root() for visualization if plot=TRUE.
#'
#' @param f The function for which the root is sought (must be continuous)
#' @param a Lower bound of the interval
#' @param b Upper bound of the interval
#' @param e Tolerance for convergence (default = 1e-6)
#' @param NMAX Maximum number of iterations (default = 101)
#' @param plot Logical indicating whether to plot the results using plot_root (default = FALSE)
#' @param title Optional title for the plot
#'
#' @return The approximate root of the function in the given interval
#' @export
#'
#' @examples
#' f <- function(x) x^2 - 4
#' regula.falsi(f, from = 1, to = 3, plot = TRUE, title = "Finding sqrt(4)")

regula.falsi <- function(f, a, b, e = 1e-6, NMAX = 101, plot = FALSE, title = NULL) {
  # Input validation
  if (a >= b) stop("'a' must be less than 'b'")
  if (f(a) * f(b) >= 0) stop("Function must have opposite signs at endpoints")
  
  from <- a
  to <- b
  
  # Initialize variables
  fa <- f(a)
  fb <- f(b)
  root_history <- numeric(NMAX)
  
  # Main iteration loop
  for (i in 1:NMAX) {
    # Calculate new point using false position formula
    p <- b - fb * (b - a) / (fb - fa)
    root_history[i] <- p
    fp <- f(p)
    
    # Check for convergence
    if (abs(fp) < e) {
      if (plot) {
        # Use plot_root function for visualization
        plot_root(f, from, to, root = p, title = title)
        
        # Add iteration points in a different color
        points(root_history[1:i], rep(0, i), 
               col = "darkorange", pch = 1, lwd = 1)
      }
      return(p)
    }
    
    # Update interval
    if (fp * fb < 0) {
      a <- b
      fa <- fb
    }
    b <- p
    fb <- fp
  }
}

# Save function to file
dump("regula.falsi", file = "regula.falsi.R")

The root of \(\mathrm{sine}\) function on \(\left[\frac{\pi}{16}, \; \frac{31\pi}{16}\right]\)

regula.falsi(f = function(x){log(x)}, a = pi/16, b = 31*pi/16, plot = TRUE,
             title  = expression(paste("y = " * sin(x))))

## [1] 1.000001

2.1 Newton-Raphson Method

The Newton-Raphson method is an efficient iterative algorithm for finding successively better approximations to the roots of a real-valued function. It uses first-order derivative information to achieve quadratic convergence near the root.

Given a function $ f(x)$ and its derivative $ f’(x) $, the method starts with an initial guess $x_0 $ and iteratively improves the estimate using:

\[ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \]

The iteration continues until either:

  • \(|f(x_n)| < \epsilon\) (function value is sufficiently small), or
  • \(|x_{n+1} - x_n| < \epsilon\) (change between iterations is negligible)
#' Newton-Raphson Method (Simplified Robust Version)
#' Finds a root of a function using Newton-Raphson iteration.
#'
#' @param f Function whose root is to be found
#' @param f_prime Derivative of the function
#' @param x0 Initial guess (default = 1)
#' @param tol Tolerance for convergence (default = 1e-8)
#' @param max_iter Maximum number of iterations (default = 100)
#'
#' @return A list: root (or NA), number of iterations, and convergence status
#' @export
#'
#' @examples
#' f <- function(x) ((x - 2)*x + 1)*x - 3
#' newton_raphson(f, x0 = 1)

newton_raphson <- function(f, x0 = 1, tol = 1e-8, max_iter = 100) {
  # Install and load Deriv only if missing
  if (!requireNamespace("Deriv", quietly = TRUE)) {
    install.packages("Deriv")
  }
  library(Deriv)
  
  # Precompute derivative function once
  f_prime <- Deriv(f)
  
  x <- x0
  iterations <- 0
  converged <- FALSE
  
  for (i in 1:max_iter) {
    iterations <- i
    fx <- tryCatch(f(x), error = function(e) NA)
    dfx <- tryCatch(f_prime(x), error = function(e) NA)
    
    if (any(is.na(c(fx, dfx)))) {
      warning("Function or derivative evaluation failed at x = ", x)
      break
    }
    
    if (abs(dfx) < .Machine$double.eps) {
      warning("Derivative too small at x = ", x)
      break
    }
    
    x_new <- x - fx / dfx
    
    if (abs(fx) < tol || abs(x_new - x) < tol) {
      converged <- TRUE
      x <- x_new
      break
    }
    
    x <- x_new
  }
  
  list(
    root = if (converged) x else NA,
    iterations = iterations,
    converged = converged
  )
}


# save to the working directory getwd(), ls()
dump("newton_raphson", file = "newton_raphson.R")

 

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 \]

# Example usage with your specific function
f <- function(x) ((x - 2)*x + 1)*x - 3
f_prime <- function(x) (3*x - 4)*x + 1

# Find root 
newton_raphson(f, x0 = 2)
## $root
## [1] 2.174559
## 
## $iterations
## [1] 5
## 
## $converged
## [1] TRUE

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)

(root1 <- newton_raphson(f, x0 = 1))
## $root
## [1] 0.7071068
## 
## $iterations
## [1] 5
## 
## $converged
## [1] TRUE
(root2 <- newton_raphson(f, x0 = -1))
## $root
## [1] -0.7071068
## 
## $iterations
## [1] 5
## 
## $converged
## [1] TRUE
# Extract just the root values
plot_root(
  function(r) r^2 - 1/2,
  a = -1.5,
  b = 1.5,
  root = c(root1$root, root2$root),
  title = expression(paste("y = " * r^2 - frac(1,2)))
)

2.2 Secant Method

The Secant Method is a numerical technique for finding the roots of a real-valued function. It is a variation of the Newton-Raphson method but avoids explicit differentiation. Instead of requiring the derivative, the Secant Method approximates it using the slope of a secant line through two prior points:

\[ x_{n+1} = x_n - f(x_n) \cdot \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})} \] This method is especially useful when:

  • The derivative \(f′(x)\) is difficult or impossible to compute.

  • An analytical derivative is unavailable, but the function can be evaluated.

Compared to Newton-Raphson:

  • No derivative needed

  • May converge more slowly

  • May not converge if the function behaves poorly near the root

As described in (Kincaid & Cheney, 2008, p. 122):

#' Secant Method for Root Finding
#'
#' Finds a root of a continuous function within a given interval using the secant method,
#' which approximates the derivative using finite differences.
#'
#' @param f A continuous function for which the root is sought.
#' @param a Initial first guess for the root.
#' @param b Initial second guess for the root.
#' @param nmax Maximum number of iterations (default = 100).
#' @param tol Tolerance for convergence (default = 1e-6).
#' @param plot Logical. If `TRUE`, plots the function and root (default = `FALSE`).
#' @param title Optional title for the plot.
#'
#' @return The approximate root if convergence is achieved; otherwise, `NA`.
#' @export
#'
#' @examples
#' f <- function(x) x^3 - x - 2
#' root <- secant(f, a = 1, b = 2, plot = TRUE, title = "Root of x^3 - x - 2")
#' print(root)
secant <- function(f, a, b, nmax = 100, tol = 1e-6, plot = FALSE, title = NULL) {
  source("plot_root.R")  # Load plotting function if needed
  
  from <- a
  to <- b
  
  fa <- f(a)
  fb <- f(b)
  
  # Ensure 'a' is the better initial guess (smaller |f(x)|)
  if (abs(fa) > abs(fb)) {
    temp <- a; a <- b; b <- temp
    temp <- fa; fa <- fb; fb <- temp
  }
  
  for (n in 1:nmax) {
    # Secant update: x_{n+1} = x_n - f(x_n) * (x_n - x_{n-1}) / (f(x_n) - f(x_{n-1}))
    dx <- (b - a) / (fb - fa) * fa
    a_new <- a - dx
    
    # Check for convergence
    if (abs(dx) < tol) {
      if (plot) {
        plot_root(f, from, to, root = a_new, title = title)
      }
      return(a_new)
    }
    
    # Update points for next iteration
    b <- a
    fb <- fa
    a <- a_new
    fa <- f(a)
  }
  
  warning("Maximum iterations reached without convergence.")
  return(NA)  # Return NA if no convergence
}
# 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}, a = -1.5, b = 1.5, plot = TRUE, title = expression(paste("y = " * x^3 - x - 1)))

## [1] 1.324718

The arc-length of catenary by root finding

Taken from (Kincaid & Cheney, 2008, p. 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
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}, a = 100, b = 200))
## [1] 126.6324

Substitute \(\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.6187

2.3 Points of intersection of two graphs of functions

Taken from (Kincaid & Cheney, 2008, p. 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, a = 0, b = 1))
## [1] 0.6190616
(root2 <- regula.falsi(f, a = 1, b = 2))
## [1] 1.512134
plot_root(f, a = 0, b = 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
if (!require("shape")) 
  install.packages("shape", dependencies = TRUE)
library(shape)

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")

2.4 Fixed-point iteration

Fixed-point iteration is a fundamental numerical method used to find solutions to equations of the form \(x=g(x)\). These solutions, called fixed points, are points where the function gg intersects the line \(y=x\). This method is particularly useful for solving root-finding problems \(f(x)=0\) by reformulating them into fixed-point form.

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).\)

  1. 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} \]

Solving \(x^3 - x - 1 = 0\)

Reformulate as \(x = g(x)\): \[ g(x) = \sqrt[3]{x + 1} \]

Check \(|g'(x)| < 1\): \[ g'(x) = \frac{1}{3}(x + 1)^{-2/3} \]

For \(x \in [1, 2]\), \(|g'(x)| \approx 0.2 < 1 \rightarrow\) Convergence guaranteed.

  1. then an iterative method may be written as \(\quad x_{n+1} = g(x_n), n = 0, 1, 2, \dots\)
#' Fixed-Point Iteration Method
#'
#' Finds a fixed point of `g(x)`, equivalent to solving `x = g(x)`.
#'
#' @param g The function for which a fixed point is sought.
#' @param x0 Initial guess (default = 0).
#' @param tol Tolerance for convergence (default = 1e-6).
#' @param max_iter Maximum iterations (default = 100).
#' @param plot Logical. If `TRUE`, plots convergence (requires `ggplot2`).
#'
#' @return Approximate fixed point if converged; otherwise `NA`.
#' @export
#'
#' @examples
#' # Solve x^3 - x - 1 = 0 via x = (x + 1)^(1/3)
#' g <- function(x) (x + 1)^(1/3)
#' root <- fixed_point(g, x0 = 1.5, plot = TRUE)
#' print(root)  # Approx. 1.324717
#' 
fixed_point <- function(g, x0 = 0, tol = 1e-6, max_iter = 100, plot = FALSE) {
  x <- x0
  history <- data.frame(iter = 0, x = x0)
  
  for (iter in 1:max_iter) {
    x_new <- g(x)
    history <- rbind(history, data.frame(iter = iter, x = x_new))
    
    if (abs(x_new - x) < tol) {
      if (plot) {
        if (!requireNamespace("ggplot2", quietly = TRUE)) {
          warning("`ggplot2` not found. Plot disabled.")
        } else {
          p <- ggplot2::ggplot(history, ggplot2::aes(x = iter, y = x)) +
            ggplot2::geom_line(color = "blue") +
            ggplot2::geom_point(color = "red") +
            ggplot2::labs(
              title = "Fixed-Point Iteration Convergence",
              x = "Iteration",
              y = "Estimate"
            )
          print(p)
        }
      }
      return(x_new)
    }
    x <- x_new
  }
  
  warning("Did not converge in ", max_iter, " iterations.")
  return(NA)
}

Plot of \(f(x) = \sqrt[3](x + 1)\)

g <- function(x) (x + 1)^(1/3)
root <- fixed_point(g, x0 = 1.5, plot = TRUE)

print(root)  # Output: ~1.324717 (the real root)
## [1] 1.324718
# x continuum
x <- seq(-1, 1.5, by = .01)

root <- secant(f = function(x){x^(3) - x - 1}, a = -1.5, b = 1.5)

# 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}, a = -1, b = 1.5, root = root, title = expression(paste("y = ", sqrt(x+1, 3), " - x")))

3 Numerical Integration

Plotting the Areas under Curves

#' Plot a Function and Shade the Integral Area
#'
#' This function plots a continuous function \( f \) over a specified range and
#' shades the area under the curve between points \( a \) and \( b \).
#'
#' @param f A function of one variable to plot.
#' @param a Lower limit of integration (start of shaded area).
#' @param b Upper limit of integration (end of shaded area).
#' @param from Numeric, start of the x-axis plotting range. Defaults to \code{a}.
#' @param to Numeric, end of the x-axis plotting range. Defaults to \code{b}.
#' @param title Optional character string for the plot title.
#'
#' @return No return value. The function produces a plot as a side effect.
#'
#' @examples
#' f <- function(x) sin(x)
#' plot_integral(f, 0, pi, from = 0, to = 2*pi, title = "Integral of sin(x)")
#'
#' @export
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_area <- seq(a, b, length.out = 100)
  y_area <- f(x_area) 
  polygon(c(x_area, b, a, a), c(y_area, 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")

3.0.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)})))

3.1 Riemann Sums

Riemann sums are a fundamental technique to approximate the area under a curve \(y=f(x)\) over an interval \([a,b]\). The idea is to partition the interval into \(n\) subintervals, then approximate the area under the curve by summing the areas of simple shapes—most commonly rectangles—over each subinterval.

Mathematical Definition

Given a partition

\[ a = x_0 < x_1 < \cdots < x_n = b, \]

the width of the \(i-th\) subinterval is

\[ \Delta x_i = x_i - x_{i-1}. \]

We choose a sample point \(c_i \in [x_{i-1}, x_i]\) and approximate the area over that subinterval by the rectangle with height \(f(c_i)\):

\[ S_n = \sum_{i=1}^n f(c_i) \Delta x_i. \]

As \[ \to \infty\] and the partition gets finer, \[S_n\] approaches the definite integral

\[ \int_a^b f(x) \, dx. \]

Types of Riemann Sums

  • Left Riemann sum: Use \(c_i = x_{i-1}\) (left endpoint of each subinterval)

  • Right Riemann sum: Use \(c_i = x_i\) (right endpoint)

  • Midpoint Riemann sum: Use

\[ c_i = \frac{x_{i-1} + x_i}{2} \]

As described in (Kincaid & Cheney, 2008, p. 185):

#' Approximate definite integral using the Trapezoidal Rule
#'
#' @param f Function to integrate
#' @param a Lower limit of integration
#' @param b Upper limit of integration
#' @param n Number of subintervals (default 1e5)
#' @return Numeric approximation of the integral
#' @export
trapezoidal <- function(f, a, b, n = 1e5) {
  h <- (b - a) / n
  x <- seq(a, b, length.out = n + 1)
  y <- f(x)
  
  # Trapezoidal sum formula
  integral <- (h / 2) * (y[1] + 2 * sum(y[2:n]) + y[n + 1])
  
  return(integral)
}

# Save to working directory
dump("trapezoidal", file = "trapezoidal.R")

Approximate the area under \(\displaystyle{\int\limits_{0}^{\pi} \mathrm{e}^{\mathrm{cos}(x)} \; dx}\) between \(0\) and \(\pi\).

trapezoidal(f = function(x) exp(cos(x)), a = 0, b = pi)
## [1] 3.977463
#' Plot function and trapezoidal approximation
#'
#' @param f Function to plot and approximate
#' @param a Lower bound of interval
#' @param b Upper bound of interval
#' @param n Number of subintervals (default 10 for clearer visualization)
#' @param main Optional plot title
#' @export
plot_trapezoidal <- function(f, a, b, n = 10, main = NULL) {
  x <- seq(a, b, length.out = 1000)
  y <- f(x)
  
  # Plot the function curve
  plot(x, y, type = "l", lwd = 2, col = "blue",
       xlab = "x", ylab = "f(x)", main = main)
  
  # Partition points for trapezoids
  xs <- seq(a, b, length.out = n + 1)
  ys <- f(xs)
  
  # Draw trapezoids
  for (i in 1:n) {
    polygon(
      x = c(xs[i], xs[i], xs[i + 1], xs[i + 1]),
      y = c(0, ys[i], ys[i + 1], 0),
      col = adjustcolor("orange", alpha.f = 0.5),
      border = "red"
    )
  }
  
  # Overlay points on partition
  points(xs, ys, pch = 19, col = "red")
}

# Save to working directory
dump("plot_trapezoidal", file = "plot_trapezoidal.R")
plot_trapezoidal(function(x) exp(cos(x)), 0, pi, n = 8, main = "Trapezoidal Approximation")

Approximate the area under \(\displaystyle{\int\limits_{0}^{1}\mathrm{e}^{-x^2} \; dx}\) between \(0\) and \(1\).

trapezoidal(f = function(x) 1/exp(x^2), a = 0, b = 1)
## [1] 0.7468241
plot_trapezoidal(f = function(x) 1/exp(x^2), 
                 a = 0, b = 1, n = 4,
                 main = expression(paste("y = " * frac(1, e^{x^2}))))

3.1.1 The Error Function

The error function, denoted by $ (x) $, arises in probability, statistics, and partial differential equations. It is defined as:

\[ \mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} \, dt \] This function represents the probability that a normally distributed random variable with mean 0 and variance $ $ falls between 0 and $ x $.

In statistics, $ (x) $ is used in computing the cumulative distribution function (CDF) of the standard normal distribution, which is a normalized form of the Gaussian function.

# erf(1) = 0.8427
trapezoidal(function(x) 2/(sqrt(pi)*exp(x^2)), a = 0, b = 1)
## [1] 0.8427008
plot_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})))))

plot_trapezoidal(f = function(x) 2/(sqrt(pi)*exp(x^2)), 
                 a = 0, b = 1, n = 4,
                 main = 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

Arc-length of catenary by numerical integration

Taken from (Kincaid & Cheney, 2008, p. 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} \]

trapezoidal(function(theta){sqrt(1 + (sinh(theta/lambda))^2)}, a = -50, b = 50)
## [1] 102.6187

3.2 Romberg Algorithm

The Romberg algorithm applies Richardson extrapolation to the trapezoidal rule to improve the accuracy of numerical integration. Richardson extrapolation is a technique for accelerating the convergence of a sequence, allowing for more accurate approximations.

Romberg integration constructs a triangular array of estimates for a definite integral. The method proceeds as follows:

  • The first column consists of estimates using the composite trapezoidal rule, where the step size is halved with each row:

\[R(k,1)=T(hk),\] where \(hk=b−a2k−1\)

  • Each subsequent column is computed using Richardson extrapolation, which eliminates lower-order error terms:

\[ R[k, j] = R[k, j-1] + (1 / (4^(j-1) - 1)) * (R[k, j-1] - R[k-1, j-1]) \] Where:

  • \(R(k,j)\) is the Romberg estimate at row kk, column jj,

  • \(T(hk)\) is the trapezoidal approximation with step size \(hk\),

The array is built iteratively until a desired level of convergence is achieved or a maximum depth is reached.

Advantages:

  • Converges more rapidly than the trapezoidal or Simpson’s rule for sufficiently smooth functions.

  • Produces highly accurate integral approximations using fewer function evaluations.

\[ 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 & Cheney, 2008, p. 206):

#' Romberg Integration
#'
#' Numerically approximates the definite integral of a function using the Romberg method.
#' This algorithm applies Richardson extrapolation to improve accuracy over the trapezoidal rule.
#'
#' @param f Function to integrate. Should be a function of one variable (e.g., \code{function(x) sin(x)}).
#' @param a Lower limit of integration.
#' @param b Upper limit of integration.
#' @param n Number of rows in the Romberg integration table (default = 5).
#'
#' @return A matrix containing the Romberg table of estimates. The most accurate estimate
#'         is in the bottom-right corner \code{[n, n]}.
#'
#' @examples
#' f <- function(x) exp(-x^2)
#' Romberg(f, 0, 1)
#'
#' # Compare with built-in integration
#' integrate(f, 0, 1)
#'
#' @export
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)
    }
    
    r[i, 1] <- 0.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)
    }
  }
  
  return(r)
}

# Save to file for later sourcing
dump("Romberg", file = "Romberg.R")

 

3.2.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

Number of prime integers by the logarithmic integral

Taken from (Kincaid & Cheney, 2008, p. 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

 

3.3 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—that is, parabolic arcs. By fitting a parabola through every three consecutive points on the function, Simpson’s Rule achieves a higher degree of accuracy for smooth functions over an interval.

As described in (Kincaid & Cheney, 2008, p. 224):

#' Adaptive Simpson's Rule for Numerical Integration
#'
#' Approximates the definite integral of a function \code{f} over the interval \code{[a, b]}
#' using adaptive Simpson's Rule. Optionally plots the function and shaded area.
#'
#' @param f A function to integrate.
#' @param a Lower bound of integration.
#' @param b Upper bound of integration.
#' @param from Plot range lower bound (default is \code{a}).
#' @param to Plot range upper bound (default is \code{b}).
#' @param epsilon Desired accuracy (default is \code{1e-6}).
#' @param level Internal recursion level (should not be manually set by user).
#' @param level.max Maximum recursion depth to avoid infinite recursion (default is 101).
#' @param plot Logical; if \code{TRUE}, plots the function with shaded integral region (default is \code{FALSE}).
#' @param title Optional title for the plot.
#'
#' @return The approximate value of the definite integral.
#' @examples
#' Simpson(function(x) sin(x), 0, pi, plot = TRUE, title = "∫ sin(x) dx from 0 to π")
#'
#' @export
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
    }
    
    if (plot)
      plot.integral(f, a, b, from, to, title)
    
    return(simpson_result)
  }

# Save to working directory
dump("Simpson", file = "Simpson.R")

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
#' Plot Simpson's Rule Approximation with Parabolic Arcs
#'
#' @param f Function to approximate integral for
#' @param a Lower bound of integration
#' @param b Upper bound of integration
#' @param n Number of subintervals (must be even)
#' @param title Optional plot title
#'
#' @return NULL (plots graph)
plot_simpson_bins <- function(f, a, b, n = 6, title = NULL) {
  if (n %% 2 != 0) {
    stop("n must be even for Simpson's rule.")
  }
  
  x <- seq(a, b, length.out = n + 1)
  y <- f(x)
  
  plot(x, y, type = "l", lwd = 2, col = "blue", 
       xlab = "x", ylab = "f(x)", main = title)
  
  for (i in seq(1, n, by = 2)) {
    # Get three points per parabolic segment
    xs <- x[i:(i+2)]
    ys <- y[i:(i+2)]
    
    # Fit a quadratic polynomial through these three points
    coef <- coef(lm(ys ~ poly(xs, 2, raw = TRUE)))
    
    # Create fine points for smooth parabola
    x_fine <- seq(xs[1], xs[3], length.out = 100)
    y_fine <- coef[1] + coef[2]*x_fine + coef[3]*x_fine^2
    
    # Draw the parabola (arc)
    lines(x_fine, y_fine, col = "red", lwd = 3)
    
    # Shade area under parabola segment
    polygon(c(x_fine, rev(x_fine)), c(y_fine, rep(0, length(y_fine))),
            col = adjustcolor("red", alpha.f = 0.3), border = NA)
  }
}
f <- function(x){cos(2*x)/exp(x)}
plot_simpson_bins(f, 0, 5*pi/4, n = 6, title = "Simpson's Rule Approximation with Parabolic Arcs")

3.3.1 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
f <- function(x){exp(-.5*x^2)/sqrt(2*pi)}
plot_simpson_bins(f, -1, 1, n = 4, title = "Simpson's Rule Approximation with Parabolic Arcs")

# 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.5024135932

4 Systems of Linear Equations

4.1 Naive Gaussian Elimination

Gaussian elimination is the standard method for solving an algebraic linear system $ A = $ for the unknown vector \(\mathbf{x}\) when the coefficient matrix $ A $ and right-hand side vector $ $ are known.

As described in (Kincaid & Cheney, 2008, p. 252):

#' Naive Gaussian Elimination for Solving Linear Systems
#'
#' Solves the system \( A \mathbf{x} = \mathbf{b} \) using the Naive Gaussian Elimination method.
#'
#' @param A Numeric matrix representing the coefficient matrix \( A \).
#' @param b Numeric vector representing the right-hand side vector \( \mathbf{b} \).
#' @return Numeric vector \(\mathbf{x}\), the solution to the system.
#' @examples
#' A <- matrix(c(2,1,-1, -3,-1,2, -2,1,2), nrow=3, byrow=TRUE)
#' b <- c(8, -11, -3)
#' Naive.Gauss(A, b)
Naive.Gauss <- function(A, b) {
  a <- as.matrix(A)
  n <- nrow(a)
  
  x <- numeric(n)
  names(x) <- paste0("x", 1:n)
  
  # Forward elimination
  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]
      }
      
      b[i] <- b[i] - xmult * b[k]
    }
  }
  
  # Back substitution
  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]
    }
    x[i] <- sum / a[i, i]
  }
  
  return(x)
}

# Save the function to the working directory
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))
## x1 x2 x3 x4 
##  3  1 -2  1
# Verify solution from Naive Gaussian Elimination

# Compute Ax to check if it matches b
Ax <- A %*% x

# Print computed Ax
print("Computed Ax:")
## [1] "Computed Ax:"
print(Ax)
##      [,1]
## [1,]   16
## [2,]   26
## [3,]  -19
## [4,]  -34
# Print original right-hand side vector b
print("Original vector b:")
## [1] "Original vector b:"
print(b)
##      [,1]
## [1,]   16
## [2,]   26
## [3,]  -19
## [4,]  -34
# Explanation:
# If the Naive Gaussian Elimination solved the system correctly,
# then Ax should be (approximately) equal to b.
# Minor differences may occur due to floating-point arithmetic.

Electrical Network Problem

Taken from (Kincaid & Cheney, 2008, p. 245):

A simple electrical network contains a number of resistances and a single source of electromotive force (a battery) as shown below.

Electrical network
Electrical network

Using Kirchhoff’s laws and Ohm’s law, we can write a system of linear equations that govern this circuit. If \(x_1\), \(x_2\), \(x_3\), and \(x_4\) 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)
##   x1   x2   x3   x4 
## 26.5  9.4 13.3  6.1

4.2 Gaussian Elimination with Scaled Partial Pivoting

As described in (Kincaid & Cheney, 2008, p. 267):

## Gaussian Elimination with Scaled Partial Pivoting
#' Gaussian Elimination with Scaled Partial Pivoting
#'
#' Solves the linear system \(Ax = b\) using Gaussian elimination with
#' scaled partial pivoting to improve numerical stability.
#'
#' @param A Numeric matrix representing the coefficient matrix.
#' @param b Numeric vector or matrix representing the right-hand side vector.
#'
#' @return Numeric vector containing the solution vector \(x\).
#'
#' @details
#' The function performs forward elimination with scaled partial pivoting,
#' followed by back substitution to find the solution vector.
#' Scaled partial pivoting selects pivot rows based on the ratio of the
#' largest element in each row to the pivot element to minimize rounding errors.
#'
#' @examples
#' A <- matrix(c(2,1,1,4,3,3,8,7,9), 3, 3, byrow = TRUE)
#' b <- c(2,5,10)
#' Gauss(A, b)
#'
#' @export
Gauss <-
  function(A, b) {
    a <- as.matrix(A)
    b <- as.matrix(b)
    
    # forward elimination with scaled partial pivoting
    n <- nrow(a)
    s <- numeric(n)
    l <- 1:n
    
    for (i in 1:n) { 
      s[i] <- max(abs(a[i, ]))
    }
    
    for (k in 1:(n - 1)) { 
      rmax <- 0
      j <- k
      
      for (i in k:n) {
        r <- abs(a[l[i], k]) / s[l[i]]
        if (r > rmax) { 
          rmax <- r
          j <- i 
        }
      }
      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]
        }
      }
    }
    
    # Back substitution phase
    x <- numeric(n)
    
    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]]
      }
    }
    
    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]
      }
      
      x[i] <- sum / a[l[i], i]
    }
    
    names(x) <- paste0("x", 1:n)
    return(x)
  }

# save to the working directory getwd(), ls()
dump("Gauss", file = "Gauss.R")

Solve the following system

\[ \begin{cases} 2x_1 + 4x_2 - 2x_3 = \:\:\:6 \\ \:\:x_1 + 3x_2 + 4x_3 = -1 \\ 5x_1 + 2x_2 \qquad \: = \:\:\: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)
## x1 x2 x3 
##  0  1 -1

References

Burden, R. L., & Faires, J. D. (2011). Numerical analysis. Brooks/Cole.
Kincaid, D. R., & Cheney, E. W. (2008). Numerical mathematics and computing. Cengage Learning.