[Tutor] problems with numdifftools

Sat Oct 23 03:30:59 CEST 2010

```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 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.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 to
comment every single step 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

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.66*10**(-27)                    #m is the mass
of a single trapped ion in the chain

#~ print self.w
#~ print self.C
#~ print self.w**2 * self.m

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

#~ m=1
#~ #C=1
#~ w=5

#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(len(x)):
if j >i:
Vx = Vx + C * (abs(x[i]-x[j]))**(-1)    #then we

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)

print '\n'
print '--------- Results--------' , '\n'
print 'groundstate positions=', xopt
#print 'groundstate difference=' , xopt-xopt
#print 'offdiagonal elements = ',
print 'Hessian = ', hessian(xopt)
print 'eigenvalues = ', eigenvalues
print 'normal modes = ', normal_modes

I appreciate any kind of help. Thanks a lot in advance
Nicolai