Amoeba simplex


wies
12-12-2005, 12:05 PM
(I may have prematurely posted a non-finished draft of this message before. Apparently hit the wrong button. Sorry.)

I try to use amoeba.c to find the minumum of Rosenbrocks test function f(x,y)=100(x-y)^2+(1-x)^2. The minimum is given by f(1,1)=0.

This is a function in 2 dimensions, so 3 starting points are needed. I use the following:
p3 = (1,a)
p1 = p3 + e1
p2 = p3 + e2
where e1 and e2 are the 2-d unit vectors.

With e.g. a=6 there are no problems, but with a=2 or a=3 the program gets stuck, the Borland compiler gives a message "Thread stopped" and I have to reset the whole process.

Upon inspection I find that just before crashing the following two things happen:
1. one of the three simplex points takes exactly the values 1,1 (tested using the == operator)
2. and just below that, two of the three simplex points will coïncide, so the simplex is no longer n+1 dimensional (also tested using ==).

Does anybody have similar experiences? And perhaps a remedy? If not, then of course I could try and modify the routine myself, but I am only a so-so programmer, and it seems such a shame to break into that beautifully written code!

Wies.

Saul Teukolsky
12-12-2005, 09:09 PM
Dear Wies,

I'm unable to reproduce your problem using the gcc compiler. In the end position, y[ihi]=y[ilo]=0, so rtol is exactly zero and the routine exits. The 3 vectors in the simplex are indeed all (1,1).

Saul Teukolsky

wies
12-13-2005, 02:47 AM
Dear professor Teukolsky,

thank you very much for your reply. After posting the mail, I found out that I have to replace the line
rtol=2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo]))
by
rtol=2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo]) + TINY).

I had checked this site for bug reports but did not find any on the amoeba routine. Perhaps this one is not mentioned because my copy of NR, although being second edition, qualifies as version 1.99 ??

It is probably time for me to buy a new set of book&files.

Thanks again for taking the trouble to answer in person to questions on this forum. It must cost you quite some time!

Wies Akkermans

Saul Teukolsky
12-13-2005, 07:44 AM
Dear Wies,

Glad you found the problem. Bugfixes earlier than those for version 2.08 are not listed on the forum (sorry about that). Don't buy a new book - you can look at the relevant section at www.nr.com for the latest C or Fortran versions and check anything you're suspicious about.

Best wishes,
Saul Teukolsky