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("install.load")

# you will need to have GNU Octave (check the RcppOctave CRAN page for more details) installed on your system, then you can use the code below:

# install and/or load the packages and their dependencies, including the extra system dependencies (this process may take a while depending on the number of dependencies)

install.packages("import")
# install the import package

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

# 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 through RcppOctave to obtain the answer for Chapra 381:

library(RcppOctave)

o_source(text = "
x = [2 9]

y = [1 6]

z = [60 57.5; 55 70]

interp2(x, y, z, 5.25, 4.8)
")
## x =
##
##    2   9
##
## y =
##
##    1   6
##
## z =
##
##    60.000   57.500
##    55.000   70.000
##
## ans =  61.214

# Example 3

“Determine the point of maximum deflection of a beam using the golden-section search.” (Chapra 184) / Problem 7.16

install.load::load_package("ggplot2", "ggpmisc")
# package (it is assumed that you have already installed these packages)

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

source("https://raw.githubusercontent.com/ggrothendieck/gsubfn/master/R/list.R")

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 create the figure
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 through RcppOctave to obtain the answer for Chapra 184 / Problem 7.16:

library(RcppOctave)

o_source(text = "
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
")

# Example 4

“Determine the critical travel time and concentration” using the Streeter-Phelps model. (Chapra 185) / Problem 7.25

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

source("https://raw.githubusercontent.com/ggrothendieck/gsubfn/master/R/list.R")

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 through RcppOctave to obtain the answer for Chapra 184 / Problem 7.16:

library(RcppOctave)

o_source(text = "
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
")
## 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

## Sources used in the R code

Source 1
anonymous function exercise in Wickham’s Advanced R book (integration of 3 functions in list) - Stack Overflow answered by jlhoward on Dec 19 2014 & answered and edited by Khashaa on Dec 19 2014. See https://stackoverflow.com/questions/27572408/anonymous-function-exercise-in-wickhams-advanced-r-book.

Source 2
Using R for photobiology: R tips 3: ggpmisc adds new stats to ‘ggplot2’ by Pedro J. Aphalo, 2016-01-30. See http://www.r4photobiology.info/2016/01/ggpmisc-adds-new-stats-to-ggplot2/.