On Wed, Oct 27, 2010 at 12:59 AM, Pauli Virtanen wrote:

Tue, 26 Oct 2010 14:24:39 -0700, Nicolai Heitz wrote:

...

...

 http://mail.scipy.org/mailman/listinfo/scipy-user

I contacted them already but they didn't responded so far and I was forwarded to that list which was supposed to be more appropriated.

I think you are thinking here about some other list -- scipy-user is the correct place for this discussion (and I don't remember seeing your mail there).

[clip]

...

1) Can I make it run/fix it, so that it is also going to work for the SI scaling?

Based on a brief look, it seems that uniform scaling will not help you, as you have two very different length scales in the problem,

       1/sqrt(m w^2) >> C

If you go to CM+relative coordinates you might be able to scale them separately, but that's fiddly and might not work for larger N.

In your problem, things go wrong when the ratio between the length scales approaches 1e-15 which happens to be the machine epsilon. This implies that the algorithm runs into some problems caused by the finite precision of floating-point numbers.

What exactly goes wrong and how to fix it, no idea --- I didn't look into how Numdifftools is implemented.

...

2) How can I be sure that increasing the number of ions or adding a somehow more complicated term to the potential energy is not causing the same problems even in natural units?

3) In which range is numdifftools working properly.

That depends on the algorithm and the problem. Personally, I wouldn't trust numerical differentiation if the problem has significantly different length scales, it is important to capture all of them accurately, and it is not clear how to scale them to the same size. Writing ND software that works as expected all the time is probably not possible even in theory.

Numerical differentiation is not the only game in the town. I'd look into automatic differentiation (AD) -- there are libraries available for Python also for that, and it is numerically stable.

E.g.

http://en.wikipedia.org/wiki/Automatic_differentiation#Software

has a list of Python libraries. I don't know which of them would be the best ones, though.

they all have their pro's and con's. Being (co-)author of some of these tools, my personal and very biased advice is: if you are on Linux, I would go for PYADOLC. it provides bindings to a feature-rich and well-tested C++ library. However, the installation is a little tricker than a "setup.py build" since you will need to compile ADOL-C and get Boost::Python to work. PYADOLC can also differentiate much more complicated code than your example in a relatively efficient manner. For your example the code looks like: --------------------- code --------------------------- .... c=classicalHamiltonian() xopt = optimize.fmin(c.potential, c.initialposition(), xtol = 1e-10) import adolc; import numpy # trace the computation adolc.trace_on(0) x = adolc.adouble(c.initialposition()) adolc.independent(x) y = c.potential(x) adolc.dependent(y) adolc.trace_off() hessian = adolc.hessian(0, xopt) eigenvalues = numpy.linalg.eigh(hessian)[0] normal_modes = c.normal_modes(eigenvalues) print 'hessian=\n',hessian print 'eigenvalues=\n',eigenvalues print 'normal_modes=\n',normal_modes --------------------- code --------------------------- and you get as output Optimization terminated successfully. Current function value: 0.000000 Iterations: 81 Function evaluations: 153 hessian= [[ 5.23748399e-12 -2.61873843e-12] [ -2.61873843e-12 5.23748399e-12]] eigenvalues= [ 2.61874556e-12 7.85622242e-12] normal_modes= [ 6283185.30717959 10882786.30440101] Also, you should use an eigenvalue solver for symmetric matrices, e.g. numpy.linalg.eigh. regards, Sebastian

-- Pauli Virtanen

_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion