The numerical solution to the time-independent Schrödinger equation in Problem 2.61 can be extended to solve the time-dependent Schrödinger equation. When we discretize the variable x, we obtain the matrix equation
H \Psi=i \hbar \frac{d}{d t} \Psi . (11.149)
The solution to this equation can be written
\Psi (t+\Delta t)=\bigcup(\Delta t) \Psi(t) (11.150).
If H is time independent, the exact expression for the time-evolution operator is
\bigcup(\Delta t)=e^{-i H \Delta t / \hbar} (11.151).
and for Δt small enough, the time-evolution operator can be approximated as
\cup(\Delta t) \approx 1-i H \frac{\Delta t}{\hbar} (11.152).
While Equation 11.152 is the most obvious way to approximate U, a numerical scheme based on it is unstable, and it is preferable to use Cayley’s form for the approximation:
U (\Delta t) \approx \frac{1-\frac{1}{2} i \frac{\Delta t}{\hbar} H }{1+\frac{1}{2} i \frac{\Delta t}{\hbar} H } (11.153).
Combining Equations 11.153 and 11.150 we have
\left(1+\frac{1}{2} i \frac{\Delta t}{\hbar} H \right) \Psi (t+\Delta t)=\left(1-\frac{1}{2} i \frac{\Delta t}{\hbar} H \right) \Psi(t) (11.154).
This has the form of a matrix equation Mx = b which can be solved for the unknown x = \Psi (t+\Delta t) . Because the matrix M =1+\frac{1}{2} i \frac{\Delta t}{\hbar} H is tridiagonal, efficient algorithms exist for doing so.
(a) Show that the approximation in Equation 11.153 is accurate to second order. That is, show that Equations 11.151 and 11.153, expanded as power series in Δt , agree up through terms of order (\Delta t)^{2} . Verify that the matrix in Equation 11.153 is unitary.
As an example, consider a particle of mass m moving in one dimension in a simple harmonic oscillator potential. For the numerical part set m=1, \omega=1, \text { and } \hbar=1 (this just defines the units of mass, time, and length).
(b) Construct the Hamiltonian matrix H for N+1 =100 spatial grid points. Set the spatial boundaries where the dimensionless length is \xi=\pm 10 (far enough out that we can assume that the wave function vanishes there for low-energy states). By computer, find the lowest two eigenvalues of H, and compare the exact values. Plot the corresponding eigenfunctions. Are they normalized? If not, normalize them before doing part (c).
(c) Take \Psi(0)=\left(\psi_{0}+\psi_{1}\right) / \sqrt{2} (from part (b)) and use Equation 11.154 to evolve the wave function from time t=0 \text { to } t=4 \pi / \omega . Create a movie (Animate, in Mathematica) showing \operatorname{Re}(\Psi(t)), \operatorname{Im}(\Psi(t)) , and |\Psi(t)| , together with the exact result. Hint: You need to decide what to use for Δt. In terms of the number of time steps N_{t}, N_{t} \Delta t=4 \pi / \omega . In order for the approximation of the exponential to hold, we need to have E \Delta t / \hbar \ll 1 . The energy of our state is of order \hbar \omega , and therefore N_{t} \gg 4 \pi . So you will need at least (say) 100 time steps.
