4th order Runge-Kutta and decay function solutions with plots

Schrick-Noah_Homework-3.R

@@ -39,8 +39,24 @@ plot(abs(log(central.diff.table[,"h"],10)), central.diff.table[,"error"],
 
 ## 2.
 # a) Runge-Kutta
-my_rk4 <- function(f, y0, t0, tfinal, dt){
+my_rk4 <- function(f, y0, t0, tfinal){
+  n <- (tfinal-t0)/h # Static step size
   
+  tspan <- seq(t0,tfinal,n)
+  npts <- length(tspan)
+  y<-rep(y0,npts) # initialize, this will be a matrix for systems
+  
+  for (i in seq(2,npts)){
+    k1 <- f(t[i-1], y[i-1])
+    k2 <- f(t[i-1] + 0.5 * h, y[i-1] + 0.5 * k1 * h)
+    k3 <- f(t[i-1] + 0.5 * h, y[i-1] + 0.5 * k2 * h)
+    k4 <- f(t[i-1] + h, y[i-1] + k3*h)
+    # Update y
+    y[i] = y[i-1] + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
+    # Update t
+    #t0 = t0 + h
+  }
+  return(list(t=tspan,y=y))
 }
 
 # b) Numerically solve the decay ode and compare to Euler and RK error
@@ -65,6 +81,27 @@ my_euler <- function(f, y0, t0, tfinal, dt){
   return(list(t=tspan,y=y))
 }
 
+# Solve
+t0 <- 0
+tfinal <- 100
+h <- 10
+k <- 0.03
+y0 <- 100
+
+decay.euler.sol <- my_euler(f=function(t,y){decay.f(t,y,k=k)},
+                            y0, t0, tfinal, dt=1)
+plot(decay.euler.sol$t,decay.euler.sol$y)
+par(new=T)
+lines(seq(t0,tfinal,len=50),
+      y0*exp(-k*seq(t0,tfinal,len=50)), col="red")
+
+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)
+par(new=T)
+lines(seq(t0,tfinal,len=50),
+      y0*exp(-k*seq(t0,tfinal,len=50)), col="red")
+
 ## 3. Use library function ode45 to solve the decay numerically
 if (!require("pracma")) install.packages("pracma")
 library(pracma)