Accueil | Cast3M

* fichier :  identifi.dgibi
************************************************************************
************************************************************************
*
* soit le polynome du 3eme degre y=a*x*x*x + b*x*x + c*x + d
*
* on suppose connue la valeur  d'une fonction G pour les
* abscisses x= 1, 2, 3, 4 et 5 .On recherche les coeffiecients a b c d
* de la fonction Y (polynome du troisieme degré) qui ajuste au mieux la
* fonction G sur les points connus.
 
 
* poly est une procédure calculant la fonction Y pour des valeurs de
* a b c d données en argument et pour les valeurs expérimentales de x
 
debp poly a*flottant b*flottant c*flottant  d*flottant;
  x=1 ;
  f1=(a*x*x*x) + (x*x*b) + ( c * x ) + d  ;
  x=2;
  f2=(a*x*x*x) + (x*x*b) + ( c * x ) + d  ;
  x=3;
  f3=(a*x*x*x) + (x*x*b) + ( c * x ) + d  ;
  x=4;
  f4=(a*x*x*x) + (x*x*b) + ( c * x ) + d  ;
  x=5;
  f5=(a*x*x*x) + (x*x*b) + ( c * x ) + d  ;
  ll=prog f1 f2 f3 f4 f5;
finpro ll;
 
* Nous choisissons la base connue de G corerspondant à un polynome du
* 3eme degré avec a=2 b=3 c=4 d=5 . Nous devons donc retrouver ces
* valeurs pour a b c d.
fG = poly 2. 3. 4.  5.;
 
* on cherche par l'opérateur MOCA les valeurs de a b c d
* on suppose un point de départ a=1 b=-1. c=10; d=8.
a=1.;
b=-1.;
c=10.;
d = 8.;
fdep = poly a b c d ;
* calcul de la fonction y pour ces paramètres et pour les abscisses
depa = prog a b c d ;
* calcul des "dérivées partielles"  par rapport à : a , b , c , d par
* différence finie
debproc deriv a*flottant b*flottant c*flottant d*flottant fdep*listreel;
  deri_a= ((poly (a*1.01) b c  d) - fdep)/( a*0.01);
  deri_b= ((poly a (b*1.01) c  d) - fdep)/ (b*0.01);
  deri_c =((poly a b (c*1.01)  d) - fdep)/(C*0.01);
  deri_d =((poly a b c  (d*1.01)) - fdep)/(D*0.01);
finproc deri_a deri_b deri_c deri_d ;
 
deri_a deri_b deri_c deri_d = deriv a b c d fdep;
ss = moca depa fG fdep deri_a deri_b deri_c deri_d ;
 
list ss;
sol = prog 2. 3. 4. 5. ;
er = abs ( ss - sol);
ermax1= maxi er;
 
* remarque : l'opérateur MOCA cherche à minimiser l'écart en moyenne
* quadratique entre la base expérimentale et la fonction F(a,b,c,d),
* supposée linéaire en a b c d. Le polynome Y étant bien linéaire
* en a b c d, il trouve la "bonne" solution.
 
 
*
* envisageons maintenant que la fonction Y soit y= a + bx + cx**d
*
* on peut  utiliser la procedure ajuste
 
* parametres lineaires a b c
* parametre non lineaire d
* on peut ecrire la fonction sous la forme:
 
* y = a * f1 + b* f2 + c* f3;
* avec f1= 1    f2 = X     f3 = x**d
 
debp fct xtab*table p*listreel;
    tab1=table;
    tab2=table;
    x = xtab . 1;
    n=dime x;
*   fonction f1
    tab1.1 = prog n*1.;
    tab1.2 = x;
    tab1.3 = x**( extr p 1);
    tab2.'F'= tab1;
finproc tab2;
debproc deri xtab*table  p*listreel;
    tab1=table;
    tab2=table;
    tabg=table;
    tabf=table;
    tab=table;
    x = xtab . 1;
    n=dime x;
*   df1/dp1
    tab1 . 1=prog n*0.;
*   df2/dp1
    tab1. 2 = prog n*0.;
*   df3/dp1
    tab1. 3= ( log x) * (x** ( extr p 1));
    tabf . 1 = tab1;
    tab. 'F' = tabf;
finproc tab;
 
* fabriquons un cas "experimental" en
* posant a =1 b=-6 c=4 d=3 e et en prenant pour base des abscisses x
sol2= prog 1. -6.  4.  3.;
x = prog   0.1 0.3 0.5 1 2.5 3.6 6. 7.8 9 11; nexp=10;
* on va calculer nos "resultats experimentaux" : fobj
ct= prog nexp*1.;
bx = -6. * X;
cxpd= 4*( X ** 3);
fobj= ct + bx + cxpd;
 
* preparation des données pour appel ajuste
k=1;
L=3;
xtab=table;
xtab. 1= x;
tab1=table;
tab1.'X' = xtab;
tab1. 'F' = fobj;
tab1. 'K'= k;
tab1.'L'= L;
tab1.'PMIN'=prog 0.  -12   -20  -2   ;
tab1.'PMAX' = prog 10.  26.  12.  4.     ;
 
temps;
p  q = ajuste tab1 ;
temps;
mess ' sortie de la procedur ajuste valeur trouvée :';
mess ' a = ' ( extr p 1)  ' valeur attendue   1.';
mess ' b = ' ( extr p 2)  ' valeur attendue  -6.';
mess ' c = ' ( extr p 3)  ' valeur attendue   4.';
mess ' d = ' ( extr q 1)  ' valeur attendue   3.';
ltrou= p et q;
era= abs (ltrou - sol2);
ermax2= maxi era;
*
* on peut aussi tenter de resoudre le probleme à l'aide de MOCA
* Comme la fonction n'est pas linéaire en fonction du parametre d
* la solution fournie ne sera pas la bonne. on utilise le résultat
* trouvé pour nous donner une direction de descente le long de laquelle
* on cherche un minimun de la "distance" entre Fcalculée et Fobjectif
* une approche par itération est necessaire
*
 
para = prog  2. 2. 2. 2.;
 
* procedure pour calculer la fonction
debproc fcal para*listreel x*listreel;
  n = dime x;
  fa=  prog n*( extr para 1) ;
  fb= ( extr para 2)*x ;
  fc= ( extr para 3) * ( x ** ( extr para 4));
  ff= fa + fb + fc;
finproc ff;
* procedur pour calculer les derivees partielles
debproc fderiv para*listreel x*listreel ;
  n=dime x;
  fdera= prog n*1.;
  fderb= x * 1;
  fderc= x**(extr para 4);
  fderd = (extr para 3)*( log x) * (x** ( extr para 4));
finproc fdera fderb fderc fderd;
* procedure pour calculer le critère de distance
debproc criter fcalc*listreel fobjec*listreel ;
  cri= 0.;
  repe aa ( dime fcalc);
    cri = cri + (  ((extr fcalc &aa) - (extr fobjec &aa) )**2);
  fin aa;
finproc cri;
 
* schema iteratif
temps;
repeter iter 500;
 
  fdep = fcal para x; crii = criter fdep fobj;
  mess ' debut iteration ' &iter ' critère ' crii;
 
*
* on teste la convergence ( coutnouv - coutanc ) / coutanc < 1.e-4
*
  si ( &iter. ega 1) ;
    crianc = crii;
  sinon;
    cri_conv= ( crianc - crii ) / crianc;
    mess ' critère de convergence ' cri_conv;
    crianc= crii;
    si ( cri_conv < 1.e-4) ;
      mess  ' convergence à l itération ' &iter ' cout ' crii;
      quitter iter;
    finsi;
  finsi;
 
  fdera fderb fderc fderd = fderiv para x ;
*  list fobj;list fdep; list fdera ; list fderb; list fderc; list fderd;
  nouvpara= moca para fobj fdep fdera fderb fderc fderd;
*  list nouvpara;
*  
*   recherche le long de la direction donnée par nouvpara - para
*  
  desc= (nouvpara - para ) * 0.333333333;
 
  imu=1;
  repe ide 30;
    nouvpara= para + (desc * imu);
    fnou= fcal nouvpara x;
    cria= criter fnou fobj;
    si ( cria > crii );
      si ( imu ega 1) ;
        si (&ide < 6) ;
          desc = desc / 10.;
          iterer ide;
        sinon;
          mess ' pas de longueur trouvée le long de la descente';
          quitter iter;
        finsi;
      finsi;
*     recherche d'un min par approximation parabilique
      aa = cri2 + cria - (2.*crii) / 2.;
      bb= cria - cri2 / 2.;
      xx = bb / -2. / aa ;
      iies=   xx - 1. ;
      para= nouvpara + ( iies * desc);
      quitter ide;
    sinon;
      si ((imu ega &ide ) et  (imu ega 9)) ;
        para= nouvpara;quitter ide;
      finsi;
        imu=imu+1 ;
        cri2=crii;
        crii=cria;
    finsi;
  fin ide;
fin iter;
temps;
mess ' résutats par utilisation de moca itératif';
mess ' a = ' ( extr para 1)  ' valeur attendue   1.';
mess ' b = ' ( extr para 2)  ' valeur attendue  -6.';
mess ' c = ' ( extr para 3)  ' valeur attendue   4.';
mess ' d = ' ( extr para 4)  ' valeur attendue   3.';
erb= abs ( para - sol2);
ermax3= maxi erb;
*
*
* utilisation de l'opérateur levmar
*
* il faut définir la procédure feval qui evalue la fonction
* à minimiser et qui calcule les derivées partielles
*
debproc feval x*listreel   para*listreel;
  dy=prog;
  n=dime x;
  m = dime para;
  a1= extr para 1 ;
  b1= extr para 2;
  c1= extr para 3;
  d1= extr para 4;
  aa1 = prog n*a1;
 
  y= aa1 + ( b1 * x ) + ( c1 * ( x ** d1)) ;
 
  fdera= prog n*1.;
  fderb= x * 1;
  fderc= x**(extr para 4);
  fderd = (extr para 3)*( log x) * (x** ( extr para 4));
  l_dy=prog;
  ia=0;
  repe baa n;
    l_dy= l_dy et (prog ( extr fdera &baa));
    l_dy= l_dy et (prog ( extr fderb &baa));
    l_dy= l_dy et (prog ( extr fderc &baa));
    l_dy= l_dy et (prog ( extr fderd &baa));
  fin baa;
finp y l_dy;
 
aa = prog 2. 2. 2. 2.;sis = prog nexp*1. ;
temps;
a0 chi2 = 'LEVM' ABSC x ORDO fobj 'PARA' aa SIGM sis PROC feval  ;
temps;
mess ' résultats par utilisation de levm ' ;
mess ' a = ' ( extr a0 1)  ' valeur attendue   1.';
mess ' b = ' ( extr a0 2)  ' valeur attendue  -6.';
mess ' c = ' ( extr a0 3)  ' valeur attendue   4.';
mess ' d = ' ( extr a0 4)  ' valeur attendue   3.';
erc = abs ( a0 - sol2);
ermax4= maxi erc;
 
message ' erreur pour moca ' ermax1;
message ' erreur pour ajuste ' ermax2;
message ' erreur pour moca '  ermax3;
message ' erreur pour levm ' ermax4;
si ( (ermax1 + ermax2 + ermax3 + ermax4 ) > 1.e-4);
   erreur  5;
'SINON';
   'ERREUR' 0 ;   
finsi;
 
fin;