# Project 3 for the University of Tulsa's CS-7863 Sci-Stat Course
# Numerical Ordinary Differential Equations
# Professor: Dr. McKinney, Spring 2023
# Noah L. Schrick - 1492657

## 1. Approximate deriv. of sin(x) at x=pi from h=10^-1 -> 10^-20
forward.approx <- function(func, x, h){
  approx <- (func(x+h)-func(x))/h
}

forward.approx.table <- matrix(nrow = 0, ncol = 3)  
colnames(forward.approx.table) <- c("h", "approximation", "error")
for(h in 10^(seq(-1,-20,-1))){
  approx <- forward.approx(sin, pi, h)
  error <- abs(cos(pi)-approx)
  forward.approx.table <- rbind(forward.approx.table, c(h, approx, error)) 
}

plot(abs(log(forward.approx.table[,"h"],10)), forward.approx.table[,"error"],
     xlab="h [1e-x]", ylab="error (logscale)", type="o", log="y",
     main = "Forward Approximation of sin(x) at x=pi")


# Repeat with central difference approx
central.diff.approx <- function(func, x, h){
  approx <- (func(x+h)-func(x-h))/(2*h)
}
central.diff.table <- matrix(nrow = 0, ncol = 3)  
colnames(central.diff.table) <- c("h", "approximation", "error")
for(h in 10^(seq(-1,-20,-1))){
  approx <- central.diff.approx(sin, pi, h)
  error <- abs(cos(pi)-approx)
  central.diff.table <- rbind(central.diff.table, c(h, approx, error)) 
}

plot(abs(log(central.diff.table[,"h"],10)), central.diff.table[,"error"],
     xlab="h [1e-x]", ylab="error (logscale)", type="o", log="y",
     main = "Central Difference Approximation of sin(x) at x=pi")

## 2.
# a) 4th-order Runge-Kutta
my_rk4 <- function(f, y0, t0, tfinal, h){
  tspan <- seq(t0,tfinal,h)
  npts <- length(tspan)
  y<-rep(y0,npts) # initialize, this will be a matrix for systems
  
  for (i in seq(2,npts)){
    k1 <- f(t[i-1], y[i-1])
    k2 <- f(t[i-1] + 0.5 * h, y[i-1] + 0.5 * k1 * h)
    k3 <- f(t[i-1] + 0.5 * h, y[i-1] + 0.5 * k2 * h)
    k4 <- f(t[i-1] + h, y[i-1] + k3*h)
    # Update y
    y[i] = y[i-1] + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
    # Update t
    #t0 = t0 + h
  }
  return(list(t=tspan,y=y))
}

# b) Numerically solve the decay ode and compare to Euler and RK error
decay.f <- function(t,y,k){
  # define the ode model
  # k given value when solver called
  dydt <- -k*y
  as.matrix(dydt) 
}

# From class docs:
my_euler <- function(f, y0, t0, tfinal, dt){
  # follow ode45 syntax, dt is extra
  # y0 is initial condition
  tspan <- seq(t0,tfinal,dt)
  npts <- length(tspan)
  y<-rep(y0,npts) # initialize, this will be a matrix for systems
  for (i in 2:npts){
    dy <- f(t[i-1],y[i-1])*dt # t is not used in this case
    y[i] <- y[i-1] + dy
  }
  return(list(t=tspan,y=y))
}

# Solve
t0 <- 0
tfinal <- 100
h <- 10
k <- 0.03
y0 <- 100

decay.euler.sol <- my_euler(f=function(t,y){decay.f(t,y,k=k)},
                            y0, t0, tfinal, dt=h)
# Exact sol
exact.y <- y0*exp(-k*seq(t0,tfinal,len=100))
exact.y.slice <- exact.y[seq(1,length(exact.y),h)]
# Add in extra point to make lengths even to euler sol
exact.y.slice <- append(exact.y.slice,exact.y[length(exact.y)])

plot(decay.euler.sol$t,decay.euler.sol$y, xlab="t", ylab="y",
     main="Euler Numerical Solution for the decay ode dy/dt=-ky")
par(new=T)
lines(seq(t0,tfinal,len=100),
      exact.y, col="red")
# Error
euler.err <- abs(decay.euler.sol$y - exact.y.slice)
euler.log.err <- log(euler.err, 10)
# Unsophisticated way to account for the -Inf error in the logscale
euler.log.err[which(!is.finite(euler.log.err))] <- NA

plot(seq(t0+h,tfinal,len=11), euler.log.err,
     xlab="t", ylab="error (logscale)", type="o",
     main = "Error from Euler Solution vs Exact")

# Repeat with RK4
decay.rk4.sol <- my_rk4(f=function(t,y){decay.f(t,y,k=k)},
                          y0, t0, tfinal, h)
plot(decay.rk4.sol$t,decay.rk4.sol$y, xlab="t", ylab="y", 
     main="RK4 Numerical Solution for the decay ode dy/dt=-ky")
par(new=T)
lines(seq(t0,tfinal,len=100),
      exact.y, col="red")
# Error
rk4.err <- abs(decay.rk4.sol$y - exact.y.slice)
rk4.log.err <- log(rk4.err, 10)
# Unsophisticated way to account for the -Inf error in the logscale
rk4.log.err[which(!is.finite(rk4.log.err))] <- NA

plot(seq(t0+h,tfinal,len=11), rk4.log.err,
     xlab="t", ylab="error (logscale)", type="o",
     main = "Error from 4th-Order Runge-Kutta Solution vs Exact")

## 3. Use library function ode45 to solve the decay numerically
if (!require("pracma")) install.packages("pracma")
library(pracma)
# specify intial conds and params
k <- 1.21e-4
tmin <- 0
tmax <- 10000
y.init <- 10000

decay.ode45.sol <- ode45(f=function(t,y){decay.f(t,y,k=k)},
                         y=y.init, t0=tmin, tfinal=tmax)
plot(decay.ode45.sol$t,decay.ode45.sol$y, xlab="t", ylab="y",
     main="Decay Model Solution Using Pracma ODE45")
par(new=T)
lines(seq(tmin,tmax,len=50),
      y.init*exp(-k*seq(tmin,tmax,len=50)), col="red")
abline(h=y.init/2, lty="dashed", lwd=3)

#Solve for half-life with bisection
findZeroBisect <- function(func, xl, xr, tol, maxsteps=1e6){
  steps <- 0
  if (func(xl) * func(xr) >= 0){
    stop("Bad interval: try again with different endpoints")
  }
  
  repeat{
    xm <- (xl + xr)/2
    if(abs(func(xm)) < tol || steps >= maxsteps){
      break
    }
    
    ifelse(func(xl) * func(xm) > 0, xl<-xm, xr<-xm)
    steps <- steps + 1
  }
  return(c(xm, func(xm), steps))
}
fun.hl <- function(x) {-5000+(y.init)*exp(-k*x)}

halflife <- findZeroBisect(fun.hl, 5000, 6500, 1e-8)[1]
abline(v=halflife, lty="dotted", lwd=2, col="blue")

## 4. Solve the predator-prey model numerically
# pred-prey ode system with ode45
pp.f <- function(t,y,k){
  # k is a vector of model params
  # y is a vector of length equal to the number of odes
  prey <- y[1]
  pred <- y[2]
  dPrey <- -k[1]*prey*pred + k[2]*prey
  dPred <- k[3]*prey*pred - k[4]*pred
  return(as.matrix(c(dPrey, dPred)))
}

# a) k1=0.01, k2=0.1, k3=0.001, k4=0.05, prey(0)=50, pred(0)=15. t=0 -> 200
tmin <- 0; tmax<-200;
pp.params <- c(.01, .1, .001, .05)
prey0 <- 50
pred0 <- 15

pp.sol <- ode45(f = function(t,y){pp.f(t,y,k=pp.params)},
                  y = c(prey0, pred0),
                  t0 = tmin, tfinal = tmax)

# plot
if (!require("reshape")) install.packages("reshape")
library(reshape)
if (!require("ggplot2")) install.packages("ggplot2")
library(ggplot2)

plot.pred.prey <- function(sol, method){
  # - pred/prey columns melted to 1 column called "value" 
  sol.melted <- melt(data.frame(time=sol$t,
                                   Prey=sol$y[,1],
                                   Pred=sol$y[,2]), 
                        id = "time") # melt based on time
  # New column created with variable names, called “variable”
  colnames(sol.melted)[2] <- "Group" # used in legend  
  
  g <- ggplot(data = sol.melted, 
              aes(x = time, y = value, color = Group))
  g <- g + geom_point(size=2) 
  g <- g + xlab("time") +  ylab("Population") 
  g <- g + ggtitle(paste("Predator-Prey Model Using", method))  
  show(g)
}

plot.pred.prey(pp.sol, "ODE45")

# b) Use Euler and RK to solve with h=10
h <- 1

my_euler2 <- function(f, y0, t0, tfinal, dt){
  # follow ode45 syntax, dt is extra
  # y0 is a vector of initial conditions
  # this y determines the number of dependent variables to solve for
  tspan <- seq(t0,tfinal,dt)
  npts <- length(tspan)
  nvars <- length(y0)
  # intialize, this will be a matrix for systems
  y <- matrix(0,nrow=npts, ncol=nvars)
  y[1,] <- y0
  for (i in seq(2,npts)){
    y_prev <- y[i-1,]    # vector at previous time
    dy <- f(t[i-1],y_prev)*dt   # also a vector/matrix
    y[i,] <- y_prev + dy        # Euler update
  }
  return(list(t=tspan, y=y))
}

my_rk4_2 <- function(f, y0, t0, tfinal, h){
  tspan <- seq(t0,tfinal,h)
  npts <- length(tspan)
  nvars <- length(y0)
  # intialize, this will be a matrix for systems
  y <- matrix(0,nrow=npts, ncol=nvars)  
  y[1,] <- y0
  for (i in seq(2,npts)){
    y_vec_prev <- y[i-1,] # both y values at previous time point
  
    k1 <- f(t[i-1], y_vec_prev)
    k2 <- f(t[i-1] + 0.5 * h, y_vec_prev + 0.5 * k1 * h)
    k3 <- f(t[i-1] + 0.5 * h, y_vec_prev + 0.5 * k2 * h)
    k4 <- f(t[i-1] + h, y_vec_prev + k3*h)
    
    y[i,] <- y_vec_prev + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
  }
  return(list(t=tspan,y=y))
}

pp.euler.sol <- my_euler2(f = function(t,y){pp.f(t,y,k=pp.params)},
                         y = c(prey0, pred0),
                         t0 = tmin, tfinal = tmax, dt=h)
plot.pred.prey(pp.euler.sol, "Euler")

pp.rk4.sol <- my_rk4_2(f = function(t,y){pp.f(t,y,k=pp.params)},
                     y = c(prey0, pred0),
                     t0 = tmin, tfinal = tmax, h)
plot.pred.prey(pp.rk4.sol, "RK4")


# plot comparing Prey solutions to ode45
plot(pp.sol$t, pp.sol$y[,1], type="l", col="green", lwd=3, ylim=c(20,120),
     xlim=c(0,200), main="Comparing Prey solutions")
par(new=T)
lines(pp.euler.sol$t,pp.euler.sol$y[,1], col="red", lwd=2)
par(new=T)
lines(pp.rk4.sol$t,pp.rk4.sol$y[,1], col="blue", lty="dashed", lwd=2)
legend(x = "topleft",          
       legend = c("ODE45", "Euler", "RK4"),  
       lty = c(1, 1, 2),           
       col = c("green", "red", "blue"),           
       lwd = 2) 

# c) Use k3=0.02
pp.params.c <- c(.01, .1, .02, .05)
pp.sol.c <- ode45(f = function(t,y){pp.f(t,y,k=pp.params.c)},
                y = c(prey0, pred0),
                t0 = tmin, tfinal = tmax)
plot.pred.prey(pp.sol.c, "ODE45 with k3=0.02")

pp.params.c2 <- c(.01, .1, .02, .75)
pp.sol.c2 <- ode45(f = function(t,y){pp.f(t,y,k=pp.params.c2)},
                  y = c(prey0, pred0),
                  t0 = tmin, tfinal = tmax)
plot.pred.prey(pp.sol.c2, "ODE45 with k3=0.02, k4=0.75")

## 5. Solve SIR model numerically from t=0 -> 20
sir.f <- function(t, y, k) {
  dS <- -k[1] * y[1] * y[2]
  dI <-  k[1] * y[1] * y[2] - k[2] * y[2]
  dR <-  k[2] * y[2]
  return(as.matrix(c(dS, dI, dR)))
}

# a) a=0.5, b=1, S(0)=0.9, I(0)=0.1, R(0)=0
sir.params <- c(0.5, 1)
S0 <- 0.9
I0 <- 0.1
R0 <- 0.0
tmin <- 0
tmax <- 20

sir.ode.sol <- ode45(f = function(t,y){sir.f(t,y,k=sir.params)},
                y = c(S0, I0, R0),
                t0 = tmin, tfinal = tmax)

# plot
plot.sir <- function(sol, method){
  # - pred/prey columns melted to 1 column called "value" 
  sol.melted <- melt(data.frame(time=sol$t,
                                   S=sol$y[,1],
                                   I=sol$y[,2],
                                   R=sol$y[,3]), 
                        id = "time") # melt based on time
  # New column created with variable names, called “variable”
  colnames(sol.melted)[2] <- "Group" # used in legend  
  
  g <- ggplot(data = sol.melted, 
              aes(x = time, y = value, color = Group))
  g <- g + geom_point(size=2) 
  g <- g + xlab("time") +  ylab("Population") 
  g <- g + ggtitle(paste("SIR Model Using", method))  
  show(g)
}

plot.sir(sir.ode.sol, "ODE45")
# b) a=3
sir.params.b <- c(3, 1)
sir.ode.sol.b <- ode45(f = function(t,y){sir.f(t,y,k=sir.params.b)},
                     y = c(S0, I0, R0),
                     t0 = tmin, tfinal = tmax)
plot.sir(sir.ode.sol.b, "ODE45 with a=3")

## 6. Decomp of dinitrogen pentoxygen into nitrogen dioxide and molecular oxygen
decomp.f <- function(t, y, k) {
  dA <- -k[1]*y[1] + k[2]* y[2]*y[4]
  dB <- k[1]*y[1] - k[2]*y[2]*y[4] + 2*k[4]*y[4]*y[5]
  dC <- k[3]*y[2]*y[4]
  dD <- k[1]*y[1] - k[2]*y[2]*y[4] - k[3]*y[2]*y[4] - k[4]*y[4]*y[5]
  dE <- k[3]*y[2]*y[4] - k[4]*y[4]*y[5]
  return(as.matrix(c(dA, dB, dC, dD, dE)))
}

# a) k1=1.0, k2=0.5, k3=0.2,k4=1.5, [N2O5]o=1, all other IC's=0, t=0 -> 10
decomp.params <- c(1.0, 0.5, 0.2, 1.5)
A0 <- 1
B0 <- 0; C0 <- 0; D0 <- 0; E0 <- 0;
tmin <- 0
tmax <- 10

decomp.ode.sol <- ode45(f = function(t,y){decomp.f(t,y,k=decomp.params)},
                     y = c(A0, B0, C0, D0, E0),
                     t0 = tmin, tfinal = tmax)

# plot
plot.decomp <- function(sol, method){
  # - pred/prey columns melted to 1 column called "value" 
  sol.melted <- melt(data.frame(time=sol$t,
                                N2O5=sol$y[,1],
                                NO2=sol$y[,2],
                                O2=sol$y[,3],
                                NO3=sol$y[,4],
                                NO=sol$y[,5]),
                     id = "time") # melt based on time
  # New column created with variable names, called “variable”
  colnames(sol.melted)[2] <- "Group" # used in legend  
  
  g <- ggplot(data = sol.melted, 
              aes(x = time, y = value, color = Group))
  g <- g + geom_point(size=2) 
  g <- g + xlab("time [seconds]") +  ylab("Concentration") 
  g <- g + ggtitle(paste("Decomposition Model Using", method))  
  show(g)
}

plot.decomp(decomp.ode.sol, "ODE45")

# b) Increase k4 to make the intermediate species -> 0
k4.mono.table <- matrix(nrow = 0, ncol = 2)  
for(k in 10^(seq(-1, 4, 1))){
  k4 <- k
  decomp.params.b <- c(1.0, 0.5, 0.2, k4)
  decomp.ode.sol.b <- ode45(f = function(t,y){decomp.f(t,y,k=decomp.params.b)},
                        y = c(A0, B0, C0, D0, E0),
                        t0 = tmin, tfinal = tmax)
  #plot.decomp(decomp.ode.sol.b, paste("ODE45 Using k4 =", k4))
  # Is NO monotonically increasing?
  mono.check <- all(decomp.ode.sol.b$y[,5] == cummax(decomp.ode.sol.b$y[,5]))
  k4.mono.table <- rbind(k4.mono.table, c(k4, mono.check)) 
}

k4.mono.table