library(graphics) #load the graphics library needed for plotting
lx <- 6.0 #length of the computational domain
lt <- 5.0 #length of the simulation time interval
nx <- 300 #number of discrete lattice points
nt <- 200 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell
dt <- lt/nt #length of timestep
V = c(1:nx) #potential energies at discrete points
for(j in c(1:nx)) {
V[j] = as.complex(2*(j*dx-3)*(j*dx-3)) #Harmonic oscillator potential with k=4
}
kappa1 = (1i)*dt/(2*dx*dx) #an element needed for the matrices
kappa2 <- c(1:nx) #another element
for(j in c(1:nx)) {
kappa2[j] <- as.complex(kappa1*2*dx*dx*V[j])
}
psi = as.complex(c(1:nx)) #array for the wave function values
for(j in c(1:nx)) {
psi[j] = as.complex(exp(-(j*dx-2)*(j*dx-2))) #Gaussian initial wavefunction, displaced from equilibrium
}
xaxis <- c(1:nx)*dx #the x values corresponding to the discrete lattice points
for(m in c(1:nt)) {
A = matrix(nrow=nx,ncol=nx) #matrix for forward time evolution
B = matrix(nrow=nx,ncol=nx) #matrix for backward time evolution
for(j in c(1:nx)) {
kappa2[j] <- as.complex(kappa1*2*dx*dx*V[j])
for(k in c(1:j)) {
kappa2[j] = kappa2[j]+2*kappa1*2*dx*dx*dx*Im(psi[k]) # Add a nonlinear integral damping term to the potential
}
}
for(j in c(1:nx)) {
for(k in c(1:nx)) {
A[j,k]=0
B[j,k]=0
if(j==k) {
A[j,k] = 1 + 2*kappa1 + kappa2[j]
B[j,k] = 1 - 2*kappa1 - kappa2[j]
}
if((j==k+1) || (j==k-1)) {
A[j,k] = -kappa1
B[j,k] = kappa1
}
}
} #main time stepping loop
sol <- solve(A,B%*%psi) #solve the system of equations
for (l in c(1:nx)) {
psi[l] <- sol[l]
}
if(m %% 3 == 1) { #make plots of psi(x) on every third timestep
jpeg(file = paste("plot_",m,".jpg",sep=""))
plot(xaxis,abs(psi)^2, xlab="position (x)", ylab="Im(Psi)",ylim=c(-1.5,1.5))
title(paste("Abs(psi(x,t))^2 at t =",m*dt))
lines(xaxis,abs(psi)^2)
dev.off()
}
}