м. Ньютона с модиф-ей Исаева-Сонина с нормировкой Федоренко

[arman]

Помогите найти ошибку. Решаю задачу оптимизации, применяю медод стрельбы, м. Ньютона с модификацией Исаева-Сонина(т.е. с дроблением шага) с использванием в условии сходимости нормировки Федоренко
Проблемы: неправильно работает, не выходит из цикла.
Конечно, здесь важное значение на успех программы влияет удачность выбора начальных параметров пристелки alpha0[i].
вот неполный листинг программы:

void FCN(int , double , double *, double *);
void Newton();
void Gauss(double **, double *, double *);
void OutPut1(double, double *);
double norma(int n, double *F);
double norma(int n, double *F){
    int i;
    double res=0;
    
    for(i=0; i<n; i++)   res+=pow(F[i],2);
    return sqrt(res);
}
    
double norma_fedor(int n, double *F, double *kappa);    
double norma_fedor(int n, double *F, double *kappa){    
    int i;
    double res=0;
    for(i=0; i<n; i++)
       {            
            res+=pow(F[i],2)/kappa[i];
       }
    return sqrt(res);
}
double calculate(double t_0, double *alpha, double *delta, double *delta_alpha, double *F, double *kappa);
double calculate(double t_0, double *alpha, double *delta, double *delta_alpha, double *F, double *kappa)
{
   int i,j;
   double   *Y = (double*)malloc(15*sizeof(double));
   double DELTA_B[6][6],DELTA_A[6][6];
   double **J;
   J=(double**)malloc(6*sizeof(void*));
   J[0]=(double*)malloc(36*sizeof(double));
   for(i=1; i<6; i++) J[i]=J[i-1]+6;
        
   func(t_0_earth, t_0, &x_earth, &y_earth, &z_earth, &velocity_x_earth,&velocity_y_earth,&velocity_z_earth, &r_e, a_e, e_e, I_e, OMEGA_e, w_e, M_0_e);
   func(t_0_ast, alpha[0], &x_ast, &y_ast, &z_ast, &velocity_x_ast,&velocity_y_ast, &velocity_z_ast, &r, a, e, I, OMEGA, w, M_0);
   nach_uslovia(14, Y, alpha);
   DOPRI8(14, FCN, t_0, alpha[0], Y, 1e-12, 50, 1);
   for(i=0; i<6; i++)
   {
      for(j=0;j<6;j++)
         {
            DELTA_B[i][j]=Y[i+1];
         }
   }      
   //printf("Pravaya chast sistemy ur-ii, v-r F: \n");
//F-вектор невязки!!!

   F[0]=Y[1]-x_ast; F[1]=Y[2]-y_ast; F[2]=Y[3]-z_ast; F[3]=Y[4]; F[4]=Y[5]; F[5]=Y[6];

   //izmneyaem alfa dlya vychislenia F(...,alpha[i]+delta[i], ...)
   for(j=0;j<6;j++)
   {
      alpha[j]=alpha[j]+delta[j];
      nach_uslovia(14, Y, alpha);   
      DOPRI8(14, FCN, t_0, alpha[0], Y, 1e-12, 50, 1);
      for(i=0; i<6; i++)
         {
            DELTA_A[i][j]=Y[i+1];
         }
      alpha[j]=alpha[j]-delta[j]; //vozvrashaem prezgnee znachenie alpha[j]
   }
   printf("J=Matrica Jacobi\n");
   for(i=0;i<6;i++){
      for(j=0; j<6; j++)
         {
            J[i][j]=((DELTA_A[i][j]-DELTA_B[j][i])/delta[j]);
      //      printf("%le\t", J[i][j]);
         }
      //   printf("\n");   
      }
// NORMIROVKA FEDORENKO!!!!
   for(i=0; i<6; i++){
         kappa[i]=0;
         for(int j=0; j<6; j++)
         {   
            kappa[i]=+pow(J[i][j],2);
         }
      }
   Gauss(J,F,delta_alpha); //решаем сист. лин. ур-й метода Гаусса J*(delta_alpha)=-F
                                            //найдем векотр delta_alpha[i]-элементы обратной матрицы к матр. Якоби J
   double norma_fedorenko;
   norma_fedorenko=norma_fedor(6, F, kappa);
   free(J[0]);
   free(J); 
   return norma_fedorenko;
}
void Newton()
{
   double t_0=55703.251575;   //t_0-nach. moment[MJD]
      int i,j;
   double eps_newton=1.0e-8;   //Tochnost' v metode Newtona
   double gamma;
   double NORMA=1.e-2; 
   double norma_f, err;
   double norma_f_old;
   double norma_f_new;
   double delta[6],kappa[6],alpha[6], alpha_next[6], F[6],delta_alpha[6];
   double F_old[6], F_next[6], F_zap[6], zap[6];
   bool fl;     
   for(i=0;i<6;i++)
   {
      delta[i]=1.0e-6;
   }
   alpha[0]=55753.251575;   //=T   //perelet za 170(libo 40 dnei) dnei//nach. priblizhenie
   alpha[1]=-10*velocity_y_earth;
   alpha[2]=-10*velocity_z_earth;
   alpha[3]=-0.00001;
   alpha[4]=-0.000001;
   alpha[5]=-0.0001;
   printf("==============NEWTON==================\n");
   int q=0; //q-число итераций
gamma=1.0;
do
{   
   q++;
   printf("================q=%d============\n",q);
while(fl==0)
{
   norma_f_old=calculate(t_0, alpha, delta, delta_alpha, F, kappa);
   for(i=0; i<6; i++) 
   {   zap[i]=alpha[i];
      F_zap[i]=F[i];
      alpha_next[i]=alpha[i]-gamma*delta_alpha[i];
   }   
   norma_f_new=calculate(t_0, alpha_next, delta, delta_alpha, F, kappa);
   if(norma_f_old<norma_f_new)
   {
      alpha[i]=zap[i];
      F[i]=F_zap[i];
      gamma=0.5*gamma;
      fl=1;
   }
   printf("gamma=%.20le\n", gamma);
   NORMA=norma(6, F);
   printf("NORMA=%le\n", NORMA);
   printf("norma_f_new=%le\n", norma_f_new);
   fl=0;
}
}while(NORMA>eps_newton && ++q<MAX_IT);
}
void FCN(int N, double T, double *Y, double *F) 
{
   N=14;
   double rot=P_opt/(sqrt(Y[4]*Y[4]+Y[5]*Y[5]+Y[6]*Y[6])*Y[7]);   //P_opt/((ro)*m)
   double val=(niu_s/pow(r,5))*(Y[4]*Y[1]+Y[5]*Y[2]+Y[6]*Y[3]);
   
   F[1]=Y[8];
   F[2]=Y[9];
   F[3]=Y[10];
     F[4]=-Y[11];
     F[5]=-Y[12];
     F[6]=-Y[13];
 
     F[7]=P_max/c;//-(sqrt(P_x*P_x+P_y*P_y+P_z*P_z))/c;
  
     F[8]=-niu_s/pow(r,3)*Y[1]+(Y[4]*rot);                  //-niu_s/pow(r,3)*Y[0]+P_x/Y[6];
     F[9]=-niu_s/pow(r,3)*Y[2]+(Y[5]*rot);                  //-niu_s/pow(r,3)*Y[1]+P_y/Y[6];
     F[10]=-niu_s/pow(r,3)*Y[3]+(Y[6]*rot);                  //-niu_s/pow(r,3)*Y[2]+P_z/Y[6];
  
     F[11]=-niu_s/pow(r,3)*Y[4]-3*Y[1]*val;
     F[12]=-niu_s/pow(r,3)*Y[5]-3*Y[2]*val;
     F[13]=-niu_s/pow(r,3)*Y[6]-3*Y[3]*val;

     F[14]=(1/Y[7]) *(rot)*(Y[4]*Y[4]+Y[5]*Y[5]+Y[6]*Y[6]);//(1/pow(Y[6],2))*(Y[3]*P_x+Y[4]*P_y+Y[5]*P_z);
}
void nach_uslovia(int N, double *Y, double *alpha)
{
   Y[1] = x_earth;                   
   Y[2] = y_earth;                
   Y[3] = z_earth;    
   Y[4]=-10*velocity_x_earth;
   Y[5]=alpha[1];
   Y[6]=alpha[2];
   Y[7]=1000.;
   Y[8]=velocity_x_earth;               //u_earth(t_0);
   Y[9]=velocity_y_earth;               //v_earth(t_0);
   Y[10]=velocity_z_earth;               //w_earth(t_0);
   Y[11]=alpha[3];
   Y[12]=alpha[4];
   Y[13]=alpha[5];
   Y[14]=-pow(10,20);//-pow(10,20);
}

[p51x]

arman в сообщении #286801 писал(а):

   if(norma_f_old<norma_f_new)   {      alpha[i]=zap[i];      F[i]=F_zap[i];      gamma=0.5*gamma;      fl=1;   }

А дальше вы $fl$ обнуляете и опять в цикле... Т.е. $fl=1$ не имеет смысла.

[arman]

p51x
А что предлагаешь? убрать fl=1;

[p51x]

Тогда будет вечный цикл....

Если нужны вывод и вычисления после if'a, тогда подумайте над алгоритмом. Иначе break или return.

[arman]

Может кто напишет правильный алгоритм модифицированного м. Ньютона с нормировкой федоренко, а то я че-то запарился)