Click here for other works of the author on RPubs

Load package

Load knitr for presenting tables in report and purrr for manipulating list

library(knitr)
library(purrr)

1. Write a function to do root-finding based on Newton-Ralphson method and solve for \(\frac{sin(x)}{x} - 0.6\). (Note: write your own codes, set tol=0.000001, try different initial values)

Define function \(y = \frac{sin(x)}{x} - 0.6\) and its derivative \(\frac{x*cos(x) - sin(x)}{x^2}\)

y = function(x) sin(x) / x - 0.6
y_dx = function(x) (x*cos(x) - sin(x))/(x^2)

Plot the objective function \(\frac{sin(x)}{x} - 0.6\) and see for which values of x does it equal to 0

#plot the function y
x = seq(-10, 10, length = 100)
plot(x, y(x), type = "l")

#add lines to indicate the location of roots
abline(h = 0, col = "red")
lines(x = rep(-1.660035, 100), y = seq(-0.9, 0, length = 100), lty = 2)
lines(x = rep( 1.660035, 100), y = seq(-0.9, 0, length = 100), lty = 2)
axis(1, c(-1.66, 1.66))

It can be observed that the function has two roots: -1.66 and 1.66

Use the Newton-Ralphson method to find the root(s)

Define function NR that conducts a Newton-Ralphson method search to find a root for the function y

NR <- function(ini = 0.1, tol = 0.000001){
    x_i = ini
    improve = 1
    
    #keep searching until improvement is lower than the specified tolerance
    while(improve > tol){
        x_ori = x_i
        x_i = x_ori - y(x_ori)/y_dx(x_ori)
        improve = abs((x_i - x_ori) / x_ori)
        
        #break out of the loop if NaN occurs
        if(is.nan(improve)){
            x_i = NA
            break
        }
    }
    return(x_i)
}

#try different initial values
ini_values = seq(-5, 5, by = 0.5)
NR_result = sapply(ini_values, NR)
kable(cbind(ini_values, NR_result), col.names = c("Initial value", "Root found"), caption = "Root found using Newton-Ralphson method with different initial values")
Root found using Newton-Ralphson method with different initial values
Initial value Root found
-5.0 1.660035
-4.5 1.660035
-4.0 1.660035
-3.5 -1.660035
-3.0 -1.660035
-2.5 -1.660035
-2.0 -1.660035
-1.5 -1.660035
-1.0 -1.660035
-0.5 -1.660035
0.0 NA
0.5 1.660035
1.0 1.660035
1.5 1.660035
2.0 1.660035
2.5 1.660035
3.0 1.660035
3.5 1.660035
4.0 -1.660035
4.5 -1.660035
5.0 -1.660035

Depending on the initial value given, the Newton-Ralphson method finds one of the two roots, -1.66 or 1.66. However, the method fails at 0, the inflection point(\(\frac{d}{dx}f(x)=0\)) of the function \(\frac{sin(x)}{x} - 0.6\)

2. 2. Use data from Vidal (1980) and find the Belehradek’s equation for C2, C3, C4, C5 by minimizing the least square error, and set b = -2.05. Plot the data and fitted curves.

Import data

data <- read.table("VidalTvsDuration.txt", header = T)
kable(data, col.names = c("Temperature", "C2", "C3", "C4", "C5"), caption = "Development time data")
Development time data
Temperature C2 C3 C4 C5
15.5 1.93 4.40 7.86 12.68
12.0 2.63 5.81 9.75 16.00
8.0 5.11 10.54 16.70 25.62

Define function D that uses the Belehradek’s equation to calculate the zooplankton development time and function LSE to calculate the sum of least square error of our estimation

b = -2.05
D = function(t, para) para[1]*(t - para[2])^b
LSE = function(obs, ...) sum((D(...) - obs)^2)

Find the optimal parameters that minimizes the least square error for C2, C3, C4 and C5 using optim

par=list()
for(i in 1:4){
    par[[i]] = optim(c(0.5, 0.5), LSE, t = data[, 1], obs = data[, 1 + i])
}

Report the esimated parameters for C2, C3, C4 and C5 in the Belehradek’s equation

#extract parameters
coef = t(matrix(unlist(transpose(par)$par), 2, 4))
rownames(coef) = c("C2", "C3", "C4", "C5")

#show results
kable(coef, caption = "Esimated parameters for C2, C3, C4 and C5", col.names = c("a", expression(alpha)))
Esimated parameters for C2, C3, C4 and C5
a alpha
C2 752.1069 -3.440511
C3 2051.1450 -5.114860
C4 4515.2148 -7.418156
C5 8756.0723 -9.268106

Plot the data and fitted curves

#plot original data
x = seq(min(data[, 1]), max(data[, 1]), length = 100)
plot(rep(data[, 1], 4), unlist(data[-1]), xlab = "Temperature", ylab = "Stage duration", pch = 16, col = "red")

#add fitted curves for C2, C3, C4 and C5
lines(x, D(x, par[[1]][[1]]), type = "l", col = "blue")
lines(x, D(x, par[[2]][[1]]), type = "l", col = "blue")
lines(x, D(x, par[[3]][[1]]), type = "l", col = "blue")
lines(x, D(x, par[[4]][[1]]), type = "l", col = "blue")