* 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;
© Cast3M 2003 - All rights reserved.
Disclaimer