First version of predator-prey model

Schrick-Noah_Homework-3.R

@@ -39,7 +39,7 @@ plot(abs(log(central.diff.table[,"h"],10)), central.diff.table[,"error"],
 
 ## 2.
 # a) 4th-order Runge-Kutta
-my_rk4 <- function(f, y0, t0, tfinal){
+my_rk4 <- function(f, y0, t0, tfinal, h){
   n <- (tfinal-t0)/h # Static step size
   
   tspan <- seq(t0,tfinal,n)
@@ -113,7 +113,7 @@ plot(seq(t0+h,tfinal,len=11), euler.log.err,
 
 # Repeat with RK4
 decay.rk4.sol <- my_rk4(f=function(t,y){decay.f(t,y,k=k)},
-                          y0, t0, tfinal)
+                          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)
@@ -171,16 +171,115 @@ 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 <- 10
+
+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 2:npts){ # rows are time
+    for (j in 1:nvars){
+      y_vec_prev <- y[i-1,] # both y values at previous time point
+    }
+    dy <- f(t[i-1],y_vec_prev)*dt # f returns a vector, t not used
+    y[i,] <- y_vec_prev + dy
+  }
+  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 2:npts){ # rows are time
+    for (j in 1:nvars){
+      y_vec_prev <- y[i-1,] # both y values at previous time point
+    }
+    k1 <- f(t[i-1], y_vec_prev[i-1])
+    k2 <- f(t[i-1] + 0.5 * h, y_vec_prev[i-1] + 0.5 * k1 * h)
+    k3 <- f(t[i-1] + 0.5 * h, y_vec_prev[i-1] + 0.5 * k2 * h)
+    k4 <- f(t[i-1] + h, y_vec_prev[i-1] + k3*h)
+    y[i,] <- y_vec_prev[i-1] + (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
 
 # 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)},
+                y = c(prey0, pred0),
+                t0 = tmin, tfinal = tmax)
+plot.pred.prey(pp.sol.c, "ODE45 with k3=0.02")
 
 ## 5. Solve SIR model numerically from t=0 -> 20