Код ошибки 27 GNU GSL C ++ для многомерной минимизации: итерация не делает прогресса в направлении решения - PullRequest
Новая
       21

Код ошибки 27 GNU GSL C ++ для многомерной минимизации: итерация не делает прогресса в направлении решения

0 голосов
/ 23 марта 2012

Я пытаюсь найти минимум функции, используя научную библиотеку GNU, пакет «Многомерная минимизация».Я использую метод Бройдена-Флетчера-Голдфарба-Шанно (BFGS), который реализован в функции 

gsl_multimin_fdfminimizer_vector_bfgs2

К сожалению, после первой итерации с 

gsl_multimin_fdfminimizer_iterate

я получаюошибка с кодом 27: итерация не продвигается к решению.

Я пробовал несколько начальных точек и различные комбинации допусков / шагов, но все равно появляется та же ошибка.Не могли бы вы сказать мне, что может быть виновником здесь?

Идея состоит в том, чтобы найти минимум потенциальной энергии для кластера из 13 атомов.Энергия вычисляется в функции my_f, где X, Y и Z - координаты атомов.Используются следующие параметры: 

#define uA 0.1221
#define uXi 1.316
#define uP  8.612
#define uQ  2.516
#define uR0 2.8637

Функция, которую я пытаюсь минимизировать, - это Гупта-Потенциал для кластера из атомов алюминия 13:

// ЗДЕСЬ МОЙ ГУПТА-ПОТЕНЦИАЛ // 

double
      my_f (const gsl_vector *v, void *params)
     {
       double *p = (double *)params;
        int ii,jj;
        double vv = 0;
        double sum1 = 0;
        double sum2 = 0;

          for(ii=1;ii<=NN;ii++){
          X[ii]=gsl_vector_get (v, 3*ii-3);
          Y[ii]=gsl_vector_get (v, 3*ii-2);
          Z[ii]=gsl_vector_get (v, 3*ii-1 );
         }

        for(ii = 1; ii <= NN; ii++){
                sum1 = 0;
                for(jj = 1; jj <= NN; jj++){

                        if(jj!=ii){
                                sum1 += uA*exp(0 - uP*((rr(ii,jj)/uR0) - 1.));
                        }
                }
                sum2 = 0;
                for(jj = 1; jj <= NN; jj++){
                        if(jj!=ii){
                                sum2 += uXi*uXi*exp(0 - 2*uQ*((rr(ii,jj)/uR0) - 1.));
                        }
                }
                vv += sum1 - sqrt(sum2);
        }
        return vv;
 }

// ЗДЕСЬ МОИ ГРАДИЕНТНЫЕ ФУНКЦИИ // 

double dVdr(int ii, int jj){
        double vr;
        double sum1 = 0;
        double sum2 = 0;
        int kk;

        vr = 0. - 2.*uA*(uP/uR0)*exp(0. - uP*((rr(ii,jj)/uR0) - 1.));

        for(kk = 1; kk <= NN; kk++){
                if(kk != ii){
                        sum1 += uXi*uXi*exp(0 - 2*uQ*((rr(ii,kk)/uR0) - 1.));
                }
                if(kk != jj){
                        sum2 += uXi*uXi*exp(0 - 2*uQ*((rr(kk,jj)/uR0) - 1.));
                }
        }

        vr -= 0.5*uXi*uXi*(0. - 2.*uQ/uR0)*exp(0. - 2.*uQ*((rr(ii,jj)/uR0) - 1.))*( (1./sqrt(sum1)) + (1./sqrt(sum2)) );

        return vr;
}

//-----------------------------------------------
 void Force(int jj, double & fx, double & fy, double & fz){
        int mm;
        double dxx, dyy, dzz, r;
        fx = 0;
        fy = 0;
        fz = 0;
   for(mm = 1; mm <= NN; mm++){

    if(mm != jj){
                        dxx = X[mm] - X[jj];
                        dyy = Y[mm] - Y[jj];
                        dzz = Z[mm] - Z[jj];

     r = sqrt(dxx*dxx + dyy*dyy + dzz*dzz);
     fx += dVdr(mm,jj)*(dxx)/r;
     fy += dVdr(mm,jj)*(dyy)/r;
     fz += dVdr(mm,jj)*(dzz)/r;
   }
   }
}
//--------------------------------------------
 void
     my_df (const gsl_vector *v, void *params,
            gsl_vector *df)
     {
       int ii;
       double fx,fy,fz;
       double *p = (double *)params;

          for(ii=1;ii<=NN;ii++){

          X[ii]=gsl_vector_get (v, 3*ii-3);
          Y[ii]=gsl_vector_get (v, 3*ii-2);
          Z[ii]=gsl_vector_get (v, 3*ii-1 );
         }

          for(ii=1;ii<=NN;ii++){
          Force(ii,fx,fy,fz);
          gsl_vector_set(df, 3*ii-3, fx );
          gsl_vector_set(df, 3*ii-2, fy );
          gsl_vector_set(df, 3*ii-1, fz );
         }
     }
//----------------------------------
...