From 5297dafedb7ea9ea1d3ca85924782a9f29ff808c Mon Sep 17 00:00:00 2001
From: noah <noah@theschricks.com>
Date: Mon, 27 Mar 2023 15:46:10 -0500
Subject: [PATCH] Quantum harmonic oscillator

---
 Schrick-Noah_Homework-5.R | 60 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 60 insertions(+)

diff --git a/Schrick-Noah_Homework-5.R b/Schrick-Noah_Homework-5.R
index 1014a21..2437e18 100644
--- a/Schrick-Noah_Homework-5.R
+++ b/Schrick-Noah_Homework-5.R
@@ -81,9 +81,69 @@ i.b <- solve(A.3b,b.3vec)
 i.b
 
 ## 4. Schrodinger eq
+if (!require("pracma")) install.packages("pracma")
+library(pracma)
+
 #a. Solve and plot the 1d quantum harmonic oscillator wave function
+quantumHO <- function(x,y,E){
+  psi     <- y[1]
+  y1p <- y[2]   # y1’ fill in blank
+  y2p <- (x^2)*y[1]-2*E*y[1]  # y2’ fill in blank
+
+  as.matrix(c(y1p,y2p)) 
+}
+
+if (!require("reshape")) install.packages("reshape")
+library(reshape)
+if (!require("ggplot2")) install.packages("ggplot2")
+library(ggplot2)
+plotQHO <- function(E){
+  # plot the wavefunction for a given energy
+  QHO.sol <- ode45(function(x,y){quantumHO(x,y,E)}, 
+                   y0=c(1,0),t0=0, tfinal=5)
+  QHO.melted <- melt(data.frame(position=QHO.sol$t,
+                                psi=QHO.sol$y[,1],
+                                psi_p=QHO.sol$y[,2]), 
+                     id = "position")
+  colnames(QHO.melted)[2] <- "Solutions"  # rename "variable" 
+  
+  g <- ggplot(data = QHO.melted, aes(x=position, y=value, color=Solutions))
+  g <- g + geom_line(size=1) + geom_abline(slope=0, intercept=0)
+  g <- g + xlab("position") +  ylab("Amplitude") 
+  g <- g + ggtitle("Quantum Harmonic Oscillator")  
+  show(g)
+}
+
+plotQHO(0.5)
 
 # b. Solve with Shooting Method
+## Shooting Method Objective Function
+QHO_shoot <- function(E){
+  # E is an energy guess
+  ho.sol <- ode45(function(x,y){quantumHO(x,y,E)}, 
+                  y0=c(1,0),t0=0, tfinal=5)
+  last_idx <- length(ho.sol$t)
+  psi_right <- ho.sol$y[,1][last_idx] # psi at last x
+  # We want the wave function at the right boundary to be 
+  # exponentially small. So, the log-abs value will be large 
+  # and negative for the correct solution.
+  return(log(abs(psi_right)))  
+}
+
+## The objective function plot informs the energy guesses below.
+energy_vals <- seq(0,5.,len=100)
+objective_vals <- sapply(energy_vals, QHO_shoot)
+plot(energy_vals,objective_vals, xlab="Energy", ylab="Objective")
+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
+
+# Find first excited state E1.
+E1<-optimize(QHO_shoot,interval=c(2,3),tol=1e-8)$minimum
+E1
+plotQHO(E1)
 
 ## 5. Solve harmonic oscillator Schrodinger with finite difference
 # a. Solve with R and Julia