Note: If you wish to replicate the R code below, then you will need to copy and paste the following commands in R first (to make sure you have all the packages and their dependencies):

install.packages(c("install.load", "import", "ramify", "ggplot2", "ggpmisc", "gsubfn", "sessioninfo", "pracma", "units"))
# install the packages and their dependencies, including the extra system dependencies (this process may take a while depending on the number of dependencies)



# load the required packages and provide the session information

install.load::load_package("ramify", "ggplot2", "ggpmisc", "gsubfn", "units")
# load needed packages using the load_package function from the install.load
# package (it is assumed that you have already installed these packages)


import::from(pracma, golden_ratio, fminbnd)
# import golden_ratio, fminbnd from the pracma package

sessioninfo::session_info()$packages
##  package      * version date       source        
##  assertthat     0.2.0   2017-04-11 CRAN (R 3.5.1)
##  backports      1.1.2   2017-12-13 CRAN (R 3.5.1)
##  base64enc      0.1-3   2015-07-28 CRAN (R 3.5.1)
##  bindr          0.1.1   2018-03-13 CRAN (R 3.5.1)
##  bindrcpp       0.2.2   2018-03-29 CRAN (R 3.5.1)
##  clisymbols     1.2.0   2017-05-21 CRAN (R 3.5.1)
##  colorspace     1.3-2   2016-12-14 CRAN (R 3.5.1)
##  crayon         1.3.4   2017-09-16 CRAN (R 3.5.1)
##  digest         0.6.17  2018-09-12 CRAN (R 3.5.1)
##  dplyr          0.7.6   2018-06-29 CRAN (R 3.5.1)
##  evaluate       0.11    2018-07-17 CRAN (R 3.5.1)
##  formatR        1.5     2017-04-25 CRAN (R 3.5.1)
##  ggplot2      * 3.0.0   2018-07-03 CRAN (R 3.5.1)
##  ggpmisc      * 0.3.0   2018-07-16 CRAN (R 3.5.1)
##  glue           1.3.0   2018-07-17 CRAN (R 3.5.1)
##  gsubfn       * 0.7     2018-03-16 CRAN (R 3.5.1)
##  gtable         0.2.0   2016-02-26 CRAN (R 3.5.1)
##  htmltools      0.3.6   2017-04-28 CRAN (R 3.5.1)
##  import         1.1.0   2015-06-22 CRAN (R 3.5.1)
##  install.load   1.2.1   2016-07-12 CRAN (R 3.5.1)
##  knitr          1.20    2018-02-20 CRAN (R 3.5.1)
##  lazyeval       0.2.1   2017-10-29 CRAN (R 3.5.1)
##  magrittr       1.5     2014-11-22 CRAN (R 3.5.1)
##  munsell        0.5.0   2018-06-12 CRAN (R 3.5.1)
##  pillar         1.3.0   2018-07-14 CRAN (R 3.5.1)
##  pkgconfig      2.0.2   2018-08-16 CRAN (R 3.5.1)
##  plyr           1.8.4   2016-06-08 CRAN (R 3.5.1)
##  pracma         2.1.5   2018-08-25 CRAN (R 3.5.1)
##  proto        * 1.0.0   2016-10-29 CRAN (R 3.5.1)
##  purrr          0.2.5   2018-05-29 CRAN (R 3.5.1)
##  R6             2.2.2   2017-06-17 CRAN (R 3.5.1)
##  ramify       * 0.3.3   2016-12-17 CRAN (R 3.5.1)
##  Rcpp           0.12.18 2018-07-23 CRAN (R 3.5.1)
##  rlang          0.2.2   2018-08-16 CRAN (R 3.5.1)
##  rmarkdown    * 1.10    2018-06-11 CRAN (R 3.5.1)
##  rprojroot      1.3-2   2018-01-03 CRAN (R 3.5.1)
##  scales         1.0.0   2018-08-09 CRAN (R 3.5.1)
##  sessioninfo    1.0.0   2017-06-21 CRAN (R 3.5.1)
##  stringi        1.2.4   2018-07-20 CRAN (R 3.5.1)
##  stringr        1.3.1   2018-05-10 CRAN (R 3.5.1)
##  tibble         1.4.2   2018-01-22 CRAN (R 3.5.1)
##  tidyselect     0.2.4   2018-02-26 CRAN (R 3.5.1)
##  units        * 0.6-0   2018-06-09 CRAN (R 3.5.1)
##  withr          2.1.2   2018-03-15 CRAN (R 3.5.1)
##  yaml           2.2.0   2018-07-25 CRAN (R 3.5.1)


Example 1 (modified from Source 1 and Wickham)



# Source 1 begins

# option 1
integrate_vector <- Vectorize(integrate, vectorize.args = c("f", "lower", "upper"))

funcs <- list(function(x) x^5 - 7 * x^4 + 9 * x^2 - 3 * x + 1, function(x) 8 * 
    sin(x) + 3 * cos(x), function(x) 8 * exp(x)/x^4)

integrate_vector(funcs, lower = c(0, -4 * pi, 5), upper = c(100, 4 * pi, 30))
##              [,1]         [,2]         [,3]      
## value        152669651767 -5.58059e-15 122584887 
## abs.error    0.001694974  3.323246e-08 128.3941  
## subdivisions 1            2            1         
## message      "OK"         "OK"         "OK"      
## call         Expression   Expression   Expression
# option 2
f1 <- c(f = function(x) x^5 - 7 * x^4 + 9 * x^2 - 3 * x + 1, l = 0, u = 100)

f2 <- c(f = function(x) 8 * sin(x) + 3 * cos(x), l = -4 * pi, u = 4 * pi)

f3 <- c(f = function(x) 8 * exp(x)/x^4, l = 5, u = 30)

funs <- list(f1, f2, f3)

lapply(funs, function(x) integrate(x[["f"]], x[["l"]], x[["u"]]))
## [[1]]
## 152669651767 with absolute error < 0.0017
## 
## [[2]]
## -5.58059e-15 with absolute error < 3.3e-08
## 
## [[3]]
## 122584887 with absolute error < 128
# option 3
functions <- list(function(x) x^5 - 7 * x^4 + 9 * x^2 - 3 * x + 1, function(x) 8 * 
    sin(x) + 3 * cos(x), function(x) 8 * exp(x)/x^4)

lower <- c(0, -4 * pi, 5)

upper <- c(100, 4 * pi, 30)

mapply(integrate, functions, lower, upper)
##              [,1]         [,2]         [,3]      
## value        152669651767 -5.58059e-15 122584887 
## abs.error    0.001694974  3.323246e-08 128.3941  
## subdivisions 1            2            1         
## message      "OK"         "OK"         "OK"      
## call         Expression   Expression   Expression
# Source 1 ends



Example 2

“Suppose you have measured temperatures at a number of coordinates on the surface of a rectangular heated plate”:

“T(2, 1) = 60
T(9, 1) = 57.5
T(2, 6) = 55
T(9, 6) = 70”

“Use bilinear interpolation to estimate the temperature at xi = 5.25 and yi = 4.8.” (Chapra 381)



library("ramify")

import::from(pracma, interp2)
# import interp2 from the pracma package


x <- mat("2, 9")
x
##      [,1] [,2]
## [1,]    2    9
y <- mat("1, 6")
y
##      [,1] [,2]
## [1,]    1    6
z <- mat("60, 57.5; 55, 70")
z
##      [,1] [,2]
## [1,]   60 57.5
## [2,]   55 70.0
T <- interp2(x, y, z, 5.25, 4.8)

T
## [1] 61.21429

The “temperature at xi = 5.25 and yi = 4.8 is 61.”



Using GNU Octave to obtain the answer for Chapra 381:



x = [2 9]

y = [1 6]

z = [60 57.5; 55 70]

interp2(x, y, z, 5.25, 4.8)
## GNU Octave, version 3.8.1
## Copyright (C) 2014 John W. Eaton and others.
## This is free software; see the source code for copying conditions.
## There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE.  For details, type 'warranty'.
## 
## Octave was configured for "x86_64-pc-linux-gnu".
## 
## Additional information about Octave is available at http://www.octave.org.
## 
## Please contribute if you find this software useful.
## For more information, visit http://www.octave.org/get-involved.html
## 
## Read http://www.octave.org/bugs.html to learn how to submit bug reports.
## For information about changes from previous versions, type 'news'.
## 
## warning: function /usr/share/octave/packages/strings-1.1.0/strjoin.m shadows a core library function
## x =
## 
##    2   9
## 
## y =
## 
##    1   6
## 
## z =
## 
##    60.000   57.500
##    55.000   70.000
## 
## ans =  61.214



Example 3

“The deflection of a uniform beam subject to a linearly increasing distributed load can be computed as

y = w_0 / 120EIL * (-x^5 + 2L2x3 - L^4x)

Given that L = 600 cm, E = 50,000 kN/cm^2, I = 30,000 cm^4, and w_0 = 2.5 kN/cm, determine the point of maximum deflection (a) graphically, (b) using the golden-section search until the approximate error falls below epsilon_s = 1% with initial guesses of x_l = 0 and x_u = L." (Chapra 184) / Problem 7.16



install.load::load_package("ggplot2", "ggpmisc", "gsubfn")
# load needed packages using the load_package function from the install.load
# package (it is assumed that you have already installed these packages)

import::from(pracma, golden_ratio)
# import golden_ratio from the pracma package



y <- function(x) {
    y1 <- (w0/(120 * E * I * L)) * (-x^5 + 2 * L^2 * x^3 - L^4 * x)
    return(y1)
}

L <- 600  # cm

E <- 50000  # kN/cm^2

I <- 30000  # cm^2

w0 <- 2.5  # kN/cm

epsilons <- 1/100  # %

xl <- 0

xu <- L

list[xmin, fmin, iter, ea] <- golden_ratio(y, xl, xu, tol = epsilons)

xmin
## [1] 268.3303
fmin
## [1] -0.5151901
iter
## [1] 20
ea
## [1] 0.009363442
# use ggplot2 to graphically determine the maximum deflection
ggplot(data.frame(x = xl:xu, y = y(xl:xu)), aes(x, y)) + geom_line() + labs(x = "Length (cm)", 
    y = "Deflection (cm)", title = "Beam Deflection") + stat_valleys(colour = "purple") + 
    stat_valleys(y.label.fmt = "%1.2f", geom = "text", colour = "red", vjust = 1.5, 
        aes(label = paste(..y.label.., "cm at", ..x.label.., "cm"))) + theme_bw()  # Source 2



Using GNU Octave to obtain the answer for Chapra 184 / Problem 7.16:



L = 600 % cm

E = 50000 % kN/cm^2

I = 30000 % cm^2

w0 = 2.5 % kN/cm

epsilons = 1 / 100 % %

xl = 0

xu = L


z = @(x) (w0 / (120 * E * I * L)) * (-x ^ 5 + 2 * L ^ 2 * x ^ 3 - L ^ 4 * x)


function [x, k] = golden_section( f, a, b, tol, maxiter, verbose )
  %---------------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % Written May 3, 1999
  % Revised May 4, 2005
  % Revised December 18, 2008 - Works with both Matlab and Octave
  %
  % Implements the golden section search as described on pages 555-556 of
  % Numerical Mathematics and Computing, 4th Edition, by Cheney & Kincaid,
  % Brooks/Cole, 1999 except that y is calculated as an offset from b rather
  % than from a.
  %---------------------------------------------------------------------------

  if ( nargin <= 5 )
    verbose = 0;
  end

  % Set things up to begin the iteration
  r  = ( sqrt( 5 ) - 1 ) / 2;   % Golden ratio
  x = a + r * ( b - a );
  y = b - r * ( b - a );
  u = f( x );
  v = f( y );

  % Main loop
  k = 0;
  if ( verbose ~= 0 )
    fprintf( '%4d: Interval is [%f, %f]; y = %f, x = %f\n', k, a, b, y, x );
  end
  while ( abs( b - a ) >= 2 * tol && k < maxiter )
    k = k + 1;
    if ( u > v )        % minimum must be in [a, x]
      b = x;
      x = y;
      u = v;
      y = b - r * ( b - a );
      v = f( y );
    else                % minimum must be in [y, b]
      a = y;
      y = x;
      v = u;
      x = a + r * ( b - a );
      u = f( x );
    end
    if ( verbose ~= 0 )
      fprintf( '%4d: Interval is [%f, %f]; y = %f, x = %f\n', k, a, b, y, x );
    end
  end

  % All done; best estimate of solution is midpoint of final interval
  x = 0.5 * ( a + b );

endfunction % needed for GNU Octave only



[xmin, iter] = golden_section(z, 0, 8, tol = 1e-08, maxiter = 100)


xmin

iter
## GNU Octave, version 3.8.1
## Copyright (C) 2014 John W. Eaton and others.
## This is free software; see the source code for copying conditions.
## There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE.  For details, type 'warranty'.
## 
## Octave was configured for "x86_64-pc-linux-gnu".
## 
## Additional information about Octave is available at http://www.octave.org.
## 
## Please contribute if you find this software useful.
## For more information, visit http://www.octave.org/get-involved.html
## 
## Read http://www.octave.org/bugs.html to learn how to submit bug reports.
## For information about changes from previous versions, type 'news'.
## 
## warning: function /usr/share/octave/packages/strings-1.1.0/strjoin.m shadows a core library function
## L =  600
## E =  50000
## I =  30000
## w0 =  2.5000
## epsilons =  0.010000
## xl = 0
## xu =  600
## z =
## 
## @(x) (w0 / (120 * E * I * L)) * (-x ^ 5 + 2 * L ^ 2 * x ^ 3 - L ^ 4 * x)
## 
## xmin =  8.0000
## iter =  42
## xmin =  8.0000
## iter =  42



Example 4

“The Streeter-Phelps model can be used to compute the dissolved oxygen concentration in a river below a point discharge of sewage,

where o = dissolved oxygen concentration (mg/L), o_s = oxygen saturation concentration (mg/L), t = travel time (d), L_o = biochemical oxygen demand (BOD) concentration at the mixing point (mg/L), k_d = rate of decomposition of BOD (d^-1), k_s = rate of settling of BOD (d^-1), k_a = reaeration rate (d^-1), and S_b = sediment oxygen demand (mg/L/d)."

The given equation “produces an oxygen”sag" that reaches a critial minimum level o_c, some travel time t_c below the point discharge. This point is called “critical” because it represents the location where biota that depend on oxygen (like fish) would be the most stressed. Determine the critical travel time and concentration, given the following values:

o_s = 10 mg/L k_d = 0.1 d^-1 k_o = 0.6 d^-1

k_s = 0.05 d^-1 L_o = 50 mg/L S_b = 1 mg/L/d"

(Chapra 185) / Problem 7.25



import::from(pracma, fminbnd)
# import fminbnd from the pracma package


library("gsubfn")

o <- function(t) {
    o1 <- os - ((kd * L0/(kd + ks - ka)) * (exp(-ka * t) - exp((-kd + ks) * 
        t))) - Sb/ka * (1 - exp(-ka * t))
    return(o1)
}

os <- 10  # mg/L

kd <- 0.1  # d^-1

ka <- 0.6  # d^-1

ks <- 0.05  # d^-1

L0 <- 50  # mg/L

Sb <- 1  # mg/L/d

list[oc, tc] <- fminbnd(o, 0, 5)

oc
## [1] 4.772125
tc
## [1] 0.3102083



Using GNU Octave to obtain the answer for Chapra 184 / Problem 7.16:



function o1 = o(t)
global os kd L0 ks ka Sb
o1 = os - ((kd * L0 / (kd + ks - ka)) * (exp(-ka * t) - exp((-kd + ks) * t))) - Sb / ka * (1 - exp(-ka * t));
endfunction % only needed for GNU Octave


global os kd L0 ks ka Sb
os = 10 % mg/L

kd = 0.1 % d^-1

ka = 0.6 % d^-1

ks = 0.05 % d^-1

L0 = 50 % mg/L

Sb = 1 % mg/L/d

[oc, tc] = fminbnd(@o, 0, 5)

oc

tc
## GNU Octave, version 3.8.1
## Copyright (C) 2014 John W. Eaton and others.
## This is free software; see the source code for copying conditions.
## There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE.  For details, type 'warranty'.
## 
## Octave was configured for "x86_64-pc-linux-gnu".
## 
## Additional information about Octave is available at http://www.octave.org.
## 
## Please contribute if you find this software useful.
## For more information, visit http://www.octave.org/get-involved.html
## 
## Read http://www.octave.org/bugs.html to learn how to submit bug reports.
## For information about changes from previous versions, type 'news'.
## 
## warning: function /usr/share/octave/packages/strings-1.1.0/strjoin.m shadows a core library function
## os =  10
## kd =  0.10000
## ka =  0.60000
## ks =  0.050000
## L0 =  50
## Sb =  1
## oc =  4.7721
## tc =  0.31021
## oc =  4.7721
## tc =  0.31021



Example 5

“Strength of Materials - Elastic Deformation”

“A 1 in (25 mm) diameter soft steel rod carries a tensile load of 15,000 lbf (67 kN). The elongation is 0.158 in (4 mm). The modulus of elasticity is \(2.9 x 10^7\) lbf/in^2 (200 GPa). What is the total length of the rod?” (Lindeburg 44-1) / Problem 44-1



library("units")

mm <- as_units("mm")

m <- as_units("m")

delta <- 4/1000 * mm
delta
## 0.004 mm
kN <- with(ud_units, 1000 * kg * m/s^2)

GPa <- with(ud_units, 1 * 10^9 * kg/m * s^2)

F <- 67 * kN
F
## 67000 kg*m/s^2
E <- 200 * GPa
E
## 2e+11 kg*s^2/m
d <- 25/1000 * mm
d
## 0.025 mm
A <- (pi * d^2)/4
A
## 0.0004908739 mm^2
delta <- as.numeric(delta)

F <- as.numeric(F)

E <- as.numeric(E)

d <- as.numeric(d)

A <- as.numeric(A)


Lofun <- function(Lo) {
    delta - (Lo * F)/(E * A)
}

Louse <- uniroot(Lofun, interval = c(1e-07, 200), extendInt = "yes")

Lo <- Louse$root


L <- Lo + delta

L <- L * m

L
## 5.86518 m

The total length (L) of the rod is 5.8651803."



Example 6

“Numerical Analysis - Newton’s method”

“Newton’s method is being used to find the roots of the equation f(x) = (x - 2) ^ 2 - 1. What is the third approximation of the root if 9.33 is chosen as the first approximation?” () / Problem 15



library("animation")

newton.method(FUN = function(x) (x - 2)^2 - 1, init = 9.33, rg = c(-1, 10), 
    tol = 0.001, interact = FALSE, col.lp = c("blue", "red", "green"))