* fichier : testalea.dgibi
GRAPH = faux ;
********************* CAS TEST : testalea.dgibi ************************
*
* Cas-Test vérifiant le bon fonctionnement de l'opérateur 'ALEA'
* avec différentes options.
* Test de la bonne statistique du champ résultat.
*
************************************************************************
*
* Paramètres d'entrée :
* ---------------------
'OPTION' echo 1 ;
'OPTION' 'DIME' 2 'ELEM' 'QUA4' ;
'OPTION' 'TRACER' PSC ;
* Vérifie la syntaxe et les différentes possibilités d'usage ?
verisynt = vrai ;
* Vérifie la corection statistique des résultats ?
veriresu = vrai ;
* Longueur de corrélation :
portee = .7 ;
* Moyenne et Ecart-type du Logarithme népérien de perméabilité :
VALYM = -4.;
SIGMAY = 1. ;
VARIAN = sigmay * sigmay ;
*======================================================================
*== ==
*== PROCEDURES ==
*== ==
*======================================================================
'DEBPROC' FCT xtab*'TABLE' p*'LISTREEL' ;
* pour la procédure AJUSTE
* Calcul de la fonction objectif
* tbfonc . 'F' . i (fonction non-linéaire précédée d'un coef. inconnu n°i)
* tbfonc . 'G' (fonction non-linéaire non précédée d'un coef. inconnu)
* y = q^i f_i(x,p) + g(x,p)
x = xtab . 1 ;
n = 'DIME' x ;
b = 'EXTR' p 1;
tbfonc = 'TABLE';
tbfonc.'F' = 'TABLE';
tbfonc.'F' . 1 = 'EXP' (-1. / b * x);
'FINPROC' tbfonc;
'DEBPROC' DERI xtab*'TABLE' p*'LISTREEL';
* pour la procédure AJUSTE
* Calcul de la dérivée de la fonction objectif par rapport aux
* paramètres non linéaires
* tbderi . 'F' . i . k (dérivée par rapport au param k de la fonction f_i)
* tbderi . 'G' . k (dérivée par rapport au param k de la fonction g)
* dy_k / dp_k= q^i (df_i(x,p) / dp_k) + (dg(x,p) / dp_k)
x = xtab . 1;
n = 'DIME' x;
b = 'EXTR' p 1;
tbderi = 'TABLE';
tbderi . 'F' = 'TABLE';
tbderi . 'F' . 1 = 'TABLE';
tbderi . 'F' . 1 . 1 = x / b / b * ('EXP' (-1. * x / b)) ;
'FINPROC' tbderi;
*======================================================================
*== ==
*== PROGRAMME PRINCIPAL ==
*== ==
*======================================================================
* ===========================================
* Test de la syntaxe et des différents usages
* ===========================================
'SI' verisynt ;
* Création maillage (grossier) :
* ------------------------------
* Nombre de mailles pour 1 lambda :
Np = 1.;
* Dimensions :
lx = 2. * portee;
ly = 2. * portee;
lz = 2. * portee;
* Discrétisation :
* le +10d-10 sert à compenser l'imprécision machine pour ne pas avoir
* des valeurs de nx et ny trop petites de 1
nx = 'ENTIER' (lx / portee * Np + 0.0001) ;
ny = 'ENTIER' (ly / portee * Np + 0.0001) ;
nz = 'ENTIER' (lz / portee * Np + 0.0001) ;
* Points :
p1 = 0. 0. ;
p2 = lx 0. ;
p3 = lx ly ;
p4 = 0. ly ;
* Droites :
p1p2 = 'DROIT' p1 p2 nx ;
p2p3 = 'DROIT' p2 p3 ny ;
p3p4 = 'DROIT' p3 p4 nx ;
p4p1 = 'DROIT' p4 p1 ny ;
* Surface :
carre = 'SURFACE' (p1p2 'ET' p2p3 'ET' p3p4 'ET' p4p1) ;
* domaine
DOMHYB = 'CHANGER' CARRE 'QUAF';
* modèle :
MODHYB = 'MODELE' DOMHYB 'DARCY' 'ISOTROPE' ;
* points centres :
HYSOM = 'DOMA' MODHYB 'SOMMET';
HYCEN = 'DOMA' MODHYB 'CENTRE';
* -- Début tests syntaxiques --
* Syntaxe de base
YY = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
'LAMBDA1' portee (0. 1.) ;
YY = 'ALEA' 'BANDES_TOURNANTES' modhyb 'GRAVITE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
* Syntaxes alternatives
*- en 2D
* champ aléatoire dans une direction 1D
YY = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
'SI' GRAPH ;
'TITRE' 'Champ aleatoire unidirectionnel 2D, bandes tournantes';
'TRACER' MODHYB YY ;
'FINSI' ;
* champ aléatoire dans une direction 1D en spécifiant la direction
YY = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
'LAMBDA1' portee (0. 1.) ;
* champ point centre
YY = 'ALEA' 'BANDES_TOURNANTES' HYCEN
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
* champ point sommet (on fournit en test le maillage de 'QUA4')
YY = 'ALEA' 'BANDES_TOURNANTES' CARRE
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
* champ par élément sommet
YY = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'NOEUD'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
* champ par élément centre
YY = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
* sans le mot-clef 'MOYENNE'
YY = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
'EXPO' 'SIGMA' sigmay 'LAMBDA' portee ;
* avec un modèle mécanique
modmec = 'MODELE' CARRE 'MECANIQUE' ;
YY = 'ALEA' 'BANDES_TOURNANTES' modmec 'MASSE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
YY = 'ALEA' 'BANDES_TOURNANTES' modmec 'STRESSES'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
* avec un double modèle dont l'un seulement convient
mod2 = modhyb 'ET' modmec ;
YY = 'ALEA' 'BANDES_TOURNANTES' mod2 'RIGIDITE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
* anisotrope penché
portee1 = 3*portee ;
portee2 = .2*portee ;
vec1 = 2. 1. ;
vec2 = -1. 2.;
YY = 'ALEA' 'BANDES_TOURNANTES' CARRE
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
'LAMBDA1' portee1 vec1 'LAMBDA2' portee2 vec2 ;
*- en 3D
'SI' vrai ;
'OPTION' 'DIME' 3 'ELEM' 'CUB8' ;
cube = 'VOLUME' carre nz 'TRANSLATION' (0 0 lz) ;
domhyb3D = 'CHANGER' cube 'QUAF' ;
modhyb3D = 'MODELE' DOMHYB3D 'DARCY' 'ISOTROPE' ;
* plan
YY = 'ALEA' 'BANDES_TOURNANTES' MODHYB 'GRAVITE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee;
'SI' GRAPH;
'TITRE' 'Champ aleatoire plan en 3D, bandes tournantes';
'TRACER' MODHYB YY ;
'FINSI' ;
* isotrope
YY = 'ALEA' 'BANDES_TOURNANTES' MODHYB3D 'GRAVITE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee;
'SI' GRAPH;
'TITRE' 'Champ aleatoire isotrope 3D, bandes tournantes';
'TRACER' MODHYB YY ;
'FINSI' ;
* unidirectionnel
YY = 'ALEA' 'BANDES_TOURNANTES' MODHYB3D 'GRAVITE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
'LAMBDA1' portee (0 0 1) ;
'SI' GRAPH;
'TITRE' 'Champ aleatoire monodirectionnel 3D, bandes tournantes';
'TRACER' MODHYB YY ;
'FINSI' ;
* bidirectionnel
YY = 'ALEA' 'BANDES_TOURNANTES' MODHYB3D 'GRAVITE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
'LAMBDA1' portee 'LAMBDA2' portee ;
* 'SI' GRAPH;
'TITRE' 'Champ aleatoire bidirectionnel 3D, bandes tournantes';
'TRACER' MODHYB YY ;
* 'FINSI' ;
* tridirectionnel stratifié
YY = 'ALEA' 'BANDES_TOURNANTES' MODHYB3D 'GRAVITE'
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM
'LAMBDA1' portee 'LAMBDA2' portee 'LAMBDA3' (.1*portee);
'SI' GRAPH;
'TITRE' 'Champ aleatoire stratifie 3D, bandes tournantes';
'TRACER' MODHYB YY ;
'FINSI' ;
'OPTION' 'DIME' 2 ;
'FINSI';
'FINSI' ;
* =====================================
* Test de la qualité de la génération :
* =====================================
* en 2D
'SI' veriresu ;
*-- on utilise un maillage plus fin et plus vaste :
portee1 = portee ;
portee2 = portee ;
vec1 = (1. 0.) ;
vec2 = (0. 1.) ;
lx = 20. * portee1;
ly = 20. * portee2;
Np = 5.;
nx = 'ENTIER' (lx / portee * Np + 0.0001) ;
ny = 'ENTIER' (ly / portee * Np + 0.0001) ;
p1 = 0. 0. ;
p2 = lx 0. ;
p3 = lx ly ;
p4 = 0. ly ;
p1p2 = 'DROIT' p1 p2 nx ;
p2p3 = 'DROIT' p2 p3 ny ;
p3p4 = 'DROIT' p3 p4 nx ;
p4p1 = 'DROIT' p4 p1 ny ;
carre = 'SURFACE' (p1p2 'ET' p2p3 'ET' p3p4 'ET' p4p1) ;
DOMHYB = 'CHANGER' CARRE 'QUAF';
MODHYB = 'MODELE' DOMHYB 'DARCY' 'ISOTROPE' ;
* génération :
YY = 'ALEA' 'BANDES_TOURNANTES' CARRE
'EXPO' 'SIGMA' sigmay 'MOYENNE' VALYM 'LAMBDA' portee ;
* Tracé :
'SI' GRAPH;
'TITRE' ('CHAINE' 'Champ aleatoire fin 2D, bandes tournantes, '
((lx / portee) @ARR 0) 'x' ((ly / portee) @ARR 0) ', l/'
(Np @ARR 0) ) ;
'TRACER' CARRE YY ('CONTOUR' CARRE) ;
'FINSI' ;
*-- Calcul de la moyenne
moyc = ('MINIMUM' ('RESULT' yy)) / ('NBNO' ('EXTR' yy 'MAIL')) ;
mess 'moyenne : ' (moyc @ARR 3) ', objectif : ' (valym @ARR 3) ;
err1 = 'ABS' ((moyc - valym) / ('MAXIMUM' ('PROG' moyc valym) 'ABS'));
*-- Calcul de l'écart-type
sigc = ( ('MINIMUM' ('RESULT' ((yy - moyc) * (yy - moyc))))
/ ('NBEL' ('EXTR' yy 'MAIL')) ) ** .5 ;
mess 'écart-type : ' (sigc @ARR 3) ', objectif : ' (sigmay @ARR 3) ;
err2 = 'ABS' ((sigc - sigmay) / sigmay) ;
******************************************************
'SI' faux ;
* l'opérateur STAT n'est pas encore évolué.
*-- Calcul de la fonction de corrélation
* ====================================
* on opère une transformation de l'espace réel qui ramène au cas
* isotrope de longueur de corrélation 1.
'DEPLACER' carre 'AFFI' (1./portee1) (0. 0.) vec1 ;
'DEPLACER' carre 'AFFI' (1./portee2) (0. 0.) vec2 ;
portee0 = 1. ;
lidist = 'PROG' 0. 'PAS' (portee0 / Np) lx ;
YYM = 'MANU' 'CHPO' CARRE 1 'SCAL' moyc ;
'SI' (portee1 'EGA' portee2) ;
ev0 = 'STAT' 'COVA' YY YYM lidist 'REGU' ;
'SINON' ;
ev0 = 'STAT' 'COVA' YY YYM lidist ;
'FINSI' ;
evcov = 'COULEUR' ('EXTR' ev0 'COURBE' 1) 'ROUGE' ;
evnb = 'COULEUR' (('EXTR' ev0 'COURBE' 2)) 'VERT' ;
evnb = evnb * (varian / ('MAXIMUM' ('EXTR' evnb 'ORDO'))) ;
liabs = 'EXTR' evcov 'ABSC' ;
licov = 'EXTR' evcov 'ORDO' ;
* on impose la valeur de la variance calculée avant (rétroactif)
'REMPLACER' licov 1 (sigc * sigc) ;
liordth = varian * ('EXP' (-1. * liabs / portee0)) ;
evcovth = 'EVOL' 'DEFA' 'MANU' 'h' liabs 'Cyy' liordth ;
*-- Calcul de la longueur de corrélation
* ====================================
* On cale la fonction evcov sur
* y = varian * exp (- l * x)
* sur les 10 premières longueurs :
litest = 'MASQUE' liabs 'EGINFE' (10. * portee0) ;
ntest = 'ENTIER' ('RESULT' litest) ;
liind = 'LECT' 1 'PAS' 1 ntest ;
liabs2 = 'EXTR' liabs liind ;
licov2 = 'EXTR' licov liind ;
* on part d'un jeu de couple yi = f(xi)
* utilisation de la procedure AJUSTE :
* on cherche a identifier les parametres a et b de la fonction
* y = a * (exp (- x / b))
*
* on part d'un jeu de couple yj = f(xj)
* p = ensemble des paramètres k mis-en-jeu dans des expressions non linéaires
* q = ensemble des paramètres i mis-en-jeu dans des expressions linéaires
*
* Table des couples (x_j,y_j) à atteindre :
* NOTA : l'indice 'F' suivant aurait dû s'appeller plutôt 'Y'.
tab1 = 'TABLE' ;
tab1.'X' = 'TABLE' ;
tab1.'X' . 1 = liabs ;
tab1.'F' = licov;
* nombre de paramètres non linéaires
tab1.'K' = 1;
* nombre de paramètres linéaires (nb de jeux de valeurs à caler)
tab1.'L' = 1;
* espace de discrétisation pour les paramètres
pmin = 'PROG' 0.001;
pmax = 'PROG' 10.;
prc = 'PROG' 0.0001;
tab1.'MXTER' = 150;
tab1.'PMIN' = pmin ;
tab1.'PMAX' = pmax ;
tab1.'PRECISION' = prc ;
* Calcul de l'ajustement de la variance et de la longueur de
* corrélation :
q p = AJUSTE tab1;
varical = 'EXTR' q 1 ;
lambcal = 'EXTR' p 1 ;
ycal = varical * ('EXP' (-1. / lambcal * liabs));
evcal = 'EVOL' 'ROUG' 'MANU' 'X' liabs 'Y1' ycal;
portcal = 1. / lambcal ;
err3 = 'ABS' ((portcal - portee0) / portee0) ;
*-- tracé de la corrélation spatiale
'SI' GRAPH;
tabdess = 'TABLE' ;
tabdess.1 = 'MOT' ' ';
tabdess.2 = 'MOT' 'TIRR';
tabdess.3 = 'MOT' 'TIRC';
tabdess.4 = 'MOT' 'TIRL';
tabdess.titre = table ;
tabdess.titre.1 = 'CHAINE' 'Bantou';
tabdess.titre.2 = 'CHAINE' 'Calage' ;
tabdess.titre.3 = 'CHAINE' 'Theorique' ;
tabdess.titre.4 = 'CHAINE' 'Nb couples' ;
titdess = 'Correlation du champ aleatoire, bandes tournantes' ;
'DESSIN' (evcov 'ET' evcal 'ET' evcovth 'ET' evnb)
'MIMA' 'LEGE' 'TITRE' titdess tabdess 'AXES'
'XBOR' 0. ('MINIMUM' ('PROG' lx (10.*portee0))) ;
'FINSI' ;
'SINON' ;
err3 = 0;
'FINSI' ;
******************************************************
*-- Les valeurs calculées doivent être conformes aux valeurs théoriques
* à 5% pour la moyenne et l'écart-type, 20% pour la longueur de corrélation.
mess 'erreurs sur la moyenne, écart-type et covariance (%) = ';
mess (100 * err1) (100 * err2) (100 * err3) ;
'SI' ((err1 > 0.05) 'OU' (err2 > 0.05) 'OU' (err3 > 0.20));
'ERREUR' 5;
'SINON';
'ERREUR' 0;
'FINSI';
'FINSI' ; 'COMM' 'finsi veriresu' ;
'TEMPS' 'SGAC' 'IMPR' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer