Алгоритм метода Нелдера-Мида

[Imperator]

Нужно запрограммировать на C\C++ метод деформированного многогранника и прогнать его на некоторых "неэлементарных" функциях.
Метод запрограммировать совсем несложно. Но проблема в том, что, видимо, где-то что-то не учел. В общем результаты при минимизации простой функции 2-х переменных отличаются для различных начальных симплексов. А если брать начальный треугольник далеко от точки минимума, то сходимости даже для этой простой функции нет. Итог: где-то ошибся, но не могу понять, где именно.

Вопрос заключается в том, где можно найти готовый код этого метода на любом языке, чтобы можно было понять ошибку и обобщить алгоритм для многих переменных?

[Yuri Gendelman]

Код на Фортране:
Химмельблау Д.
Прикладное нелинейное программирование.
Мир, 1975.

[Imperator]

Yuri Gendelman, спасибо Вам.

Вот код на Turbo C. Может кому пригодится. Обобщить его не составит труда.

#include <iostream.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
#define n  2

double f (double x, double y) {
   //return 4 * (x - 5) * (x - 5) + (y - 6) * (y - 6);
   return (1 - x) * (1 - x) + 100 * (y - x * x) * (y - x * x);
   //return (x * x + y - 11) * (x * x + y - 11) + (x + y * y - 7) * (x + y * y - 7);
}

void main () {
   clrscr();
   const double eps = 0.000000005;
   double x[n + 5][n];
   double s;

   x[0][0] = 0;  x[0][1] = 0;
   x[1][0] = 10; x[1][1] = 11;
   x[2][0] = 8; x[2][1] = 11;


   int h, l, g;
   long int kol = 0;

   cout<<x[0][0]<<"      "<<x[0][1]<<endl;
   cout<<x[1][0]<<"      "<<x[1][1]<<endl;
   cout<<x[2][0]<<"      "<<x[2][1]<<endl;
   cout<<"_________________________________"<<endl;

   do {
   kol++;

   h = 0;
   l = 0;
   g = 0;
   double max = f(x[0][0], x[0][1]);
   double min = f(x[0][0], x[0][1]);
   for (int k = 1; k <= n; k ++) {
   if (f(x[k][0], x[k][1]) < min) {
   min = f(x[k][0], x[k][1]);
   l = k;
   }
   if (f(x[k][0], x[k][1]) > max) {
   max = f(x[k][0], x[k][1]);
   h = k;
   }
   }

   max = f(x[0][0], x[0][1]);

   for (k = 1; k <= n; k ++) {
   if ((f(x[k][0], x[k][1]) > max) && (k != h)) {
   max = f(x[k][0], x[k][1]);
   g = k;
   }
   }

   x[n + 1][0] = (x[0][0] + x[1][0] + x[2][0] - x[h][0]) / n;
   x[n + 1][1] = (x[0][1] + x[1][1] + x[2][1] - x[h][1]) / n;


   for (int i = 0; i < n; i++)
   x[n + 2][i] = 2 * x[n + 1][i] - x[h][i];

   if (f(x[n + 2][0], x[n + 2][1]) < f(x[l][0], x[l][1])) {
   x[n + 3][0] = x[n + 1][0] + 2 * (x[n + 2][0] - x[n + 1][0]);
   x[n + 3][1] = x[n + 1][1] + 2 * (x[n + 2][1] - x[n + 1][1]);

   if (f(x[n + 3][0], x[n + 3][1]) < f(x[l][0], x[l][1])) {
   x[h][0] = x[n + 3][0];
   x[h][1] = x[n + 3][1];
   }
   else {
   x[h][0] = x[n + 2][0];
   x[h][1] = x[n + 2][1];
   }
   }

   else
   if (f(x[n + 2][0], x[n + 2][1]) > f(x[g][0], x[g][1])) {
   if (!(f(x[n + 2][0], x[n + 2][1]) > f(x[h][0], x[h][1]))) {
   x[h][0] = x[n + 2][0];
   x[h][1] = x[n + 2][1];
   }
   x[n + 4][0] = x[n + 1][0] + 0.5 * (x[h][0] - x[n + 1][0]);
   x[n + 4][1] = x[n + 1][1] + 0.5 * (x[h][1] - x[n + 1][1]);
   if (!(f(x[n + 4][0], x[n + 4][1]) > f(x[h][0], x[h][1]))) {
   x[h][0] = x[n + 4][0];
   x[h][1] = x[n + 4][1];
   }
   else {
   for (k = 0; k < n; k++) {
   x[k][0] = x[k][0] + 0.5 * (x[k][0] - x[l][0]);
   x[k][1] = x[k][1] + 0.5 * (x[k][1] - x[l][1]);
   }
   }
   }
   else {
   x[h][0] = x[n + 2][0];
   x[h][1] = x[n + 2][1];
   }

   float s1 = 0, s2 = 0;
   for (k = 0; k < n + 1; k++) {
   s1 += f(x[k][0], x[k][1]);
   s2 += f(x[k][0], x[k][1]) * f(x[k][0], x[k][1]);
   }
   s = s2 - s1 * s1 / (n + 1);
   s /= (n + 1);

   cout<<x[0][0]<<"      "<<x[0][1]<<endl;
   cout<<x[1][0]<<"      "<<x[1][1]<<endl;
   cout<<x[2][0]<<"      "<<x[2][1]<<endl;
   min = f(x[0][0], x[0][1]);
   for (k = 1; k <= n; k++) {
   if (f(x[k][0], x[k][1]) < min) {
   min = f(x[k][0], x[k][1]);
   l = k;
   }
   }
   cout<<"f = " <<f(x[l][0], x[l][1])<<endl;
   cout<<"___________________"<<endl;


   } while (s > eps);

   cout << "x = " << x[n + 1][0] << endl;
   cout << "y = " << x[n + 1][1] << endl;
   cout << "minF = " << f(x[n + 1][0], x[n + 1][1]) << endl;
   cout << "Количество итераций = " << kol;

}