Télécharger la plaquette Cast3M


fr en

Télécharger hbm_vanderpol_force.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  hbm_vanderpol_force.dgibi
  2. ******************************************************************      
  3. *
  4. *  Systeme a un degre de liberte
  5. *
  6. *  Oscillateur Van der Pol force calcule en k harmoniques
  7. *
  8. *  m*x'' - c (1-x^2) *x' + k*x = f cos wt
  9. *  m=1; k=1; c=1; w=1; f = [ 0.1 ... 2 ]
  10. *  Continuation + HBM + AFT
  11. *
  12. ******************************************************************
  13.  
  14. *-----------------------------------------------------------------------
  15. * CHARMECA.PROCEDUR 
  16. * Calcul de la force non lineaire fnl(x,t) et de sa derivee dfnl(x,t)/dx
  17. *-----------------------------------------------------------------------
  18. *
  19. *  EN ENTREE : 
  20. *
  21. *  TAB1 
  22. *     parametres génériques :
  23. *  . 'DEP_NL'    : LISTCHPO DE DEPLACEMENTS TEMPORELLES
  24. *                  ou (pas encore programmé)
  25. *        . i . 'CHPOINT'  = CHPOINT unitaire des deplacements NL
  26. *        . i . 'LISTREEL' = evolution temporelle coefficient du chpoint
  27. *  . 'VIT_NL'    : VITESSES TEMPORELLES
  28. *     parametres propres a ce cas test :
  29. *  . 'COEFF_FNL' : coefficient de non linearite
  30. *  . 'POINT_FNL' : point de non linearite
  31. *  . 'COMP_FNL'  : composante de non linearite
  32. *
  33. *
  34. *  EN SORTIE :
  35. *
  36. *  TAB1
  37. *     parametres génériques :
  38. *   . 'FNL_T'  :   LISTCHPO DES FORCES NON LINEAIRES TEMPORELLES
  39. *                  ou TABLE d'indice : (pas encore programmé)
  40. *        . i . 'CHPOINT'  = CHPOINT unitaire des forces non lineaires
  41. *        . i . 'LISTREEL' = evolution temporelle coefficient du chpoint
  42. *   . 'KNL_T' :   DERIVEE DE FORCES NON LINEAIRES sous forme de TABLE :
  43. *        . i . j = 'LISTREEL' evolution temporelle de dFi/dxj 
  44. *
  45. *-----------------------------------------------------------------------
  46. DEBPROC CHARMECA TAB1*'TABLE';
  47.  
  48. ********     parametres d'entree propres a ce cas test   **********
  49.  
  50. * parametre de continuation
  51. time = TAB1 . 'TIME';
  52.  
  53. *coefficient de la force nl
  54. SI (exis TAB1 . 'COEFF_FNL');
  55.   AL = TAB1 . 'COEFF_FNL';
  56. SINON;
  57.   AL = 0.;
  58.   MESS 'CHARMECA: PAS D EFFET NON LINEAIRE !';
  59. FINSI;
  60. * LAMBDA = IPOL AL time;
  61. LAMBDA = AL ;
  62.  
  63. SI (EXIS TAB1 . 'POINT_FNL');
  64.   pNL = TAB1 . 'POINT_FNL';
  65. SINON;
  66.   MESS 'CHARMECA: IL MANQUE LA GEOMETRIE NL dans POINT_FNL'; ERRE 21;
  67. FINSI;
  68.  
  69. SI (EXIS TAB1 . 'COMP_FNL');
  70.   COMP_FNL = TAB1 . 'COMP_FNL';  
  71.   n_CNL = DIME COMP_FNL;
  72. SINON;
  73.   MESS 'CHARMECA: IL MANQUE LA COMPOSANTE NL dans COMP_FNL'; ERRE 21;
  74. FINSI;
  75.  
  76. ********     parametres d'entree génériques   **********
  77.  
  78. *DEPNL : list de champ par points, deplacements temporels du point nl
  79. *
  80. SI (EXIS TAB1 'DEP_NL');  
  81.   DEPNL = TAB1 . 'DEP_NL';
  82.   IFCHPO = ega (type DEPNL) 'LISTCHPO';
  83.   si (non IFCHPO);
  84.     si (neg  (type DEPNL) 'TABLE');
  85.       mess 'CHARMECA: DEP_NL doit etre de type LISTCHPO ou TABLE'; 
  86.       ERRE 21;
  87.     sinon; 
  88.       si (neg (type DEPNL . 1) 'LISTREEL');
  89.         mess 'CHARMECA: DEP_NL . i doit etre un LISTREEL'; 
  90.         ERRE 21;
  91.       finsi;
  92.     finsi;
  93.   finsi;
  94. SINON; 
  95.   MESS 'CHARMECA: IL MANQUE LA LISTE DES DEPLACEMENTS dans DEP_NL'; 
  96.   ERRE 21;
  97. FINSI;
  98.  
  99. * on a aussi besoin de la vitesse VITNL
  100. VITNL = TAB1 . 'VIT_NL';
  101.  
  102. * calcul de fnl ET dfnl/dx ? dans le doute on calcule tout...
  103. SI (NEG (TYPE IFFNL) 'LOGIQUE'); IFFNL = VRAI; FINSI;
  104. SI (NEG (TYPE IFKNL) 'LOGIQUE'); IFKNL = VRAI; FINSI; 
  105.  
  106.  
  107. ********  Calcul de FNL  ********
  108.  
  109. *---------------- Cas DEPNL LISTCHPO
  110. SI IFCHPO;
  111.  
  112.   erre 21;  
  113.  
  114. *---------------- Cas DEPNL TABLE de LISTREEL
  115. SINON; 
  116.  
  117. * DEPNL : table de listreel, indicée par les ddls NL
  118.  
  119. * nbt  : nombre de pas de temps
  120.   nbt = DIME DEPNL . 1; 
  121.   ib = 0;
  122.   FNLTEMP = TABLE; FNLTEMP . 1 = PROG;
  123.   KNLp = prog;
  124.   CNLp = prog;
  125.   REPE bfor nbt; ib = ib + 1; 
  126.     DEP_i = EXTR DEPNL . 1 ib;
  127.     VIT_i = EXTR VITNL . 1 ib;
  128.     FNL = LAMBDA * (1.-(DEP_i**2)) * VIT_i;
  129.     FNLTEMP . 1 = FNLTEMP . 1 et FNL; 
  130.     SI (IFKNL);
  131.         KNLp = KNLp et ( 2. * LAMBDA * DEP_i * VIT_i);
  132.         CNLp = CNLp et (-1. * LAMBDA * (1.-(DEP_i**2)));
  133.     FINSI;
  134.   FIN bfor;
  135.  
  136. FINSI;
  137.  
  138.  
  139. ******  stockage  ******
  140.  
  141. TAB1 . 'FNL_T'  = FNLTEMP;
  142. TAB1 . 'KNL_T' = tabl;
  143. TAB1 . 'CNL_T' = tabl;
  144.  
  145.  
  146. SI (IFKNL);
  147. * uniquement sous forme de LISTREELs
  148.    TAB1 . 'KNL_T' . 1 = tabl;
  149.    TAB1 . 'KNL_T' . 1 . 1 = KNLp;
  150.    TAB1 . 'CNL_T' . 1 = tabl;
  151.    TAB1 . 'CNL_T' . 1 . 1 = CNLp;
  152. FINSI;
  153.  
  154.  
  155. FINPROC;
  156. *-----------------------------------------------------------------------
  157.  
  158.  
  159.  
  160. ************************************************************************      
  161. *
  162. *   dynamique non lineaire
  163. *   reponse forcee d'un oscillateur de type Van der Pol
  164. *   reponse par HBM 
  165. *
  166. ************************************************************************
  167. *opti echo 0 ;
  168. opti debug 1;
  169. OPTI DIME 2 ELEM SEG2 MODE PLAN ;
  170.  
  171. GRAPH = FAUX ; SAVE = FAUX; TNR = VRAI;
  172. * GRAPH = VRAI ; SAVE = VRAI; TNR = FAUX;
  173.  
  174. *-------- Definition des parametres du calcul :
  175.  
  176. *nombre d'harmoniques utilises
  177. * nhbm = 1;
  178.  nhbm = 3;
  179. *  nhbm = 5;
  180. *  nhbm = 7;
  181. *  nhbm = 9;
  182.   nhbm = 11;
  183. * nhbm = 21;
  184.  
  185. *parametre de l'adimensionnement 
  186. u1max = 1.;
  187.  
  188. *des iterations
  189. * itmoy = 3.8;
  190. itmoy = 6;
  191. *ds = 0.1;
  192. ds = 1.;
  193. dsmax = 1.;
  194. dsmin = 1.E-3 ;
  195.  
  196. * si AFT,  OTFR: nb de points 2**OTFR pour TFR
  197. * OTFR = 7;
  198. OTFR = ENTI 'SUPERIEUR' (log (2 * (2 * nhbm + 1)) / (log 2)); 
  199. mess 'OTFR = ' OTFR;
  200.  
  201.  
  202. *-------- Definition des donnees d'oscillateur ---------------
  203.  
  204. * m*x'' - c (1-x^2) *x' + k*x = 0
  205. * m=1; k=1; c=[0..3]
  206.  
  207. *mass
  208. pmass = 1.;
  209. *rigidite
  210. prigi = 1.;
  211. * amortissement NL
  212. lambda = 1. ;
  213.  
  214. * prefix :
  215. prefix = chai 'hbm_vanderpol_force_n' nhbm ;
  216. si GRAPH ;
  217.   ficps = chai prefix '.ps';
  218.   opti TRAC PSC EPTR 6 POTR 'HELVETICA_16' 'FTRA' ficps;
  219. finsi;
  220.  
  221. *plage d'etude = plage du coefficient de la force non lineaire
  222. tcalc = prog 0.1 pas 0.1 1.9;
  223.  
  224. *------------- geometrie  ---------------    
  225. *
  226. P1 = 0. 0.  ; 
  227. P2 = 1. 0.  ;
  228. ST = P1 D 1  P2 ;                                                                                                    
  229.  
  230. *----------- construit des matrices caracteristiques ----------       
  231. *                                                                               
  232. LIM1 = MOTS UY ;
  233. MASS1 = MASS 'UY' pmass P2 ;
  234. RI1 = MANU 'RIGI' P2 LIM1 (prog prigi); 
  235.  
  236.  
  237.  
  238. *------- on définit une unique table pour tout (HBM et CONTINU) -------
  239. TAB1 = TABLE ;
  240.  
  241. *--------- remplissage des matrices "temporelles" ---------
  242.  
  243. * ATTENTION AUX UNITES : ON VA TRAVAILLER EN rad/s ==> pas de 2pi
  244. TAB1 . 'MASSE_CONSTANTE'        = MASS1 ;
  245. TAB1 . 'RIGIDITE_CONSTANTE'     = RI1 ; 
  246. TAB1 . 'N_HARMONIQUE'           = nhbm;
  247.  
  248. *------- remplissage des resultats attendus sur ddls "temporels" -------
  249. mycoul12 = mots VIOL BLEU BLEU TURQ TURQ OCEA OCEA VERT VERT OLIV OLIV 
  250. JAUN JAUN ORAN ORAN ROUG ROUG ROSE ROSE AZUR AZUR GRIS GRIS DEFA DEFA ;
  251. mycoul12= mycoul12 et (mots DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA 
  252. DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA 
  253. DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA);
  254. TAB1 . 'RESULTATS' = tabl; 
  255. ires1 = 1;
  256. TAB1 . 'RESULTATS' . ires1 = tabl; 
  257. TAB1 . 'RESULTATS' . ires1 . 'POINT_MESURE' = P2; 
  258. TAB1 . 'RESULTATS' . ires1 . 'COMPOSANTE'  = mot 'UY'; 
  259. TAB1 . 'RESULTATS' . ires1 . 'TITRE'       = mot 'U'; 
  260. TAB1 . 'RESULTATS' . ires1 . 'COULEUR'     = mot 'BLEU'; 
  261. TAB1 . 'RESULTATS' . ires1 . 'COULEUR_HBM' = mycoul12; 
  262.  
  263.  
  264. *--------- passage en frequentiel ---------
  265.  
  266. HBM TAB1;
  267.  
  268.  
  269. MESS '>>>>>>>>>>>>>>> COMPOSANTES DEP TEMPORELLES';
  270. LIST TAB1 .'COMPOSANTES' . 'DEPLACEMENT';
  271. MESS '>>>>>>>>>>>>>>> COMPOSANTES DEP FREQUENTIELLES';
  272. LIST TAB1 .'COMPOSANTES' . 'DEPLACEMENT_HBM'; 
  273.  
  274. * MESS 'RIGIDITE_HBM';            LIST TAB1 . 'RIGIDITE_HBM'       ; 
  275. * MESS 'BLOCAGES_HBM';            LIST TAB1 . 'BLOCAGES_HBM'       ; 
  276. * MESS 'AMORTISSEMENT_HBM';       LIST TAB1 . 'AMORTISSEMENT_HBM'  ; 
  277. * MESS 'MASSE_HBM';               LIST TAB1 . 'MASSE_HBM'          ; 
  278. * MESS 'MASSE_HBM_1';             LIST TAB1 . 'MASSE_HBM_1'        ; 
  279. * MESS 'MASSE_HBM_2';             LIST TAB1 . 'MASSE_HBM_2'        ; 
  280. * MESS 'AMORTISSEMENT_HBM_1';     LIST TAB1 . 'AMORTISSEMENT_HBM_1'; 
  281. list TAB1 . 'RESULTATS_HBM' ;
  282. TAB1 . 'AFFICHAGE_RESULTATS' = lect 2 3;
  283.  
  284.  
  285. *----------- definition du chargement ----------------------------------
  286.  
  287. *!ici, on definit directement le chargement en frequentiel
  288. *
  289. * FP11 = MANU 'CHPO' p2 1 'F2' 0. 'NATURE' 'DISCRET' ;
  290. FP11 = MANU 'CHPO' p2 1 'F2' 1. 'NATURE' 'DISCRET' ;
  291. *                                                         
  292. LIX1 = PROG 0. pas 0.1 5. ;  
  293. LIY1 = PROG (dime LIX1)*1.;   
  294. EV1 = EVOL MANU 't' LIX1 'F(t)' LIX1 ; 
  295. * dess EV1 POSX CENT POSY CENT;
  296. CHA1 = CHAR MECA FP11 EV1 ; 
  297.  
  298.  
  299. * frequence associee (en rad/s) = constante = 1
  300. Wev = evol manu 't' LIX1 'w' (LIY1);
  301. TAB1 . 'FREQUENCE' = Wev;
  302.  
  303.  
  304. *----------- OPTIONS DE CALCUL ----------------------
  305.  
  306. TAB1 . 'GRANDS_DEPLACEMENTS' = FAUX ;
  307. TAB1 . 'CHARGEMENT' = CHA1;
  308. TAB1 . 'PRECISION' = 1.E-5;
  309. TAB1 . 'TEMPS_CALCULES' = tcalc;
  310. TAB1 . 'MAXIPAS' = 300;
  311. TAB1 . 'ACCELERATION' = 4;
  312. TAB1 . 'NB_ITERATION' = itmoy;
  313. TAB1 . 'MAXITERATION' = 30;
  314. TAB1 . 'WTABLE' = tabl;
  315. TAB1 . 'WTABLE' . 'DS' = ds;
  316. TAB1 . 'WTABLE' . 'DSMAX' = dsmax;
  317. TAB1 . 'WTABLE' . 'DSMIN' = dsmin;
  318. TAB1 . 'MAXI_DEPLACEMENT' = u1max;
  319. TAB1 . 'PAS_SAUVES' = 1;
  320. TAB1 . 'NORME_RESIDU' = 1.;
  321.  
  322. * on donne la solution analytique pour lambda=0
  323. u0 = manu CHPO p2 1 'U2' 0. 'NATURE' 'DIFFUS';
  324. TAB1 . 'DEPLACEMENTS' = TABL;
  325. TAB1 . 'DEPLACEMENTS' . 0 = u0;
  326.  
  327. * calcul des forces NL
  328. TAB1 . 'PROCEDURE_CHARMECA'= VRAI ;
  329. TAB1 . 'CALC_CHPO' = FAUX;
  330. TAB1 . 'PROCEDURE_FREQUENCE_TEMPS' = mot 'AFT';
  331. TAB1 . 'N_PT_TFR' = OTFR;
  332. * 2 manieres de calculer la Jacobienne : 
  333. * directement ou par differentiation numerique
  334. * bp, 2018-10-01 : on passe a la differentiation numerique
  335. * car sinon, point initial singulier (nouvelle TFR trop precise!)
  336. TAB1 . 'JACOBIENNE' = mot 'DIFFERENTIATION';
  337.  
  338. * parametres du cas test lié a AFT et CHARMECA
  339. TAB1 . 'COEFF_FNL' = lambda;
  340. TAB1 . 'POINT_FNL' = P2;
  341. TAB1 . 'COMP_FNL'  = mots 'UY' ;
  342. TAB1 . 'CALC_VITE' = VRAI;
  343.  
  344. * recalcul de K apores la prediction
  345. * TAB1 . 'RECALCUL_K' = vrai;
  346.  
  347. * stabilité 
  348. TAB1 . 'STABILITE' = MOTS 'DIAG' 'HILL';
  349.  
  350.  
  351. *----------- resolution par la procedure CONTINU -----------------------  
  352.  
  353. *calcul avec continuation
  354. CONTINU TAB1;
  355.  
  356. TEMP 'IMPR' 'MAXI';
  357.  
  358. *------------------------- post-traitement ----------------------------- 
  359.  
  360. tit1 = mot 'Van der Pol - \l=1 -\w=1';
  361. TITRE tit1;
  362.  
  363. ********** tracé de toutes les harmoniques individuellement **********
  364.  
  365. Fprog = TAB1  . 'TEMPS_PROG';
  366. evtot = TAB1 . 'RESULTATS_HBM' . 'RESULTATS_EVOL';
  367.  
  368. TDESS1 = TABL;
  369. TDESS1 . 'TITRE' = TAB1 . 'RESULTATS_HBM' . 'TITRE';
  370.  
  371. si GRAPH;
  372.   dess evtot 
  373.      TITX 'F' POSX CENT  
  374.      TITY 'U_{k}' POSY CENT TDESS1 LEGE NO;
  375. finsi;
  376.  
  377. *******************   uy2p  : amplitude en norme 2  ******************
  378.  
  379. ihbm = 0;
  380. uy2p =  TAB1 . 'RESULTATS_HBM' . 1 . 'RESULTATS' ** 2;
  381. repe BUYHBM nhbm; ihbm = ihbm + 1;
  382.   ikcos = 2*ihbm;
  383.   iksin = 2*ihbm + 1;
  384.   uy2p = uy2p 
  385.   + (0.5* (  (TAB1 . 'RESULTATS_HBM' . ikcos . 'RESULTATS' ** 2)
  386.            + (TAB1 . 'RESULTATS_HBM' . iksin . 'RESULTATS' ** 2)) );
  387. fin  BUYHBM;
  388.  
  389. *** evolution + tracé *** 
  390. evuy2   = evol VIOL MANU Fprog uy2p;
  391. si GRAPH;
  392.   dess evuy2 
  393.      TITX 'F' POSX CENT  
  394.      TITY '|U|_{2}' POSY CENT ;
  395. finsi;
  396.  
  397.  
  398. ********** traduction des resultats frequentiels en temporels **********
  399. * on met 512 pas de temps pour le post traitement
  400. TAB1 . 'N_PT_POST' = 512 ;
  401. HBM_POST TAB1 (MOTS 'MAXI' 'TEMP' 'VITE');
  402.  
  403. ************* tracé du max d'amplitude max_t(u(t)) *********************
  404.  
  405. evuyamp = TAB1 . 'RESULTATS' . 'RESULTATS_EVOL';
  406.  
  407. si GRAPH;
  408.   dess evuyamp 
  409.     TITX 'F'    POSX CENT  
  410.     TITY 'max|U(t)|' POSY CENT ;
  411. finsi;
  412.  
  413. ************* tracé du max d'amplitude avec la stabilite ! *************
  414.   istab = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'STABILITE';
  415. Tdess2 = TABL;
  416. Tdess2 . 'LIGNE_VARIABLE' = TABL;
  417. Tdess2 . 'LIGNE_VARIABLE' . 1 = istab;
  418. si GRAPH;
  419.   dess evuyamp 
  420.     TITX 'F'    POSX CENT  
  421.     TITY 'max|U(t)|' POSY CENT Tdess2;
  422. finsi ;
  423.  
  424.  
  425. ************* tracé des plans de phase *************
  426. * pour F=0.5
  427. F05 = 0.5;
  428. TRES1T = TAB1 . 'RESULTATS' . 1 . 'RESULTATS_TEMPORELS';
  429.  
  430. evt05 = vide 'EVOLUTIO' ;
  431. evxy05 = vide 'EVOLUTIO' ;
  432. coul05 = mots 'BLEU' 'VERT' 'ROUG'; ic05 = 0;
  433. Upas05 = lect 5 29 55;
  434. REPE BPAS ((DIME TRES1T)/2 - 1);
  435.   ipas = &BPAS;
  436.   Upas1 = TAB1 . 'TEMPS' . ipas;
  437.   Upas2 = TAB1 . 'TEMPS' . (ipas + 1);
  438.   dU1 = Upas1 - F05;
  439.   dU2 = Upas2 - F05;
  440.   si (dU1 * dU2 > 0); iter BPAS; finsi;
  441. * comme u(t)=\sumt f_i(t) U_i, 
  442. * l'interpolation lineaire entre 2 pas U^1 et U^2 est equivalente 
  443. * a l'interpolation entre u^1(t) et u^2(t)
  444.   yp1 = extr TRES1T . ipas 'ORDO';
  445.   zp1 = extr TRES1T . (-1*ipas) 'ORDO';
  446.   yp2 = extr TRES1T . (ipas + 1) 'ORDO';
  447.   zp2 = extr TRES1T . (-1*(ipas + 1)) 'ORDO';
  448.   dU12 = Upas2 - Upas1;
  449.   yp = ((dU2 / dU12) * yp1) - ((dU1 / dU12) * yp2);
  450.   zp = ((dU2 / dU12) * zp1) - ((dU1 / dU12) * zp2);
  451.   evyz = EVOL 'VIOL' 'MANU' 'u' (yp)  'du/dt' (zp);
  452. *trop de tracés!   cha1pas = chai tit1 '- F=' FORMAT '(F5.3)' Upas;
  453. *trop de tracés!  dess (evyz) 'CARR' 'TITRE' cha1pas; 
  454. *trop de tracés! ==> on trace seulement pour F=0.5 :
  455. *   si (exis Upas05 ipas); 
  456.     ic05 = ic05 + 1;
  457.     evt05 = evt05  et (TRES1T . ipas coul (extr coul05 ic05));
  458.     evxy05= evxy05 et (evyz coul (extr coul05 ic05));
  459. *   finsi;
  460. FIN  BPAS;
  461.  
  462. si GRAPH;
  463.   dess (evt05) 'TITRE' (chai tit1 '- F=0.5');  
  464.   dess (evxy05) 'CARR' 'TITRE' (chai tit1 '- F=0.5');  
  465. finsi;   
  466.  
  467.  
  468.  
  469. si SAVE;
  470. *** recup des resultats d'integration temporelle de Matlab ***
  471. * ~/dynamique/oscillateurs/oscil_vanderpol/oscil4_2016_f05.m
  472.  
  473. opti 'ACQU' 'vdp05_1a';
  474. na = 3*2050;
  475. ACQU vdp_1a*'LISTREEL' na;
  476. * on ne recupere que les 512 derniers point = la derniere periode
  477. ideb_1a = (na - (3*512));
  478. t_1a = extr vdp_1a (lect (ideb_1a + 1) pas 3 (na - 2));
  479. u_1a = extr vdp_1a (lect (ideb_1a + 2) pas 3 (na - 1));
  480. v_1a = extr vdp_1a (lect (ideb_1a + 3) pas 3 (na - 0));
  481. t_1a = (t_1a - (extr t_1a 1)) / (2.*pi);
  482. opti 'ACQU' 'vdp05_1b';
  483. ACQU vdp_1b*'LISTREEL' na;
  484. t_1b = extr vdp_1b (lect (ideb_1a + 1) pas 3 (na - 2));
  485. u_1b = extr vdp_1b (lect (ideb_1a + 2) pas 3 (na - 1));
  486. v_1b = extr vdp_1b (lect (ideb_1a + 3) pas 3 (na - 0));
  487.  
  488. opti 'ACQU' 'vdp05_2a';
  489. na = 3*10000;
  490. ACQU vdp_2a*'LISTREEL' na;
  491. t_2a = extr vdp_2a (lect 1 pas 3 29998);
  492. u_2a = extr vdp_2a (lect 2 pas 3 29999);
  493. v_2a = extr vdp_2a (lect 3 pas 3 30000);
  494.  
  495. opti 'ACQU' 'vdp05_2b';
  496. ACQU vdp_2b*'LISTREEL' na;
  497. t_2b = extr vdp_2b (lect 1 pas 3 29998);
  498. u_2b = extr vdp_2b (lect 2 pas 3 29999);
  499. v_2b = extr vdp_2b (lect 3 pas 3 30000);
  500.  
  501. * petite reduction pour n'a'
  502. evt_1a  = evol 'ROSE' manu 't' t_1a 'u'     u_1a;
  503. evxy_1a = evol 'ROSE' manu 'u' u_1a 'du/dt' v_1a;
  504. evt_1b  = evol 'VIOL' manu 't' t_1b 'u'     u_1b;
  505. evxy_1b = evol 'VIOL' manu 'u' u_1b 'du/dt' v_1b;
  506.  
  507. evt_2a  = evol 'ROSE' manu 't' t_2a 'u'     u_2a;
  508. evxy_2a = evol 'ROSE' manu 'u' u_2a 'du/dt' v_2a;
  509. evt_2b  = evol 'VIOL' manu 't' t_2b 'u'     u_2b;
  510. evxy_2b = evol 'VIOL' manu 'u' u_2b 'du/dt' v_2b;
  511.  
  512. * dess (evt_1a et evt_1b);
  513. * dess (evt_2a et evt_2b);
  514. trevb = tabl; 
  515. trevb . 1 = mot 'TIRL';
  516. trevb . 2 = mot 'TIRR';
  517. dess (evxy_2a et evxy_2b) 'CARR' trevb;
  518. dess (evxy_1a et evxy_1b) 'CARR' trevb;
  519.  
  520. *** comparaison ***
  521. Tcomp =tabl;
  522. Tcomp . 4 = trevb . 1;
  523. Tcomp . 5 = trevb . 2;
  524. dess (evxy05 et evxy_2a et evxy_2b) 'CARR' Tcomp
  525. 'TITRE' (chai tit1 '- F=0.5');  
  526. dess (evxy05 et evxy_1a et evxy_1b) 'CARR' Tcomp
  527. 'TITRE' (chai tit1 '- F=0.5'); 
  528.  
  529. finsi;   
  530.  
  531.  
  532. * autres tracé de la stabilite : lambda_R en fonction de F
  533. mycoul = mots 'VIOL' 'BLEU' 'AZUR' 'TURQ' 'VERT' 'OLIV' 
  534.               'JAUN' 'ORAN' 'ROUG' 'ROSE';
  535. ncoul = dime mycoul;
  536. Fprog2 = enle Fprog (dime Fprog);
  537. evlrtot = VIDE 'EVOLUTIO' ;
  538. evlitot = VIDE 'EVOLUTIO' ;
  539. nl = dime TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL';
  540. il = 0;
  541. icoul = 0; coul1 = mot 'BLEU';
  542. repe bl nl; il = il + 1;  icoul = icoul + 1;
  543.   si(icoul > ncoul); icoul = 1; finsi; coul1 = extr mycoul icoul;
  544.   lr1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL' . il;
  545.   evlrtot = evlrtot  
  546.   et (evol coul1 MANU 'F' Fprog2 '\l_{R}' lr1);
  547.   li1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_IMAG' . il;
  548.   evlitot = evlitot  
  549.   et (evol coul1 MANU 'F' Fprog2 '\l_{I}' li1);
  550. fin  bl;
  551. Fprog3 = prog (mini Fprog2) (maxi Fprog2);
  552. evzero = evol manu 'F' Fprog3 '\l_{R}' (0.*Fprog3);
  553. Tdess3 = tabl;
  554. Tdess3 . 1 = mot 'TIRR';
  555. si GRAPH;
  556.   dess (evzero et evlrtot) 
  557.   'TITY' '\l_{R} (s^{-1})' Tdess3;
  558.   dess (evlitot) 'TITY' '\l_{I} (Hz)' ;
  559. finsi;
  560. * on cible les 2 premieres vp
  561. evlr12 = EXTR evlrtot 'COUR' (lect 1 2);
  562. evli12 = EXTR evlitot 'COUR' (lect 1 2);
  563. si GRAPH;
  564.   dess (evzero et evlr12) 
  565.   'TITY' '\l_{R} (s^{-1})' Tdess3;
  566.   dess (evzero et evli12) 
  567.   'TITY' '\l_{I} (s^{-1})' Tdess3;
  568. finsi;
  569.  
  570. * on met en evidence les solutions a F=0.5
  571. lr05_1 = IPOL F05 (extr evlr12 'COUR' 1) 'TOUS';
  572. lr05_2 = IPOL F05 (extr evlr12 'COUR' 2) 'TOUS';
  573. n_1 = dime lr05_1;
  574. n_2 = dime lr05_2;
  575. evr05_1 = EVOL 'VIOL' 'MANU' (prog n_1*F05) lr05_1;
  576. evr05_2 = EVOL 'BLEU' 'MANU' (prog n_2*F05) lr05_2;
  577. li05_1 = IPOL F05 (extr evli12 'COUR' 1) 'TOUS';
  578. li05_2 = IPOL F05 (extr evli12 'COUR' 2) 'TOUS';
  579. evi05_1 = EVOL 'VIOL' 'MANU' (prog n_1*F05) li05_1;
  580. evi05_2 = EVOL 'BLEU' 'MANU' (prog n_2*F05) li05_2;
  581. Tdess3 . 4 = mot 'NOLI MARQ PLUS';
  582. Tdess3 . 5 = mot 'NOLI MARQ LOSA';
  583. si GRAPH;
  584.   dess (evzero et evlr12 et evr05_1 et evr05_2) 
  585.   'TITY' '\l_{R} (s^{-1})' Tdess3;
  586.   dess (evzero et evli12 et evi05_1 et evi05_2) 
  587.   'TITY' '\l_{I} (s^{-1})' Tdess3;
  588. finsi;
  589.  
  590. * stabilite : trace dynamique pour faire une animation
  591. si GRAPH;
  592. TABORB = tabl;
  593. TABORB . 'PAS'       = 1;
  594. TABORB . 'QUEUE'     = MOT 'INFINIE';
  595. TABORB . 'EVOL_FIXE' = evzero ;
  596. TABORB . 'CARR'      = FAUX ;
  597.  
  598. ficps = chai prefix '-ORBITE-U.ps';
  599. opti 'TRAC' PSC 'FTRA' ficps;
  600. ORBITE evuyamp TABORB;
  601.  
  602. ficps = chai prefix '-ORBITE-lR.ps';
  603. opti 'TRAC' PSC 'FTRA' ficps;
  604. ORBITE evlr12 TABORB;
  605.  
  606. ficps = chai prefix '-ORBITE-lI.ps';
  607. opti 'TRAC' PSC 'FTRA' ficps;
  608. ORBITE evli12 TABORB;
  609. finsi;
  610.  
  611.  
  612.  
  613. *-------------------------- sauvegarde -----------------------------
  614.  
  615. si SAVE ;
  616.   ficxdr = chai prefix '.xdr'; 
  617.   mess 'sauvegarde dans ' ficxdr;
  618.   opti sauv ficxdr;
  620. finsi;
  621.  
  622. *--------------------- test de non regression ---------------------
  623.  
  624. si TNR ;
  625.  
  626. * reference = calcul Matlab integration temporelle directe ode23t
  627. * qq points du portrait de phase stable (= le 3eme)
  628.   tref_3 = prog 
  629.   0.0000      3.32031E-02  6.64062E-02  9.96094E-02  0.13281    
  630.   0.16602      0.19922      0.23242      0.26562      0.29883    
  631.   0.33203      0.36523      0.39844      0.43164      0.46484    
  632.   0.49805      0.53125      0.56445      0.59766      0.63086    
  633.   0.66406      0.69727      0.73047      0.76367      0.79687    
  634.   0.83008      0.86328      0.89648      0.92969      0.96289  0.99609;
  635.   uref_3 = prog 
  636.   -1.3333      -1.0770     -0.76083     -0.35992      0.15393
  637.   0.78173       1.4267       1.9024       2.1258       2.1685
  638.    2.1176       2.0190       1.8907       1.7376       1.5582
  639.    1.3469       1.0935      0.78151      0.38640     -0.12024
  640.  -0.74273      -1.3917      -1.8815      -2.1186      -2.1693
  641.   -2.1222      -2.0257      -1.8989      -1.7473      -1.5695 -1.3604;
  642.   si GRAPH;
  643.     evref_3 = evol JAUN MANU tref_3 uref_3;
  644.     tt = tabl; tt . 4 = mot 'NOLI MARQ PLUS';
  645.     dess (evt05 et evref_3) tt;
  646.   finsi;
  647. * Calcul de l'Aire des portraits de phase pour F=0.5 
  648. * (calculé avec la totalité des points)
  649.   Aref1 = 0.91267;
  650.   Aref3 = 17.520;
  651.  
  652. * Calcul de l'Aire des portraits de phase pour F=0.5
  653.   Axy = intg evxy05;
  654.   Axy1 = mini Axy;
  655.   Axy3 = maxi Axy;
  656.   mess Axy1 Axy3;
  657. * valeurs obtenues lors de la creatio ndu cas test : 0.91537  17.557 
  658. * et avec nhbm=3  : 0.91540          17.638
  659.   si ( ((ABS (Axy1 - Aref1  / Aref1)) > 0.01) 
  660.     ou ((ABS (Axy3 - Aref3  / Aref3)) > 0.01)); 
  661.     mess 'Erreur dans les portraits de phase !';
  662.     mess Axy1 Axy3;
  663.     ERRE 5; 
  664.   finsi;
  665.  
  666. * verifs ponctuelles de u(t) solution stable pour F=0.5
  667.   ev_3  = extr evt05 'COUR' 3;
  668.   n_3 =dime tref_3;
  669.   u_3 = IPOL tref_3 ev_3;  
  670.   err_3 = somm (abs (u_3 - uref_3)) / n_3;
  671.   si (err_3 > 0.05);
  672.     mess 'Erreur dans la solution stable pour F=0.5';
  673.     list ev_3;
  674.     mess err_3;
  675.     ERRE 5;     
  676.   finsi;
  677.  
  678. finsi;
  679.  
  680.  
  681. fin ;
  682.  
  683.  
  684.  
  685.  
  686.  
  687.  
  688.  
  689.  
  690.  
  691.  
  692.  
  693.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales