diff --git a/Schrick-Noah_Homework-3.R b/Schrick-Noah_Homework-3.R index 4e9baf6..c4b1d11 100644 --- a/Schrick-Noah_Homework-3.R +++ b/Schrick-Noah_Homework-3.R @@ -140,10 +140,35 @@ 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) +plot(decay.ode45.sol$t,decay.ode45.sol$y, xlab="t", ylab="y", + main="Numerical Decay Model Solution Using Pracma ODE45") par(new=T) lines(seq(tmin,tmax,len=50), - y.init*exp(-k*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