Помогите найти ошибку. Решаю задачу оптимизации, применяю медод стрельбы, м. Ньютона с модификацией Исаева-Сонина(т.е. с дроблением шага) с использванием в условии сходимости нормировки Федоренко
Проблемы: неправильно работает, не выходит из цикла.
Конечно, здесь важное значение на успех программы влияет удачность выбора начальных параметров пристелки 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);
}
обнуляете и опять в цикле... Т.е. 
не имеет смысла.