[Tutor] problems with numdifftools

[Nicolai Heitz]

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

#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.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 
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)

print '\n'
print '--------- Results--------' , '\n'
print 'groundstate positions=', xopt
#print 'groundstate difference=' , xopt[0]-xopt[1]
#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

ps: If you have any questions or need some more information please let 
me know.