# 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) Runge-Kutta
my_rk4 <- function(f, y0, t0, tfinal){
  n <- (tfinal-t0)/h # Static step size
  
  tspan <- seq(t0,tfinal,n)
  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)
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=50),
      y0*exp(-k*seq(t0,tfinal,len=50)), col="red")

decay.rk4.sol <- my_rk4(f=function(t,y){decay.f(t,y,k=k)},
                          y0, t0, tfinal)
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=50),
      y0*exp(-k*seq(t0,tfinal,len=50)), col="red")

## 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)
par(new=T)
lines(seq(tmin,tmax,len=50),
      y.init*exp(-k*seq(tmin,tmax,len=50)))

## 4. Solve the predator-prey model numerically

# a) k1=0.01, k2=0.1, k3=0.001, k4=0.05, prey(0)=50, pred(0)=15. t=0 -> 200

# plot

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

# plot comparing Prey solutions to ode45

# c) Use k3=0.02

## 5. Solve SIR model numerically from t=0 -> 20

# a) a=0.5, b=1, S(0)=0.9, I(0)=0.1, R(0)=0

# plot

# b) a=3

## 6. Decomp of dinitrogen pentoxygen into nitrogen dioxide and molecular oxygen

# a) k1=1.0, k2=0.5, k3=0.2,k4=1.5, [N2O5]o=1, all other IC's=0, t=0 -> 10

# b) Increase k4 to make the intermediate species -> 0