Adding error plots for Runge-Kutta and Euler

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")