Hey,
I am not sure if you are the right persons to contact but if not I would
appreciate a short notice and maybe an address where I can find help. I already
posted this
this message in an other python mailing list and they forwarded me to this list
and told me that I might find help here.
I am using Python 2.6 and the following packages:
1) numpy
2) scipy
3) numdifftools
I am a physicist and started to use Python 2-3 weeks ago. I want to use
Python to calculate the eigenvalues of the Hamiltonian given in the
code. Therefore I have to do the following steps:
1) define the Hamiltonian or the potential respectively (--> function:
potential)
2) minimize the potential ( I am using scipy.optimize.fmin to calculate
the minimum of the potential)
(2b) Expand the Hamiltonian around the minimum position. This step is
not done in the code because it is only necessary for the purpose that
one understand the physical background and why one have to do step 3)
3) Calculate the Hessian matrix of the potential, that means calculate
the second derivatives of the potential at the point of the minimum
position (--> numdifftools.Hessian)
4) Calculate the eigenvalues of the Hessian (-->scipy.linalg.eigvals)
The problem can be solved analytically except of the calculation of the
minima in step 2:
Now I have the following problem:
The Hessian matrix evaluated at the minimum position is wrong.
What I checked so far:
1) The potential seems to be calculated correctly and the minimum
position for two ions (N=2) is fine.
2) The Hessian matrix calculation works fine for several other functions.
3) The Hessian matrix is correct when I set the Coulomb interaction to
zero (C=0)
4) The Hessian matrix is correct when I set C to some simple arbitrary
integer numbers like 1, 5 ,8 ....
5) The Hesse matrix gives the correct number of the first part of the
diagonal elements even with C=2.3e-28 ! But in that case the offdiagonal
elements and the part of the diagonal element which refers to the second
derivative in the Coulomb term is wrong! The offdiagonal term should be
C/(abs(x_1-x_2))**3 = 2.6e-12 but it is 3.4e-24!
5) In principal Python can divide those numbers (2*2.3e-28)/(5.6e-6)**3
6) I played around with the accurateness of the calculation of the
Hessian by varying the arguments numTerms, method and metOrder but that
didn't change anything in my case.
My source code is attached. Please keep in mind that I just started
programing and Python is my first Language. I tried to do my best with the
comments to make the code readable. Here it comes:
# import a tool to use / as a symbol for normal division
from __future__ import division
#import system data
import sys, os
#Loading the required packages
import scipy as sp
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import numdifftools as nd
# and subpackages
from scipy import *
from scipy import linalg, optimize, constants
#-----------------------------------------------------------------------------------------
# Hamiltonian H=sum_i(p_i^2/(2m)+ 1/2 * m * w^2 x_i^2)+
sum_(i!=j)(a/|x_i-x_j|)
#-----------------------------------------------------------------------------------------
class classicalHamiltonian:
def __init__(self):
self.N = 2
#N is a scalar, it's the number of ions in the chain
f = 1000000
#f is a scalar, it's the trap frequency
self.w = 2*pi*f
#w is a scalar, it's the angular velocity corresponding to the trap frequency
self.C = (4*pi*constants.epsilon_0)**(-1)*constants.e**2
#C is a scalar, it's the Coulomb constant times the electronic charge in SI
self.m = 39.96*1.66e-27
#m is the mass of a single trapped ion in the chain
def potential(self, positionvector):
#Defines the potential that is going to be minimized
x= positionvector
#x is an 1-d array (vector) of lenght N that contains the positions of the N
ions
w= self.w
C= self.C
m= self.m
#First we consider the potential of the harmonic osszilator
Vx = 0.5 * m * (w**2) * sum(x**2)
for i in range(len(x)):
for j in range(i+1, len(x)):
Vx += C * (abs(x[i]-x[j]))**(-1) #then
we add the coulomb interaction
return Vx
def initialposition(self): #Defines the initial position
as an estimate for the minimize process
N= self.N
x_0 = r_[-(N-1)/2:(N-1)/2:N*1j]
return x_0
def normal_modes(self, eigenvalues): #the computed eigenvalues of
the matrix Vx are of the form (normal_modes)^2*m.
m = self.m
normal_modes = sqrt(eigenvalues/m)
return normal_modes
#C=(4*pi*constants.epsilon_0)**(-1)*constants.e**2
c=classicalHamiltonian()
#print c.potential(array([-0.5, 0.5]))
xopt = optimize.fmin(c.potential, c.initialposition(), xtol = 1e-10)
hessian = nd.Hessian(c.potential)
eigenvalues = linalg.eigvals(hessian(xopt))
normal_modes = c.normal_modes(eigenvalues)
I appreciate any kind of help. Thanks a lot in advance
Nicolai
ps: If you have any questions or need some more information please let
me know.
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
