R Code for Climate Mathematics:
Theory and Applications
A Cambridge University Press book By
SSP Shen and RCJ Somerville
Version 1.0 released in July 2019 and coded by Dr. Samuel Shen, Distinguished Professor
San Diego State University, California, USA
Email: sshen@sdsu.edu
Version 2.0 compiled and updated by Momtaza Sayd
San Diego State University May 2021
Version 3.0 compiled and updated by Joaquin Stawsky
San Diego State University June 2022
Chapter 6: Calculus Applications to Climate Science I: Derivatives
Plot Fig. 6.1
#setEPS()
#postscript("fig0601.eps", height = 8, width = 10)
#Cess' Budyko parameters
par(mar = c(4.5,5,2,0.5))
A <- 212
B <- 1.6
ep <- 0.6
sg <- 5.670373*10^(-8)
T <- seq(-50,50, by=0.1)
S <- 1365
alf <- 0.30
plot(T,A+B*T, xlim = c(-50,50), ylim = c(0,500),
type = "l",lwd = 4, cex.lab = 1.2, cex.axis = 1.2,cex.main = 1.35,
main = "Radiation Law and Its Approximations",
xlab = expression(paste("Temperature T [",degree,"C]")),
ylab = expression(paste("Outgoing Radiation: ", E[out], " [W/",m^2,"]")))
lines(T, ep*sg*(273.15+T)^4,col = "blue", lwd = 3)
lines(T, 189.4 + 2.77*T,type = "l", col = "purple",lwd = 3)
lines(T, (1-alf)*S/4*rep(1, length(T)), col = "red",lwd = 2)
legend(-50, 510, legend = c("Budyko radiation formula",
"Stefan-Boltzmann radiation law",
"Linear Approximation of SB law",
"Observed outgoing radiation"),
col = c("black", "blue","purple","red"),
lty = 1, bty = "n",lwd = 3, cex = 1.2)
text(-40,180, "Budyko formula",cex = 1.2)
text(15,310, "Stefan-Boltzmann law", col = "blue",cex = 1.2)
text(-3,65, "Linear Approximation to Stefan-Boltzmann law",
col = "purple",cex = 1.2)
text(-37,255, expression(paste( E[out], " = 239 [W/",m^2,"]")),
col = "red", cex = 1.2)
#dev.off()
Generate Table 6.1
x <- seq(-0.3,0.3, by = 0.1)
fx <- c(1:7)
lx <- c(1:7)
ex <- c(1:7)
for(i in 1:7){fx[i] <- (1+x[i])^4
lx[i] <- 1+4*x[i]
ex[i] <- ((1+x[i])^4-(1+4*x[i]))/((1+x[i])^4)*100
}
round(cbind(x, fx, lx, ex), digits = 4)## x fx lx ex
## [1,] -0.3 0.2401 -0.2 183.2986
## [2,] -0.2 0.4096 0.2 51.1719
## [3,] -0.1 0.6561 0.6 8.5505
## [4,] 0.0 1.0000 1.0 0.0000
## [5,] 0.1 1.4641 1.4 4.3781
## [6,] 0.2 2.0736 1.8 13.1944
## [7,] 0.3 2.8561 2.2 22.9719
Solve an equation by the bisection method
#install.packages("NLRoot") #if not installed before
library(NLRoot)
func <- function(x){x^3-x-1}
#[1] 1.324716 is the solution
BFfzero(func, 0, 2)## [1] 1
## [1] 1.324716
## [1] -7.551464e-06
## [1] "finding root is successful"
Another bisection example
f1 <- function(x){(1+x)^4-(2+x)}
#[1] 0.2207428 is the solution
BFfzero(f1,0,2)## [1] 1
## [1] 0.2207428
## [1] -8.076542e-06
## [1] "finding root is successful"
Solve an EBM by the bisection method
S <- 1368
ep <- 0.6
sg <- 5.670373*10^(-8)
f2 <- function(T){ep*sg*T^4 - (1-(0.5 - 0.2*tanh((T-265)/10)))*(S/4)}
#[1] 289.3097 is the solution
BFfzero(f2,270,300)## [1] 1
## [1] 289.3097
## [1] -7.509889e-06
## [1] "finding root is successful"
R code for the Newton’s method to solve an equation
newton <- function(f, tol = 1E-12,x0 = 1,N = 20) {
h <- 0.001
i <- 1; x1 <- x0
p <- numeric(N)
while (i <= N) {
df.dx <- (f(x0+h)-f(x0))/h
x1 <- (x0 - (f(x0)/df.dx))
p[i] <- x1
i <- i + 1
if (abs(x1-x0) < tol) break
x0 <- x1
}
return(p[1:(i-1)])
}
f <- function(x) { x^3 + 4*x^2 -10 }
root <- newton(f, tol = 1E-12, x0 = 1, N = 10)
root## [1] 1.454256 1.368917 1.365238 1.365230 1.365230 1.365230 1.365230uniroot(f,c(1,2), tol = 10^(-6))## $root
## [1] 1.36523
##
## $f.root
## [1] 6.71676e-11
##
## $iter
## [1] 6
##
## $init.it
## [1] NA
##
## $estim.prec
## [1] 5e-07
Plot Fig. 6.2
#Illustration of Newton's method
#setEPS()
#postscript("fig0602.eps", height = 6, width = 8)
par(mar = c(0,1,0,0))
x <- seq(0.2,1.7, len = 30)
f <- function(x) { x^3 + 4*x^2 -10 }
g <- function(x) { 3*x^2 + 8*x }
plot(x,f(x), type = 'l', lwd = 1.5, bty = "n",
xaxt = "n",yaxt = "n", xlab = "", ylab = "",
xlim = c(0,1.8),ylim = c(-10,7))
axis(1, at = c(0, 0.5, 1.0, 1.5, 1.7), pos = 0, cex.lab = 1.3)
axis(2, at = c(-8, -6, -4, -2, 0, 2, 4), pos = 0, cex.lab = 1.3, las = 1)
arrows(0,0, 1.8,0, angle = 10, length = 0.25)
arrows(0,-11, 0,7, angle = 10, length = 0.25)
text(1.8, 0.5, expression(x), cex = 1.4)
text(0.05, 7, expression(y), cex = 1.4)
text(0.4, -8, expression(y == x^3 + 4*x^2 -10), cex = 1.4)
#Plot the initial points and then follow tangent lines
x0 <- 1.0
x1 <- 1.4543
x2 <- 1.3689
lines(x,f(x0) + g(x0)*(x - x0),type = "l", lty = 2, col = "red")
text(0.73,-8.5, expression(paste("Tangent line through ", P[0])),
col = "red", cex = 1.2, srt = 42)
points(1,0, pch = 19)
text(1.0, 1.3, "Initial guess", cex = 1.4, col = "blue")
text(1.0, 0.5, expression(x[0]), cex = 1.4, col = "blue")
arrows(1,0,1,-5, angle = 10, length = 0.15, col = "blue")
points(x0, f(x0), col = "red", pch = 20)
text(x0+ 0.02, f(x0) - 0.5, expression(P[0]), cex = 1.4, col = "red")
arrows(x0,f(x0), x1,0, angle = 10, length = 0.15, col = "blue")
points(x1,0, pch = 19)
text(x1+0.04, 0.3, expression(x[1]), cex = 1.4, col = "blue")
arrows(x1,0,x1,f(x1), angle = 10, length = 0.15, col = "blue")
points(x1, f(x1), col = "red", pch = 20)
text(x1+ 0.05, f(x1), expression(P[1]), cex = 1.4, col = "red")
lines(x,f(x1) + g(x1)*(x - x1),type = "l", lty = 2, col = "red")
text(1.05,-6.5, expression(paste("Tangent line through ", P[1])),
col = "red", cex = 1.2, srt = 56)
arrows(x1,f(x1),x2,0, angle = 10, length = 0.15, col = "blue")
points(x2,0, pch = 19)
text(1.34,0.5, expression(x[2]), cex = 1.4, col = "blue")
text(0.8, 6, "Newton's method for finding a root:", cex = 1.8)
text(0.8, 5,"Follow the blue arrows from the initial guess",
cex = 1.3, col = "blue")
text(0.8, 4, expression(x[0] %=>% P[0] %=>% x[1] %=>% P[1] %=>% x[2] %=>% ...),
cex = 1.3, col = "blue") 
#dev.off()
Use Newton’s method to solve an equation
f <- function(x) { (x+1)^4 -(2+x) }
root <- newton(f, tol = 1E-12, x0 = 1, N = 10)
root## [1] 0.5809697 0.3335992 0.2359959 0.2210877 0.2207447 0.2207441 0.2207441
## [8] 0.2207441 0.2207441
Find multiple solutions of an EBM
S <- 1365
ep <- 0.6
sg <- 5.670373*10^(-8)
#Define the function for energy
f3 <- function(T){return(ep*sg*T^4 -
(1-(0.5 - 0.2 * tanh ((T-265)/10)))*(S/4))}
root1 <- newton(f3, tol = 1E-12, x0 = 220,N = 20)
root1## [1] 235.6965 234.3860 234.3817 234.3817 234.3817root2 <- newton(f3, tol = 1E-12, x0 = 270,N = 20)
root2## [1] 262.0567 264.5071 264.3378 264.3377 264.3377 264.3377root3 <- newton(f3, tol = 1E-12, x0 = 300,N = 20)
root3## [1] 289.9086 289.1469 289.1401 289.1401 289.1401 289.1401root5 <- newton(f3, tol = 1E-12, x0 = 100, N = 20)
root5## [1] 827.2544 623.5417 474.8968 372.5619 313.3648 292.0552 289.2231 289.1402
## [9] 289.1401 289.1401 289.1401
Use R to find derivatives of a function symbolically
D(expression(x^2),"x")## 2 * x# 2 * x
D(expression(exp(-x^2)),"x")## -(exp(-x^2) * (2 * x))# -(exp(-x^2) * (2 * x))
D(expression(sin(-3*t)-2*cos(4*t-0.3*pi)),"t")## -(cos(-3 * t) * 3 - 2 * (sin(4 * t - 0.3 * pi) * 4))# -(cos(-3 * t) * 3 - 2 * (sin(4 * t - 0.3 * pi) * 4))
D(expression(-g*t^2/2 + v0*t + h0),"t")## v0 - g * (2 * t)/2# v0 - g * (2 * t)/2
# Find derivative of this result function to find the second derivative
D(expression(v0 - g * (2 * t)/2),"t")## -(g * 2/2)# -(g * 2/2)
#or simply
D(D(expression(-g*t^2/2 + v0*t + h0),"t"),"t")## -(g * 2/2)# -(g * 2/2)
#The third-order derivative
D(D(D(expression(-g*t^2/2 + v0*t + h0),"t"),"t"),"t")## [1] 0