diff --git a/Schrick-Noah_Homework-3.R b/Schrick-Noah_Homework-3.R index 11f2965..4054d05 100644 --- a/Schrick-Noah_Homework-3.R +++ b/Schrick-Noah_Homework-3.R @@ -39,10 +39,47 @@ plot(abs(log(central.diff.table[,"h"],10)), central.diff.table[,"error"], ## 2. # a) Runge-Kutta +my_rk4 <- function(f, y0, t0, tfinal, dt){ + +} # 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)) +} ## 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