{ "cells": [ { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.760567191632131 8.750311790315656\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Variational calculation with truncated Hilbert space\n", "# - numerical diagonalization of H in basis of N partocle-in-a-box states\n", "\n", " using QuadGK\n", " using LinearAlgebra\n", "\n", "# Any potential within a box from x=0 tp x=ll can be pt here\n", "\n", " function potential(x::Float64,ll::Float64)::Float64\n", " v::Float64=10\n", " return v\n", " end \n", "\n", "# Particle-in-a-box wave function k\n", "\n", " function eigenfunc(x::Float64,k::Int,ll::Float64)::Float64\n", " f=sin(k*x*pi/ll)*(2/ll)^0.5\n", " return f\n", " end\n", "\n", "# The function to be integrated for the potential-energy matrix elements\n", "\n", " function vpsikpsip(x::Float64,k::Int,p::Int,ll::Float64)::Float64\n", " psik=eigenfunc(x,k,ll)\n", " psip=eigenfunc(x,p,ll)\n", " return psik*psip*potential(x,ll)\n", " end\n", "\n", "# Integrates the function func(x) from x=x2 to x=x2\n", "# - uses n intervals dx (n+1 grid points)\n", "\n", " function simpson(func,x1::Float64,x2::Float64,n::Int)::Float64\n", " dx=(x2-x1)/n\n", " si::Float64=func(x1)+func(x2)\n", " for i=1:2:n-3\n", " x=x1+i*dx\n", " si=si+4*func(x)+2*func(x+dx)\n", " end \n", " si=si+4*func(x2-dx)\n", " return si*dx/3\n", " end\n", "\n", "# Makes the Hamiltonian matrix\n", "# The commented-out lines are for integration using quadgk()\n", "# - this turns out to be very slow, use Simpson's method instead\n", "\n", " function hamiltonian(n::Int,ll::Float64)\n", " ham=Matrix{Float64}(undef,n,n)\n", " eps=1e-4\n", " a=0.5\n", " x1::Float64=ll/2-a*(ll/2)\n", " x2::Float64=ll/2+a*(ll/2)\n", " for p=1:n\n", " for k=1:n\n", "# ham[k,p],e = quadgk(x -> vpsikpsip(x,k,p,ll),x1,x2,rtol=eps)\n", " ham[k,p] = simpson(x -> vpsikpsip(x,k,p,ll),x1,x2,1000)\n", " end\n", " ham[p,p] += (0.5*(pi*p/ll)^2)\n", " end\n", " return ham\n", " end\n", "\n", "# Transforms from the eigenvector in the box basis to real space\n", "\n", " function psi(v,x)\n", " nx=length(x)\n", " ll=x[nx]\n", " p=zeros(Float64,nx)\n", " for j=1:nx\n", " for i=1:length(v)\n", " p[j] += v[i]*eigenfunc(x[j],i,ll)\n", " end\n", " end\n", " return p\n", " end\n", "\n", " n=100 # number of basis states\n", " ll=2. # box size (between x=0 and x=ll)\n", "\n", " h=hamiltonian(n,ll)\n", "\n", "# Get the eigenvalues and eigenvectors\n", "# - there is another function eigen() that does both, may be faster\n", "\n", " e=eigvals(h)\n", " v=eigvecs(h)\n", "\n", " nx=200 # grid points (-1) for wave function\n", " dx=ll/nx # grid spacing for wave function\n", " x=zeros(Float64,nx+1) # vector holding the grid points for the wave function \n", " for j=1:nx+1\n", " x[j]=(j-1)*dx\n", " end\n", "\n", " k=1 # The eigenfunction to be graphed (k in 1,...,n) \n", " println(e[k],\" \",e[k+1])\n", " p=psi(v[:,k],x)\n", "\n", " using Plots\n", " \n", " plot(x,p)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.6.1", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.1" } }, "nbformat": 4, "nbformat_minor": 4 }