Download Worksheet “Shooting Method” and more Study notes Particle Physics in PDF only on Docsity! PHYSICS 250 Spring 2005 P. NELSON The Shooting Method In this illustration I show how to solve the harmonic oscillator energy levels numerically. We begin by defining the potential energy function. Since computers don't use units I have selected units where the mass and spring constant, and hbar, all have numerical parts equal to one. (So the expected energy levels (n+1/2)*hbar*omega are just 1/2, 3/2, ...) VHO:=x->x^2/2; := VHO →x 1 2 x 2 Next I define the Schrodinger equation for energy level ET > Schr:=-diff(psi(x),x$2)/2+(VHO(x)-ET)*psi(x)=0; := Schr =− + 1 2 ∂ ∂ 2 x 2 ( )ψ x − 1 2 x 2 ET ( )ψ x 0 Now I need initial conditions to solve this equation. It's second order, so I need two conditions, one on psi and one on psi'. It turns out (ask me if you're curious) that half the solutions will be symmetric functions; that is, psi(-x)=psi(x). This implies that psi'(0)=0. I'm not interested in the normalization of psi, so any symmetric solution can be scaled until psi(0)=1. So the initial conditions for the symmetric solutions are > ICsym:=psi(0)=1,D(psi)(0)=0; := ICsym , =( )ψ 0 1 =( )( )D ψ 0 0 I need a certain package: > with(plots): Now I call a long gobbledygook (which you don't need to understand) to solve the equation for a chosen value of energy and plot the answer. I just happen to know the exact solution to this special problem, so I put in the right ET=0.5. Note the colon at the end suppresses the output. > P1:=odeplot(dsolve({subs(ET=0.5,Schr),ICsym},psi(x),numeric),[x,psi(x)], - 1 - 0..3,0..2,title=`E=0.5`,color=blue): The solution is in blue below. For reference I'll also plot the energy (black) and a line at ET (red). I scaled them both by a factor of 3 so they'd fit nicely. The thing to notice is that these lines cross at x=1, so any nonzero psi beyond that point is classically forbidden. > P2:=plot({VHO(x)/3,0.5/3},x=0..3,y=0..2): display({P1,P2}); x0 2.21.510.5 y 0 2 1.5 1 0.5 E=0.5 So there's lots of probability outside the classical turning point -- tunneling is possible. Next I show what happens if our guess for the energy value is a bit high: > odeplot(dsolve({subs(ET=0.55,Schr),ICsym},psi(x),numeric),[x,psi(x)],0..3,- 2..2,title=`0.55`); Error, (in dsolve) invalid arguments - 2 -
