Final edits

Schrick-Noah_Homework-3.R

@@ -40,9 +40,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, h){
-  n <- (tfinal-t0)/h # Static step size
-  
-  tspan <- seq(t0,tfinal,n)
+  tspan <- seq(t0,tfinal,h)
   npts <- length(tspan)
   y<-rep(y0,npts) # initialize, this will be a matrix for systems
   
@@ -141,7 +139,7 @@ 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="Numerical Decay Model Solution Using Pracma ODE45")
+     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")
@@ -218,7 +216,7 @@ plot.pred.prey <- function(sol, method){
 plot.pred.prey(pp.sol, "ODE45")
 
 # b) Use Euler and RK to solve with h=10
-h <- 10
+h <- 1
 
 my_euler2 <- function(f, y0, t0, tfinal, dt){
   # follow ode45 syntax, dt is extra
@@ -230,14 +228,12 @@ my_euler2 <- function(f, y0, t0, tfinal, dt){
   # 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
-    }
+  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))
+  return(list(t=tspan, y=y))
 }
 
 my_rk4_2 <- function(f, y0, t0, tfinal, h){
@@ -247,17 +243,15 @@ my_rk4_2 <- function(f, y0, t0, tfinal, h){
   # 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
+  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)
     
-      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)
-    }
+    y[i,] <- y_vec_prev + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
   }
   return(list(t=tspan,y=y))
 }
@@ -274,6 +268,10 @@ plot.pred.prey(pp.rk4.sol, "RK4")
 
 
 # plot comparing Prey solutions to ode45
+plot(pp.sol$t, pp.sol$y[,1], type="l", col="blue")
+par(new=T)
+lines(pp.euler.sol$t,pp.euler.sol$y[,1], col="red")
+
 
 # c) Use k3=0.02
 pp.params.c <- c(.01, .1, .02, .05)
@@ -307,7 +305,8 @@ 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]), 
+                                   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  
@@ -341,7 +340,7 @@ decomp.f <- function(t, y, k) {
 # 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
+B0 <- 0; C0 <- 0; D0 <- 0; E0 <- 0;
 tmin <- 0
 tmax <- 10
 
@@ -373,10 +372,17 @@ plot.decomp <- function(sol, method){
 plot.decomp(decomp.ode.sol, "ODE45")
 
 # b) Increase k4 to make the intermediate species -> 0
-k4 <- 10
-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)},
+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 = 100)
-plot.decomp(decomp.ode.sol.b, paste("ODE45 Using k4 =", k4))
-decomp.ode.sol.b$y[,5]
+                        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