next up previous contents
Next: Numerical Evaluation of the Up: Finite Differences Previous: An Algorithm Based on   Contents

The GOLDBERG Algorithm

In the TROTTER-based algorithm, the position and momentum variables remain nondiscretized; only through the FFT -- which essentially is just a sophisticated formulation of the discrete FOURIER transformation using a discrete set of nodes in position and momentum space -- effectively $x$ and $p$ become discretized, too.

Now I describe a more typical finite differences approach to solving the SCHRÖDINGER equation for the free harmonic oscillator in the position representation, where both space and time are discretized in equidistant steps from the beginning. The solution $\left< x \left\vert \,\psi(t) \right> \right.$ on the interval $[x_{\mbox{\scriptsize min}},x_{\mbox{\scriptsize max}}]$ is constructed at the nodes

x_j & = & x_{\mbox{\scriptsize min}} + ...
...tsize max}}} & = & x_{\mbox{\scriptsize max}} \; ,
\end{array}\end{displaymath} (3.10)

and the time steps
t_n^{(m)} \; = \; nT + m \, \delta t, \qquad m=0,\dots,s
\end{displaymath} (3.11)

are used. For simplicity, only the dynamics in the time interval $\big[nT,(n+1)T\big)$ is considered here, such that the subscript $n$ can be dropped. The wave function at site $x_j$ and time $t^{(m)}$, and the potential contribution to $\H_{\mbox{\scriptsize free}}$ at $x_j$ are denoted as
$\displaystyle \psi_j^{(m)}$ $\textstyle :=$ $\displaystyle % \psi\!\left(x_j,t^{(m)}\right)
\left< x_j \left\vert \,\psi\left(t^{(m)}\right) \right> \right.$ (3.12)
$\displaystyle v_j \hspace*{0.45cm}$ $\textstyle :=$ $\displaystyle \frac{1}{2}x_j^2,$ (3.13)


In a first step, ${\hat{U}}_{\mbox{\scriptsize kick}}$ is applied to the discretized wave function:

\psi_j^{(0)} \; \longmapsto \; e^{ \textstyle -\frac{i}{\hbar}V_0\cos x_j } \,
\psi_j^{(0)} .
\end{displaymath} (3.14)

Then, the new $\psi_j^{(0)}$ is propagated for a period of time of length $T$ using the GOLDBERG algorithm [Koo86,KW96] that I describe now in some more detail.

With the discrete first order approximation for the second derivative [SB00]

\Big< x_j \Big\vert \,\frac{\partial^2}{\partial x^2}\, \Big...
+ {\mathcal O}\left( (\delta x)^2 \right) ,
\end{displaymath} (3.15)

the application of the free harmonic oscillator Hamiltonian to $\psi_j^{(m)}$ can be written as
\H_{\mbox{\scriptsize free}}\,\psi_j^{(...
...}\left( \left(\hbar \,\delta x \right)^2 \right) .
\end{array}\end{displaymath} (3.16)

Using this notation, the propagation over the time interval $\delta t$ is approximately
\psi_j^{(m+1)} \; = \; e^{ \textstyle -\frac{i}{\hbar}\H_{\mbox{\scriptsize free}} \,\delta t } \,
\end{displaymath} (3.17)

As in the previous subsection it remains to replace the exponential with a suitable approximant, and again using the TAYLOR formula (3.2) is not a good choice, as it does not yield the desired unitary expression. Rather, the GOLDBERG algorithm calls for the CAYLEY form [PTVF94] of the approximant:

e^{ \textstyle -\frac{i}{\hbar}\H_{\mbox{\scriptsize free}} ...
\!\left(\frac{\delta t}{\hbar}\right)^{\!\!\!3}
\bigg) \, .
\end{displaymath} (3.18)

Clearly, the quotient on the right hand side is unitary, thereby in principle allowing for an unconditionally stable algorithm. But note that the approximation (3.16) for the application of $\H_{\mbox{\scriptsize free}}$ can still spoil the stability of the algorithm. The CAYLEY approximant is correct up to second order in $\delta t/\hbar$, one power more than would have been supposed at first sight, as both the numerator and the denominator in equation (3.18) are first order approximations only. On the other hand, this very desirable feature of the approximant is accompanied by the less desirable feature that it results in a set of implicit difference equations to be solved -- see equation (3.20) below.

With the abbreviations
\tilde{v}_j & := & \frac{2 (\delta x)^2}{\...
...\frac{4 (\delta x)^2}{\hbar \, \delta t} \; ,
combination of equations (3.16-3.18) gives the iteration scheme

\psi_{j-1}^{(m+1)} +
\left( i\lambda-\tilde{v}_j \right) \p...
...da+\tilde{v}_j \right) \psi_{j }^{(m )} -
\psi_{j+1}^{(m )} .
\end{displaymath} (3.18)

It is to be solved for the $\psi_{j}^{(m+1)}$ and essentially requires the inversion of a tridiagonal system of linear equations. This is seen more easily in vectorial notation. Equation (3.20) is equivalent to
A_+\vec{\psi}^{\, (m+1)}
\; = \; A_-\vec{\psi}^{\, (m )},
\end{displaymath} (3.19)

for which the vectors
\vec{\psi}^{\, (m)}
\; := \; \left(
\psi_{j_{\mbox{\scriptsize max}}}^{(m)}
\end{array} \right)
\end{displaymath} (3.20)

and the symmetric $(j_{\mbox{\scriptsize max}},j_{\mbox{\scriptsize max}})$ matrices
\; := \; \left(
...\mp\tilde{v}_{j_{\mbox{\scriptsize max}}}
\end{array} \right)
\end{displaymath} (3.21)

have been defined. The value of the parameter $c$ depends on the (DIRICHLET) boundary conditions to be used.

% latex2html id marker 9008\begin{itemize}
...le periodicity interval
in position space.

The right hand side of equation (3.21) is easily computed from the already known $\psi_{j}^{(m)}$. Then $\vec{\psi}^{\, (m+1)}=A_+^{-1}A_-\vec{\psi}^{\, (m)}$ is determined by inversion of the (essentially) tridiagonal linear system given by the matrix $A_+$; for that purpose, several efficient algorithms are available [PTVF94,Sto99].

As in the case of the TROTTER-based method, the time increment $\delta t$ can be changed adaptively at each time step in order to control the numerical accuracy of the algorithm. Also, as in the TROTTER algorithm, changing the spatial nodes (and thus $\delta x$) in the course of the calculation is much more difficult and prone to cause additional numerical error.

Regarding the initial choice of the parameters $\delta t$ and $\delta x$, it is natural to choose these parameters small enough such that the smallest oscillation of interest (both with respect to time and space) of the system can be resolved. Furthermore, these parameters can be chosen in such a way that the errors induced by time and space discretization are balanced -- this interdependence is made plausible by equation (3.19b) which shows that reducing $\delta x$ leads to similar numerical problems as increasing $\delta t$; more information on this issue can be found in [KW96] and references therein.

next up previous contents
Next: Numerical Evaluation of the Up: Finite Differences Previous: An Algorithm Based on   Contents
Martin Engel 2004-01-01