Pendulum problem using two first order ODEs

Schrick-Noah_Homework-4.R

@@ -6,18 +6,43 @@
 ## 1. Transform the second order ODE into a system of two first order ODEs
 
 # use ode45 or rk4sys to solve for the motion of the plane pendulum 
+if (!require("pracma")) install.packages("pracma")
+library(pracma)
 tmin <- 0
 tmax <- 7
-init_angle <- pi/4
-init_ang_spd <- 0
+theta0 <- pi/4
+omega0 <- 0
+y.init <- c(theta0, omega0)
 g <- 9.81
-L <- 1
+L <- 1.0
+
+## Substitution Tricks
+# theta_1 = theta_prime
+# theta_2 = theta
+# theta_1_prime = theta_pp
+# theta_2_prime = theta_prime = theta_1
+
+## Equation Replacements
+ode_1.f <- function(t, y, k){
+  theta <- y[1]
+  omega <- y[2]
+  g <- k[1]
+  L <- k[2]
+  omega_p <- -k[1]/k[2] * sin(theta)
+  theta_p <- omega
+  as.matrix(c(theta_p, omega_p)) 
+}
+
+p.params <- c(g, L)
+
+pend.sol <- ode45(f=function(t,y){ode_1.f(t,y,k=p.params)},
+                         y=y.init, t0=tmin, tfinal=tmax)
 
 # a. Plot
-plot(pend.sol[,1],pend.sol[,2],type="l",col="blue", 
+plot(pend.sol$t,pend.sol$y[,1],type="l",col="blue", 
      ylim=c(-2.3,2.3), xlab="time",ylab="amplitude")
 par(new=T)
-lines(pend.sol[,1],pend.sol[,3],type="l",col="red")
+lines(pend.sol$t,pend.sol$y[,2],type="l",col="red")
 abline(h=0)
 legend("topleft", # coordinates, "topleft" etc
        c("angle","omega"), # label
@@ -29,6 +54,20 @@ legend("topleft", # coordinates, "topleft" etc
 )
 
 # b. Overlay a plot of the small-angle approx
+plot(pend.sol$t,pend.sol$y[,1],type="l",col="blue", 
+     ylim=c(-2.3,2.3), xlab="time",ylab="amplitude")
+par(new=T)
+lines(seq(tmin,tmax,len=50),
+      theta0*cos(sqrt(g/L)*seq(tmin,tmax,len=50)),type="l",col="red")
+abline(h=0)
+legend("topleft", # coordinates, "topleft" etc
+       c("angle","small angle approx"), # label
+       lty=c(1,1), # line 
+       lwd=c(2.5,2.5), # weight
+       #cex=.8, 
+       bty="n", # no box
+       col=c("blue","red") # color
+)
 
 ## 2. Use ode45 to solve the damped harmonic oscillator 
 tmin <- 0