# Solving Schrodinger's equation: Free particles and uncertainty # We'll be exploring the evolution of a Gaussian packet, corresponding # to the ground state of Hydrogen atom in a harmonic oscillator with # frequency one Terahertz. We start by defining some constants. hbar = 1.054571628251774e-27; omega = 1.e12; protonMass = 1.672621637e-24; m = protonMass; a0 = sqrt(...); # RMS width of Gaussian # We define psi0 on a lattice of Np points x, from -L/2 to L/2 L = ...; Np = ...; dx = ...; x = linspace(-L/2,L/2-dx,Np); # Now we define the initial wavefunction psi0 on our grid. # Remember to use x.^2 to do the square pointwise (rather than as # a matrix multiplication) psi0 = (...)^(1./4.) * exp(-...); plot(x, abs(psi0).^2) xlim([-2*a0,2*a0]) # If psi(k,t) is expressed in Fourier space, applying the kinetic energy is # pointwise multiplication. p=−i hbar k, so p^2/2m=−ℏ^2 k^2/2m. We convert # from real space into Fourier space using a FFT (Fast Fourier Transform), # which returns psi(k) at Np points separated by dk=2 pi/L: # k=[0,dk,2dk,…,(Np/2)dk,−(Np/2−1)dk,…,−2dk,−dk] # Notice that k grows until halfway down the list, and then jumps to # negative values and shrinks. We define k and k2=k2 as arrays # Note that 'zeros' creates a matrix, so for a vector we initialize # k=zeros(Np,1) # Note that Matlab/Octave arrays start at one, so k(1)=0, k(2)=dk, ... dk = ...; k = zeros(1,Np); for j = 0:(...) k(j+1) = ... * dk; end for j = (...):(Np-1) k(j+1) = (j-Np) * ...; end k2 = k.^2; # and we define the psi(t) using Fourier transforms. # I still haven't figured out how to define functions # like psi(t) without passing lots of variables like psi0, hbar, k2, ... P = ...; psiPby4 = ifft(exp(-...*(P/4.)/(...)).*fft(psi0)); plot(x,real(psiPby2), x, imag(psiPby2)) figure() psiP = ... plot(...) figure() ... # Calculate the width of our wavepackets at various times. # The width is given by the sum of x.^2 .* abs(psi).^2 * dx # Our initial wavepacket has RMS with DeltaX = a0. sqrt(sum(... .^2 .* abs(...).^2 * ...)) a0 # We calculate the widths for a series of times times = linspace(0,2*P,101); widths = 0.*times; for j = 1:length(times) psi = ... widths(j) = ... abs(psi).^2 ... end # We know that a minimum uncertainty wavepacket like ours has # Delta x Delta p=hbar/2, so we expect the packet width to grow like # vt with v given by the momentum uncertainty. We plot the comparison... deltaP = ... v = ... plot(times, widths, times, ...)