Surface Minimum.cpp

                      Probleme de minimisation multidimentionnel
	                            Gilles Texidor - Laurent Chignac
	                            1999 - 2000
	                            DESS MSRO                                           
#include <stdio.h>
#include <iostream.h>
#include <stdlib.h>
#include <fstream.h>
#include "mailleur.h"
#include <math.h>



//********************************************************************
float abs(float a){
  if(a>0) return(a); else return(-a);
}
//********************************************************************
float CalculSurface(maillage xk,float lk,float dk[N-2][N-2]){
  float s;
  maillage xtemp(2.0);
  xtemp.refresh(); 
  for(int i=0;i<N-2;i++){
    for(int j=0;j<N-2;j++){
      xtemp.z[i+1][j+1]=xk.z[i+1][j+1]+lk*dk[i][j];
    }
  }  
  s=xtemp.surface();
  return(s);
}
//**********************************************************************
float minimisation1D(maillage xk,float dk[N-2][N-2]) {
  float lk;
  float A=0.;
  float B=5.;
  float a,b,c,d,e;
  float fc,fd,fe;
  float epsilon=0.01;

  a=A;
  b=B;
  c=0.5*(a+b);
  fc=CalculSurface(xk,c,dk);
    
  while ((b-a)>epsilon){    
    d=0.5*(a+c);
    e=0.5*(c+b);
    fd=CalculSurface(xk,d,dk);
    fe=CalculSurface(xk,e,dk);
    if(fc>fe){
      a=c;
      c=e;
      fc=fe;
    } 
    else { 
      if(fd<fc){
	b=c;
	c=d;
	fc=fd;
      } 
      else {
	a=d;
	b=e;
      }
    }
  } 
  return(c);
}
//********************************************************
//permet de calculer ||gradient||^2
//********************************************************
float norme2_grad(float gradient[N-2][N-2]){
  float somme=0;
  for(int i=0;i<(N-2);i++){
    for(int j=0;j<(N-2);j++){
      somme+=gradient[i][j]*gradient[i][j];
    }
  }
  return(somme);  
}
//********************************************************
void Fletcher_Reeves(){
  float lk;
  float s1,s2,ds=2.;
  float epsilon=0.001;
  float dk[N-2][N-2];
  float bk;
  float norm2gk,norm2gkpun;
  int stop=0;

  maillage xk(2.0);
  xk.refresh_curv();


  xk.calculer_gradient();
  norm2gk=norme2_grad(xk.gradient);
  
  for(int i=0;i<N-2;i++){
    for(int j=0;j<N-2;j++){
  s1=xk.surface();
      dk[i][j]=-xk.gradient[i][j];
    }
  } 
  s1=xk.surface();
  
  while (stop<2){
    lk=minimisation1D(xk,dk);
    for(int i=0;i<N-2;i++){
      for(int j=0;j<N-2;j++){
	xk.z[i+1][j+1]=xk.z[i+1][j+1]+lk*dk[i][j];
      }
    } 
    xk.calculer_gradient();
    norm2gkpun=norme2_grad(xk.gradient);
    bk=norm2gkpun/norm2gk;
    if(stop>=1) bk=0.;
    norm2gk=norm2gkpun;
    for(int i=0;i<N-2;i++){
      for(int j=0;j<N-2;j++){
	dk[i][j]=-xk.gradient[i][j]+bk*dk[i][j];
      }
    }
    s2=xk.surface();
    ds=s1-s2;
    s1=s2;
    if (abs(ds)<epsilon) stop++; else stop=0;
  }
  xk.results();
} 
//********************************************************
//******************************************************** 
main(){
  Fletcher_Reeves();
}