Hi Ondrej,
Am 31.01.2010 um 08:43 schrieb Ondrej Certik:
...
I run you example and got a figure (attached). Nice!
So looking at the source code, here is the hamiltonian:
hamilton_operator = sum((ba.Ket(basis_states[i]) *
ba.Bra(basis_states[(i + 1) % len(basis_states)]) for i in
xrange(len(basis_states))))
#make hermitian
hamilton_operator += hamilton_operator.adjoint()
so is it:
H = sum_i |i><i+1|
Yes exactly, except that in a second step I also add the hermitian
conjugate, to make the hamilton operator hermitian.
H = sum_i ( |i><i+1| + |i+1><i| )
I don't know of any particular system where this hamiltonian applies,
it was just meant as a demonstration.
A more "realistic" Hamiltonoperator would be for example a short Ising
chain.
Define:
def sigma_z(site_index):
return ba.K('e', site_index)* ba.B('e', site_index) - ba.K('g',
site_index)* ba.B('g', site_index)
def sigma_x(site_index):
return ba.K('e', site_index)* ba.B('g', site_index) + ba.K('g',
site_index)* ba.B('e', site_index)
In the single site basis (|e>, |g>) these operators take on the form
of the pauli matrices sigma_z and sigma_x.
These operators and products or sums of them are hermitian, so they
are good building blocks for spin chain systems.
#We can then get the Hamiltonian. In code:
J = 1.
Omega = .01
H = (J / 4) * sum(( sigma_z(site_index) * sigma_z( (site_index+1) % 4)
for site_index in range(4)) + Omega * sum((sigma_x(site_index) for
site_index in range(4))
This is now the hamilton operator for a spin chain with spin spin
interaction strength J and coupling to a transverse magnetic field
with fieldstrength Omega.
It has also other interpretations, but this is to be expected, since
basically any finite system consisting of 2d-subsystems can very
easily be expressed in terms of
the (local) pauli operators
So then I start in an initial state already in coordinate
representation (1, 0, 0, 0...0) which corresponds to the first basis
state 'gggg'.
To make this more explicit one could substitute the lines in the code
"initial_state = ..."
with this:
initial_state_ket = ba.Ket('gggg')
# project the initial state onto the basis, by inserting the identity
on the full space \sum_i |i><i| == identity == 1
# |psi> = 1* |psi> = (\sum_i |i><i| ) |psi> = \sum_i |i> (< i |
psi>) = \sum_i c_i | i >, where c_i = < i | psi>
initial_state = np.array([Bra(b_i) * initial_state_ket for b_i in
basis_states ]) #calculate the c_i
# so if |psi> = |gggg> => c_0 = 1, c_k = 0 for k > 0 (if our labels
go from 0 to dim - 1)
# so the c_i are basically the basis_coefficients at time zero t = 0
? Then you take some initial state and get it evolved? Then you plot
the value of this observable:
Right, so the system state's time evolution is determined by the
Schroedinger equation.
H | psi(t) > = i d/dt |psi(t)>
By diagonalizing the Hamiltonian we get its eigenbasis:
H |n(t)> = omega_n |n(t)>
Their time evolution is quite simple as can be seen by plugging it
into the schroedinger equation:
omega_n |n(t)> = H |n(t)> = i d/dt |n(t)> ===> d/dt |n> = - i*omega_n |
n(t)> === > |n(t)> = exp(- i*omega_n * t) |n(t=0)>
So all we have to do is express our state |psi(t)> in terms of the |
n(t)>.
Since H is hermitian, the time evolution of the system is unitary,
i.e. it is enough to project |psi(t = 0)> onto the |n(t=0)>.
The coefficients given by this must then remain constant:
<n(t) | psi(t)> = <n(t = 0) | psi(t = 0)> =: a_n
because \sum_n |n(t)><n(t)| == \sum_n |n(t)><n(t)| == 1
==> |psi(t)> = \sum_n (|n(t)><n(t)| )|psi(t)> = \sum_n |n(t)> (<n(t) |
psi(t)> ) = \sum_n |n(t)> a_n = \sum_n a_n * exp(- i*omega_n * t) * |
n(t=0)>
No we can ask, what is then the time dependance of our
basis_coefficients c_i with respect to the basis_states defined at the
beginning:
c_i(t) := < i | psi(t)> = \sum_n a_n * exp(- i*omega_n * t) * < i |
n(t=0)>
This last result corresponds exactly to the following expression,
"coefficients" corresponds to (a_n) and "eigenvals" to
"omega_n" (they are both arrays).
Moreover < i | n(t=0)> is simply the i-th component of the n-th
eigenvector and this corresponds to eigenvectors[i,n] (since the
eigenvectors are returned as columns of a matrix by scipy.linalg.eigh)
The np.dot effectively sums over n. The list comprehension calculates
this for every time step. The result is a list of vectors, each vector
representing all c_i at a particular time.
[np.dot(eigenvectors, coefficients*np.exp(-1j * eigenvals * t)) for t
in times]
my_observable = ba.KetBra('gggg','gggg') #ground state projector
The expectation value of this operator gives (for every time step) the
probability of finding the system in the state |gggg>.
This can be seen from:
<my_observable(t)> := <psi(t)| (|gggg><gggg|) |psi(t)> = | <gggg|
psi(t)> |^2 = |c_0(t)|^2
I will extend this demo a little bit and add more comments and repost
it later :)
? If you could say a few words what exactly it is computing, it'd be
cool.
We can then make it an example in sympy, once we integrate it with the
code that we already have.
Ondrej
--
You received this message because you are subscribed to the Google
Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com
.
For more options, visit this group at http://groups.google.com/group/sympy?hl=en
.
<image.png>
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to
sympy+unsubscr...@googlegroups.com.
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
