[Tutor] problems with numdifftools

[Steven D'Aprano]

On Sat, 23 Oct 2010 12:30:59 pm Nicolai Heitz wrote:
> 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'd try the numpy mailing list:

http://www.scipy.org/Mailing_Lists


[...]
> 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!

Yes, I think this is definitely a question for the numpy community :)


> 5) In principal Python can divide those numbers (2.3e-28)/(5.6e-6)**3

I get half the result you suggest above:

>>> (2.3e-28)/(5.6e-6)**3
1.3096756559766764e-12

which agrees with my HP-48GX, and the result after simplifying by hand 
to 2.3e-10/175.616. So I think that is the correct result, not 2.6e-12.

> 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:

Comments are good, but too many comments hurt readability. Don't feel 
that you need to comment every single line. Comments like:

x = x+1  # increment x by 1

are just noise. But this is a good comment:

x = x+1  # add correction for wind-chill factor.


> self.m = 39.96*1.66*10**(-27)                    #m is the mass
> of a single trapped ion in the chain

I wonder why you calculate m with an explicit exponentiation operator, 
instead of a floating point exponent?

The difference is very small, so it probably doesn't matter:

>>> 39.96*1.66*10**(-27)
6.6333600000000006e-26
>>> 39.96*1.66e-27
6.6333599999999994e-26




>          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

You can simplify that nested loop and addition:

for i in range(len(x)):
    for j in range(i+1, len(x)):
        Vx += C * (abs(x[i]-x[j]))**(-1)  # add the coulomb interaction


The function range takes up to three integer arguments:

range(start, end, step)

where start and step are optional, defaulting to 0 and 1 respectively. 
By starting the inner loop at i+1, you guarantee that j is always > i 
and therefore you can skip the test.

The other change is that instead of 

Vx = Vx + something

you can write:

Vx += something





-- 
Steven D'Aprano