Final R code

Schrick-Noah_Homework-5.R

@@ -32,6 +32,7 @@ library(matlib)
 
 # verify
 bvec <- solve(A.2, yvec)
+bvec
 all.equal(as.matrix(yvec), A.2 %*% bvec) # use all.equal instead of == (tols and storage type)
 # EX: identical(as.double(8), as.integer(8)) returns FALSE
 
@@ -139,6 +140,7 @@ grid()
 # “optimize” finds min or max of function in an interval.
 # Find ground state energy E0.
 E0<-optimize(QHO_shoot,interval=c(0,1),tol=1e-8)$minimum
+E0
 
 # Find first excited state E1.
 E1<-optimize(QHO_shoot,interval=c(2,3),tol=1e-8)$minimum
@@ -147,10 +149,106 @@ plotQHO(E1)
 
 ## 5. Solve harmonic oscillator Schrodinger with finite difference
 # a. Solve with R and Julia
+if (!require("tictoc")) install.packages("tictoc")
+library(tictoc)
+
+solve_QHO <- function(n, xL, xR){
+  #n<-1000
+  #xL <- -10  # need to change for very excited state
+  #xR <- 10
+  h<- (xR-xL)/(n-2)  # step size
+  x <- seq(xL, xR, by=h) # grid of x values
+  
+  H <- diag(1/h^2 + .5*x^2) # makes the diagnonal
+  dim(H)
+  #   diag       upper diag              lower diag
+  H <- H + Diag(rep(-1/(2*h^2),n-2),1) + Diag(rep(-1/(2*h^2),n-2),-1)
+  
+  H.eigs <- eigen(H)
+  vals <- H.eigs$values
+  vecs <- H.eigs$vectors
+  
+  return(list(x=x,vecs=vecs,vals=vals))
+}
+
+print_QHO_sol <- function(x,vecs,vals,n){
+  print("Smallest eigvalue, lowest energy:")
+  print(vals[n-1])
+  #vals[n-1]  # smallest eigvalue, lowest energy
+  plot(x,vecs[,n-1], type="l", main="Ground state wave function") # ground state wave function
+  
+  print("2nd smallest eigvalue, first excited energy:")
+  print(vals[n-2])
+  #vals[n-2]  # 2nd smallest eigvalue, first excited energy
+  plot(x,vecs[,n-2], type="l", main="First excited wave function") # first excited wave function
+  
+  print("3rd smallest eigvalue, second excited energy:")
+  print(vals[n-3])
+  #vals[n-3]  # 3rd smallest eigvalue, second excited energy
+  plot(x,vecs[,n-3], type="l", main="Second excited wave function") # second excited wave function
+}
+
+# data for timing comparison
+timing.table <- matrix(nrow = 0, ncol = 4)  
+for (n in seq(500, 3500, by=500)){
+  for (x in seq(1,10)){
+    tic(quiet=TRUE)
+    QHO_sol <- solve_QHO(n, -x, x)
+    t <- toc(quiet=TRUE)
+    timing.table <- rbind(timing.table, c(-x, x, n, t$callback_msg)) 
+  }
+}
+
+# answer to main question
+QHO_sol <- solve_QHO(1000, -5, 5)
+print_QHO_sol(QHO_sol$x, QHO_sol$vecs,QHO_sol$vals,1000)
 
 # b. Solve damped oscillator with perturbation param
-# Plot
+yo <- 0
+y1 <- 1
+tmin <- 0
+tfinal <- 2
 
-## 6. Solve heat diffusion equation
+dho.pert.encap.f <- function(ep){
+  dho.pert.f <- function(t, y){
+    yn <- y[1]
+    ydot <- y[2] - yn
+    ypp <- (-(1+ep)*yn - ydot)/ep
+    as.matrix(c(ydot,ypp))
+  }
+  
+  proj.obj <- function(v0, y0=yo, tfinal){
+    # minimize w.r.t. v0
+    proj.sol <- ode45(dho.pert.f, 
+                      y=c(yo, y1), t0=tmin, tfinal=tfinal)
+    final_index <- length(proj.sol$t)
+    yf <- proj.sol$y[final_index,1]  # want equal to right boundary 
+    log(abs(yf)) # minimize this
+  }
+  
+  # user specifies tfinal and yfinal for BVP
+  ydot_best <- optimize(proj.obj, 
+                        interval=c(1,100), #bisect-esque interval
+                        tol=1e-10, 
+                        y0=0, tfinal=2) # un-optimized obj params
+  ydot_best$minimum  # best ydot
+  
+  best.sol <- rk4sys(dho.pert.f, a=0, b=2, y0=c(0, ydot_best$minimum),
+                     n=20)  # 20 integration stepstmax 
+}
 
-# EC. Solve the square plate problem
+approx <- function(x){
+  exp(1-(x/ep)) + exp(1-x)
+}
+
+ep <- 0.5
+ep.a <- dho.pert.encap.f(ep)
+plot(ep.a$x, ep.a$y[,1], type="l", main="Epsilon=0.5")
+par(new=T)
+lines(ep.a$x, approx(seq(0,2,len=21)), type="l", col="red")
+
+ep <- 0.05
+ep.b <- dho.pert.encap.f(ep)
+plot(ep.b$x, ep.b$y[,1], type="l", main="Epsilon=0.05")
+par(new=T)
+lines(ep.b$x, approx(seq(0,2,len=21)), type="l", col="red")