diff --git a/Schrick-Noah_Homework-3.R b/Schrick-Noah_Homework-3.R index 65dc567..4e9baf6 100644 --- a/Schrick-Noah_Homework-3.R +++ b/Schrick-Noah_Homework-3.R @@ -38,7 +38,7 @@ plot(abs(log(central.diff.table[,"h"],10)), central.diff.table[,"error"], main = "Central Difference Approximation of sin(x) at x=pi") ## 2. -# a) Runge-Kutta +# a) 4th-order Runge-Kutta my_rk4 <- function(f, y0, t0, tfinal){ n <- (tfinal-t0)/h # Static step size @@ -90,20 +90,44 @@ y0 <- 100 decay.euler.sol <- my_euler(f=function(t,y){decay.f(t,y,k=k)}, y0, t0, tfinal, dt=h) +# Exact sol +exact.y <- y0*exp(-k*seq(t0,tfinal,len=100)) +exact.y.slice <- exact.y[seq(1,length(exact.y),h)] +# Add in extra point to make lengths even to euler sol +exact.y.slice <- append(exact.y.slice,exact.y[length(exact.y)]) + 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") +lines(seq(t0,tfinal,len=100), + exact.y, col="red") +# Error +euler.err <- abs(decay.euler.sol$y - exact.y.slice) +euler.log.err <- log(euler.err, 10) +# Unsophisticated way to account for the -Inf error in the logscale +euler.log.err[which(!is.finite(euler.log.err))] <- NA +plot(seq(t0+h,tfinal,len=11), euler.log.err, + xlab="t", ylab="error (logscale)", type="o", + main = "Error from Euler Solution vs Exact") + +# Repeat with RK4 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") +lines(seq(t0,tfinal,len=100), + exact.y, col="red") +# Error +rk4.err <- abs(decay.rk4.sol$y - exact.y.slice) +rk4.log.err <- log(rk4.err, 10) +# Unsophisticated way to account for the -Inf error in the logscale +rk4.log.err[which(!is.finite(rk4.log.err))] <- NA + +plot(seq(t0+h,tfinal,len=11), rk4.log.err, + xlab="t", ylab="error (logscale)", type="o", + main = "Error from 4th-Order Runge-Kutta Solution vs Exact") ## 3. Use library function ode45 to solve the decay numerically if (!require("pracma")) install.packages("pracma")