numerical solution of the time-dependent Schrodinger equation

[Alexander Paul]

Hi!

I'm working on a time-dependent Schrodinger equation with an adiabatic
Hamiltonian H(t) and an initial state psi(0).

Right now I'm solving the equation numerically by discretizing the
time in small intervals of the length delta_t and calculating
(\hbar=1)

psi(t+delta_t) = exp(-i*H(t)*delta_t) \psi(t+delta_t)

I'm sure, there is a faster and more intelligent way to do this. :)

What methods (Runge-Kutta? Monte Carlo?) would you recommend?
Are there some solvers or numerical packages (e.g. for MATLAB) for
this problem?

I would be most grateful if someone could help me on this...

alex.

[jpr0]

If you can calculate \exp(-iH(t)\Delta_t) then I assume you're diagonalizing H(t)? If this is the case then your numerics will run very slowly. I'm not sure about the time-dependent case, but when H is independent of t you can use a Crank-Nicholson scheme. If you get hold of a copy of "numerical recipes in [insert programming language here]" and look in the index for Crank-Nicholson or just Schroedinger's equation, there is a description of the method there.
Basically you approximate the exponent \exp(iA) by


\exp(iA) = [\exp(-iA/2)]^{-1}\exp(iA/2)\approx [1-iA/2]^{-1}[1+iA/2]


so that (t=n\Delta_t)

[tex]
[1+iH]\psi_{n+1} = [1-iH]\psi_n
[/itex]

(I absorbed the factor of 2 into the definition of H). If you're Hamiltonian is a discretized form of the Schroedinger Equation then H should look like a tridiagonal matrix, or at least some sparse matrix which is almost tridiagonal. If this is the case then you can solve for the wavefunction at each time step using fast routines specifically for tridiagonal matrices or whatever special form H takes.

I'm guessing that when H is time dependent there will be a similar method available.

[Chris H. Fleming]

On Jul 6, 12:41 am, Alexander Paul <alexanderpau...@googlemail.com>
wrote:
> Hi!
>
> I'm working on a time-dependent Schrodinger equation with an adiabatic
> Hamiltonian H(t) and an initial state psi(0).
>
> Right now I'm solving the equation numerically by discretizing the
> time in small intervals of the length delta_t and calculating
> (\hbar=1)
>
> psi(t+delta_t) = exp(-i*H(t)*delta_t) \psi(t+delta_t)
>
> I'm sure, there is a faster and more intelligent way to do this. :)
>
> What methods (Runge-Kutta? Monte Carlo?) would you recommend?
>
> Are there some solvers or numerical packages (e.g. for MATLAB) for
> this problem?
>
> I would be most grateful if someone could help me on this...
>
> alex.


That's a pretty good method because it's unitary evolution with a
Hamiltonian slightly different from the actual one.

If H=T(p)+U(q) then you can split Exp(-i H dt) using BCH and a Fast
Fourier Transform routine you can go back and forth between position
space and momentum space to apply separated position operators and
momentum operators trivially. This is called the "split operator
method".

It's also nice if you can work with discrete operators. Say if H =
p^2/2 + q^2/2 + g q^4 , then you can write H in terms of the raising
and lowering operators of the simple harmonic oscillator and truncate.
I'm not sure how efficient this would be for a time-dependent
Hamiltonian as you would have to exponentiate a different H matrix
every time step.

[jpr0]

If you can calculate \exp(-iH(t)\Delta_t) then I assume
you're diagonalizing H(t)? If this is the case then your
numerics will run very slowly. I'm not sure about the time-dependent
case, but when H is independent of t you can
use a Crank-Nicholson scheme. If you get hold of a copy of "numerical
recipes in [insert programming language here]" and look in the index
for Crank-Nicholson or just Schroedinger's equation, there is a
description of the method there.
Basically you approximate the exponent \exp(iA) by

\exp(iA) = [\exp(-iA/2)]^{-1}\exp(iA/2)\approx [1-iA/2]^{-1}[1+iA/2]

so that (t=n\Delta_t)

[tex]
[1+iH]\psi_{n+1} = [1-iH]\psi_n
[/itex]

(I absorbed the factor of 2 into the definition of H). If you're
Hamiltonian is a discretized form of the Schroedinger Equation then H
should look like a tridiagonal matrix, or at least some sparse matrix
which is almost tridiagonal. If this is the case then you can solve for
the wavefunction at each time step using fast routines specifically for
tridiagonal matrices or whatever special form H takes.

I'm guessing that when H is time dependent there will be a similar
method available.

--
jpr0
------------------------------------------------------------------------
jpr0's Profile: http://www.physicsforums.com/forums/member.php?action=getinfo&userid=56745
View this thread: http://www.physicsforums.com/showthread.php?t=176106

[Bill Frensley]

Alexander Paul wrote:
> Hi!
>
> I'm working on a time-dependent Schrodinger equation with an adiabatic
> Hamiltonian H(t) and an initial state psi(0).
>
> Right now I'm solving the equation numerically by discretizing the
> time in small intervals of the length delta_t and calculating
> (\hbar=1)
>
> psi(t+delta_t) = exp(-i*H(t)*delta_t) \psi(t+delta_t)
>
> I'm sure, there is a faster and more intelligent way to do this. :)
>
> What methods (Runge-Kutta? Monte Carlo?) would you recommend?
> Are there some solvers or numerical packages (e.g. for MATLAB) for
> this problem?
>
> I would be most grateful if someone could help me on this...
>
> alex.
>
Look into Cayley's method (otherwise known as Crank-Nicholson for
purely anti-Hermitian operators, which -iH is). It effectively
approximates:
e^{-iH 2dt} \approx (1 - iH dt) / (1 + iH dt) ,
which is exactly unitary for Hermitian H. The nature of the error
in this approximation is that the rate of phase advance for higher
energy states is less than the exact solution, but there is usually
negligible amplitude in the higher-energy states in most problems.

This does require solution of equations of the form y = Hx. If
the matrix representation of your Hamiltonian is such that this
can be done efficiently (in particular, if H is tridiagonal or
a close relative) Cayley's method is usually far and away the most
effective technique. Note that if you have a one-dimensional
system, a finite-difference approximation of the (second-order
differential) Hamiltonian will produce a tridiagonal matrix.

- Bill Frensley

[Chris H. Fleming]

On Jul 6, 8:09 pm, "Chris H. Fleming" <chris_h_flem...@yahoo.com>
wrote:
> On Jul 6, 12:41 am, Alexander Paul <alexanderpau...@googlemail.com>
> wrote:
>
>
>
> > Hi!
>
> > I'm working on a time-dependent Schrodinger equation with an adiabatic
> > Hamiltonian H(t) and an initial state psi(0).
>
> > Right now I'm solving the equation numerically by discretizing the
> > time in small intervals of the length delta_t and calculating
> > (\hbar=1)
>
> > psi(t+delta_t) = exp(-i*H(t)*delta_t) \psi(t+delta_t)
>
> > I'm sure, there is a faster and more intelligent way to do this. :)
>
> > What methods (Runge-Kutta? Monte Carlo?) would you recommend?
>
> > Are there some solvers or numerical packages (e.g. for MATLAB) for
> > this problem?
>
> > I would be most grateful if someone could help me on this...
>
> > alex.
>
> That's a pretty good method because it's unitary evolution with a
> Hamiltonian slightly different from the actual one.
>
> If H=T(p)+U(q) then you can split Exp(-i H dt) using BCH and a Fast
> Fourier Transform routine you can go back and forth between position
> space and momentum space to apply separated position operators and
> momentum operators trivially. This is called the "split operator
> method".
>
> It's also nice if you can work with discrete operators. Say if H =
> p^2/2 + q^2/2 + g q^4 , then you can write H in terms of the raising
> and lowering operators of the simple harmonic oscillator and truncate.
> I'm not sure how efficient this would be for a time-dependent
> Hamiltonian as you would have to exponentiate a different H matrix
> every time step.


I just thought that if the time dependence can be factored out of the
operators in your Hamiltonian,
for instance
H(t) = H0 + f(t) H1
then you should definitely use BCH on Exp(-i H dt) and calculate Exp(-
i H0) and Exp(-i H1) a priori.