Physics 6730 Assignment 10

In this assignment you will generate a solution to the time-dependent Schroedinger equation, using a second-order Crank-Nicholson algorithm adapted to complex partial differential equations. You will make a movie of your solution. And you will use Fourier transforms to analyze the results.

For each of the following exercises create the specified files with your answer(s). Submit your homework using the course submit utility.

Exercise 1 Discrete Complex Fourier Transform

This exercise works with code in ~p6730/examples/a10/cfourier.cc and the corresponding executable cfourier. You do not have to write the code.

The discrete complex Fourier transform converts between N complex values in the spatial (or time) domain and N complex values in the reciprocal wavelength (or frequency) domain. (To get wavenumber from reciprocal wavelength, multiply by 2 Pi. To get angular frequency from frequency, also multiply by 2 Pi.) The discrete reciprocal wavelength (or frequency) interval is 1/(N dx) where dx is the spatial interval. So the bookkeeping is slightly simpler than that of the real Fourier transform. In terms of reciprocal wavelength, the Fourier transform is periodic with period 1/dx. The standard interval is [-1/(2 dx),1/(2 dx)], but periodicity extends this range to the infinite real line.

We shall use the discrete complex Fourier transform to analyze the momentum content of wave functions defined on a regular discrete set of points in x.

To carry out the transform, we use the Numerical Recipes FFT subroutine four1. The first argument for this routine is in the "NRVec_DP" class.

The program run as follows 

   cfourier +1|-1 < wavef > ft
where wavef is a file containing the wavefunction (not necessarily with 2k points. The file has one complex value on each line as follows: 
  Re1 Im1
  Re2 Im2
   etc
The output file ft is in the format 
  p1 Re1 Im1
  p2 Re2 Im2
   etc
where pn is the reciprocal wavelength and Ren and Imn are the real and imaginary parts of the nth component of the Fourier transform. The single command line argument specifies either a forward transform (+1) or backward (-1). The bar | above indicates a choice.

Note that cfourier gives reciprocal wavelengths in units of 1/dx on the interval [0,1], so to get the standard interval [-1/2,1/2], use periodicity.

The Gaussian wave function ~p6730/examples/a10/gauss.0 has the form of a real Gaussian in coordinate space, multiplied by a plane wave factor exp(i k x). This is the conventional choice for a wavepacket moving in the positive x direction. Carry out the Fourier transform and determine:

  1. The value of k used to construct the wavefunction. Get this from the location of the peak in the plot of the squared norm of the Fourier components. In this way you should be able to determine the reciprocal wavelength k in units of 1/dx where dx is whatever we want for the spatial interval. Be careful to use the correct sign +1 or -1 with cfourier, remembering what was said above about the convention for motion in the positive x direction. We keep this sign for all momentum analyses.
  2. The center and width of the Gaussian wavefunction amplitude in coordinate space. Get this by plotting the squared norm of the wavefunction in the spatial domain, but remember that the standard deviation of the square is different from the standard deviation of the amplitude.
Report your answers in a file cfourier.txt to be handed in.

Exercise 2 Solving the Time-Dependent Schroedinger equation

Pseudocode for the second-order Crank-Nicholson algorithm for the time-independent Schroedinger equation equation is given in the online handout Crank-Nicholson Algorithm .

Using this algorithm, write a program barrier.cc and the associated Makefile to solve the time-dependent Schroedinger equation. That is, given an initial complex wave function, defined on the range x in [0,L], find the complex wave function at a specified later time. We take zero boundary conditions at the end points x = 0 and x = L.

For the potential use a square barrier with height V0 extending from x = x0 to x = x1 (integers). The height V0 is given in the units of the handout.

All computations are to be done in double precision. You will probabily find it convenient to use the same double_complex class as in Assignment 5.

  1. The program barrier must read the following values, one line at a time, from a file that you create, called barin: 
        nt		Number of time steps to carry out the integration (integer)
    m		Mass m
    V0 x0 x1 	Potential barrier height V0 and barrier end points
		(one double precision and two integer values)
		(The value of the potential is V0 at x0, x1, and all
		points between.  It is otherwise zero.)
    You needn't hand in this file. When we test your code, we will use our own.
  2. The initial wave function is taken from standard input. There are N lines (maximum 500) with the same format as in Exercise 1: 
            psi_real  psi_imag
    The 1st line has the value of the complex wave function at x = 1, the second, at x = 2, etc. The length of the interval in x is determined by the number of input values. The boundaries of the interval are at x = 0 and x = L, where it is assumed that the wave function vanishes. Thus L = N + 1, where N is the number of lines of input. Notice that you have N+2 complex values for psi, altogether.
  3. The final wave function is written to standard output and must have exactly the same format, so that it can be used subsequently as input to another computation.

Exercise 3 Digital Movies of Wavefunctions

The final output of this exercise is an MPEG format movie to hand in. This is the movie Lucas Films decided not to release!

Although the parameters for the movie that you hand in are prescribed, you are encouraged to experiment as well, since there is a great deal of interesting physics to be discovered with your program.

You will create a movie of 1500 frames that follows the evolution of the wave function over the time interval [0,3000]. Since this takes time, you should do a quick test with only 10 frames. Run your program with the initial wave function taken from the file in ~p6730/examples/a10/gauss.0. Use the parameter values 

 
     300 		time units
     10  		mass
     0.0075 60 65	V0 x0 x1 (barrier height and sides)
from the file ~p6730/examples/a10/barin. Put the output from your program in a file gauss.1 and plot the squared norm with the same procedure as in Exercise 2 above. You should see that the wave packet has shifted toward positive x, it has spread slightly, and the leading edge of the wave packet shows a little distortion from the encounter with the barrier. If you don't see this, check your program. Another check of your program is to set the barrier height to zero and see if the wave packet continues to move smoothly to the right until it encounters the boundary on the right edge at x = 101.

You should repeat this process for another 9 frames, creating the time sequence gauss.2 through gauss.10. A check of your code is that the total probability (found by summing the squared norm over all x) never changes in time (up to roundoff errors).

Keep the files gauss.1 through gauss.10 for use in Exercise 4.

After you are satisfied your code is working properly, change the file "barin" so it integrates only two time steps. Then from the directory ~p6730/examples/a10 copy the files make_movie.sh, make_movie.awk and mpe.param to the same directory as your files barin and gauss.0.

Run the script make_movie.sh. The result should be a file gauss.mpg that can be viewed with an MPEG player. On our Unix boxes you can use the command 

   mplayer gauss.mpg
The movie should take about one minute. It can also be viewed with Windows MPEG players.

Hand in the file gauss.mpg.

CAUTION: This script creates two working directories wavef and movie. The script removes these directories when it is finished. So please do not put any files you want to save in these directories.

Optional: You might try running this exercise with the potential turned off, with a potential well (negative barrier) and a much higher barrier. You could also try a smaller mass (say 1) and a shorter interval between frames to see the higher velocity and the difference in dispersion. Have fun! You might also change the initial wave function.

Exercise 4 Scattering in momentum space

It isn't always easy to understand all the physics by just staring at plots of squares of wavefunctions in coordinate space. So we use the discrete complex Fourier transform of Exercise 1 to analyze the scattering process in momentum space. In scattering from a barrier we expect that the reflected wave has momentum opposite the incident wave. So the objective here is to determine at what stage we start seeing reflected momenta. Also when the wave passes over the barrier, we might see some component with a reduced momentum, corresponding to the lower kinetic energy over the barrier.

Use the code cfourier from Exercise 1 to study the process in momentum space and describe your findings in a file barrier.txt to be handed in.

This page is maintained by:
Carleton DeTar Mail Form
Last modified 6 May 2004