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. OPTI EPSI LINEAIRE;
  171.  
  172. GRAPH = FAUX ; SAVE = FAUX; TNR = VRAI;
  173. * GRAPH = VRAI ; SAVE = VRAI; TNR = FAUX;
  174.  
  175. *-------- Definition des parametres du calcul :
  176.  
  177. *nombre d'harmoniques utilises
  178. * nhbm = 1;
  179. nhbm = 3;
  180. * nhbm = 5;
  181. * nhbm = 7;
  182. * nhbm = 9;
  183. nhbm = 11;
  184. * nhbm = 21;
  185.  
  186. *parametre de l'adimensionnement
  187. u1max = 1.;
  188.  
  189. *des iterations
  190. * itmoy = 3.8;
  191. itmoy = 6;
  192. *ds = 0.1;
  193. ds = 1.;
  194. dsmax = 1.;
  195. dsmin = 1.E-3 ;
  196.  
  197. * si AFT, OTFR: nb de points 2**OTFR pour TFR
  198. * OTFR = 7;
  199. OTFR = ENTI 'SUPERIEUR' (log (2 * (2 * nhbm + 1)) / (log 2));
  200. mess 'OTFR = ' OTFR;
  201.  
  202.  
  203. *-------- Definition des donnees d'oscillateur ---------------
  204.  
  205. * m*x'' - c (1-x^2) *x' + k*x = 0
  206. * m=1; k=1; c=[0..3]
  207.  
  208. *mass
  209. pmass = 1.;
  210. *rigidite
  211. prigi = 1.;
  212. * amortissement NL
  213. lambda = 1. ;
  214.  
  215. * prefix :
  216. prefix = chai 'hbm_vanderpol_force_n' nhbm ;
  217. si GRAPH ;
  218. ficps = chai prefix '.ps';
  219. opti TRAC PSC EPTR 6 POTR 'HELVETICA_16' 'FTRA' ficps;
  220. finsi;
  221.  
  222. *plage d'etude = plage du coefficient de la force non lineaire
  223. tcalc = prog 0.1 pas 0.1 1.9;
  224.  
  225. *------------- geometrie ---------------
  226. *
  227. P1 = 0. 0. ;
  228. P2 = 1. 0. ;
  229. ST = P1 D 1 P2 ;
  230.  
  231. *----------- construit des matrices caracteristiques ----------
  232. *
  233. LIM1 = MOTS UY ;
  234. MASS1 = MASS 'UY' pmass P2 ;
  235. RI1 = MANU 'RIGI' P2 LIM1 (prog prigi);
  236.  
  237.  
  238.  
  239. *------- on définit une unique table pour tout (HBM et CONTINU) -------
  240. TAB1 = TABLE ;
  241.  
  242. *--------- remplissage des matrices "temporelles" ---------
  243.  
  244. * ATTENTION AUX UNITES : ON VA TRAVAILLER EN rad/s ==> pas de 2pi
  245. TAB1 . 'MASSE_CONSTANTE' = MASS1 ;
  246. TAB1 . 'RIGIDITE_CONSTANTE' = RI1 ;
  247. TAB1 . 'N_HARMONIQUE' = nhbm;
  248.  
  249. *------- remplissage des resultats attendus sur ddls "temporels" -------
  250. mycoul12 = mots VIOL BLEU BLEU TURQ TURQ OCEA OCEA VERT VERT OLIV OLIV
  251. JAUN JAUN ORAN ORAN ROUG ROUG ROSE ROSE AZUR AZUR GRIS GRIS DEFA DEFA ;
  252. mycoul12= mycoul12 et (mots DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA
  253. DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA
  254. DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA DEFA);
  255. TAB1 . 'RESULTATS' = tabl;
  256. ires1 = 1;
  257. TAB1 . 'RESULTATS' . ires1 = tabl;
  258. TAB1 . 'RESULTATS' . ires1 . 'POINT_MESURE' = P2;
  259. TAB1 . 'RESULTATS' . ires1 . 'COMPOSANTE' = mot 'UY';
  260. TAB1 . 'RESULTATS' . ires1 . 'TITRE' = mot 'U';
  261. TAB1 . 'RESULTATS' . ires1 . 'COULEUR' = mot 'BLEU';
  262. TAB1 . 'RESULTATS' . ires1 . 'COULEUR_HBM' = mycoul12;
  263.  
  264.  
  265. *--------- passage en frequentiel ---------
  266.  
  267. HBM TAB1;
  268.  
  269.  
  270. MESS '>>>>>>>>>>>>>>> COMPOSANTES DEP TEMPORELLES';
  271. LIST TAB1 .'COMPOSANTES' . 'DEPLACEMENT';
  272. MESS '>>>>>>>>>>>>>>> COMPOSANTES DEP FREQUENTIELLES';
  273. LIST TAB1 .'COMPOSANTES' . 'DEPLACEMENT_HBM';
  274.  
  275. * MESS 'RIGIDITE_HBM'; LIST TAB1 . 'RIGIDITE_HBM' ;
  276. * MESS 'BLOCAGES_HBM'; LIST TAB1 . 'BLOCAGES_HBM' ;
  277. * MESS 'AMORTISSEMENT_HBM'; LIST TAB1 . 'AMORTISSEMENT_HBM' ;
  278. * MESS 'MASSE_HBM'; LIST TAB1 . 'MASSE_HBM' ;
  279. * MESS 'MASSE_HBM_1'; LIST TAB1 . 'MASSE_HBM_1' ;
  280. * MESS 'MASSE_HBM_2'; LIST TAB1 . 'MASSE_HBM_2' ;
  281. * MESS 'AMORTISSEMENT_HBM_1'; LIST TAB1 . 'AMORTISSEMENT_HBM_1';
  282. list TAB1 . 'RESULTATS_HBM' ;
  283. TAB1 . 'AFFICHAGE_RESULTATS' = lect 2 3;
  284.  
  285.  
  286. *----------- definition du chargement ----------------------------------
  287.  
  288. *!ici, on definit directement le chargement en frequentiel
  289. *
  290. * FP11 = MANU 'CHPO' p2 1 'F2' 0. 'NATURE' 'DISCRET' ;
  291. FP11 = MANU 'CHPO' p2 1 'F2' 1. 'NATURE' 'DISCRET' ;
  292. *
  293. LIX1 = PROG 0. pas 0.1 5. ;
  294. LIY1 = PROG (dime LIX1)*1.;
  295. EV1 = EVOL MANU 't' LIX1 'F(t)' LIX1 ;
  296. * dess EV1 POSX CENT POSY CENT;
  297. CHA1 = CHAR MECA FP11 EV1 ;
  298.  
  299.  
  300. * frequence associee (en rad/s) = constante = 1
  301. Wev = evol manu 't' LIX1 'w' (LIY1);
  302. TAB1 . 'FREQUENCE' = Wev;
  303.  
  304.  
  305. *----------- OPTIONS DE CALCUL ----------------------
  306.  
  307. TAB1 . 'GRANDS_DEPLACEMENTS' = FAUX ;
  308. TAB1 . 'CHARGEMENT' = CHA1;
  309. TAB1 . 'PRECISION' = 1.E-5;
  310. TAB1 . 'TEMPS_CALCULES' = tcalc;
  311. TAB1 . 'MAXIPAS' = 300;
  312. TAB1 . 'ACCELERATION' = 4;
  313. TAB1 . 'NB_ITERATION' = itmoy;
  314. TAB1 . 'MAXITERATION' = 30;
  315. TAB1 . 'WTABLE' = tabl;
  316. TAB1 . 'WTABLE' . 'DS' = ds;
  317. TAB1 . 'WTABLE' . 'DSMAX' = dsmax;
  318. TAB1 . 'WTABLE' . 'DSMIN' = dsmin;
  319. TAB1 . 'MAXI_DEPLACEMENT' = u1max;
  320. TAB1 . 'PAS_SAUVES' = 1;
  321. TAB1 . 'NORME_RESIDU' = 1.;
  322.  
  323. * on donne la solution analytique pour lambda=0
  324. u0 = manu CHPO p2 1 'U2' 0. 'NATURE' 'DIFFUS';
  325. TAB1 . 'DEPLACEMENTS' = TABL;
  326. TAB1 . 'DEPLACEMENTS' . 0 = u0;
  327.  
  328. * calcul des forces NL
  329. TAB1 . 'PROCEDURE_CHAR_MECA'= VRAI ;
  330. TAB1 . 'CALC_CHPO' = FAUX;
  331. TAB1 . 'PROCEDURE_FREQUENCE_TEMPS' = mot 'AFT';
  332. TAB1 . 'N_PT_TFR' = OTFR;
  333. * 2 manieres de calculer la Jacobienne :
  334. * directement ou par differentiation numerique
  335. * TAB1 . 'JACOBIENNE' = mot 'DIFFERENTIATION';
  336.  
  337. * parametres du cas test lié a AFT et CHARMECA
  338. TAB1 . 'COEFF_FNL' = lambda;
  339. TAB1 . 'POINT_FNL' = P2;
  340. TAB1 . 'COMP_FNL' = mots 'UY' ;
  341. TAB1 . 'CALC_VITE' = VRAI;
  342.  
  343. * recalcul de K apores la prediction
  344. * TAB1 . 'RECALCUL_K' = vrai;
  345.  
  346. * stabilité
  347. TAB1 . 'STABILITE' = MOTS 'DIAG' 'HILL';
  348.  
  349.  
  350. *----------- resolution par la procedure CONTINU -----------------------
  351.  
  352. *calcul avec continuation
  353. CONTINU TAB1;
  354.  
  355. TEMP 'IMPR' 'MAXI';
  356.  
  357. *------------------------- post-traitement -----------------------------
  358.  
  359. tit1 = mot 'Van der Pol - \l=1 -\w=1';
  360. TITRE tit1;
  361.  
  362. ********** tracé de toutes les harmoniques individuellement **********
  363.  
  364. Fprog = TAB1 . 'TEMPS_PROG';
  365. evtot = TAB1 . 'RESULTATS_HBM' . 'RESULTATS_EVOL';
  366.  
  367. TDESS1 = TABL;
  368. TDESS1 . 'TITRE' = TAB1 . 'RESULTATS_HBM' . 'TITRE';
  369.  
  370. si GRAPH;
  371. dess evtot
  372. TITX 'F' POSX CENT
  373. TITY 'U_{k}' POSY CENT TDESS1 LEGE NO;
  374. finsi;
  375.  
  376. ******************* uy2p : amplitude en norme 2 ******************
  377.  
  378. ihbm = 0;
  379. uy2p = TAB1 . 'RESULTATS_HBM' . 1 . 'RESULTATS' ** 2;
  380. repe BUYHBM nhbm; ihbm = ihbm + 1;
  381. ikcos = 2*ihbm;
  382. iksin = 2*ihbm + 1;
  383. uy2p = uy2p
  384. + (0.5* ( (TAB1 . 'RESULTATS_HBM' . ikcos . 'RESULTATS' ** 2)
  385. + (TAB1 . 'RESULTATS_HBM' . iksin . 'RESULTATS' ** 2)) );
  386. fin BUYHBM;
  387.  
  388. *** evolution + tracé ***
  389. evuy2 = evol VIOL MANU Fprog uy2p;
  390. si GRAPH;
  391. dess evuy2
  392. TITX 'F' POSX CENT
  393. TITY '|U|_{2}' POSY CENT ;
  394. finsi;
  395.  
  396.  
  397. ********** traduction des resultats frequentiels en temporels **********
  398. * on met 512 pas de temps pour le post traitement
  399. TAB1 . 'N_PT_POST' = 512 ;
  400. HBM_POST TAB1 (MOTS 'MAXI' 'TEMP' 'VITE');
  401.  
  402. ************* tracé du max d'amplitude max_t(u(t)) *********************
  403.  
  404. evuyamp = TAB1 . 'RESULTATS' . 'RESULTATS_EVOL';
  405.  
  406. si GRAPH;
  407. dess evuyamp
  408. TITX 'F' POSX CENT
  409. TITY 'max|U(t)|' POSY CENT ;
  410. finsi;
  411.  
  412. ************* tracé du max d'amplitude avec la stabilite ! *************
  413. istab = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'STABILITE';
  414. Tdess2 = TABL;
  415. Tdess2 . 'LIGNE_VARIABLE' = TABL;
  416. Tdess2 . 'LIGNE_VARIABLE' . 1 = istab;
  417. si GRAPH;
  418. dess evuyamp
  419. TITX 'F' POSX CENT
  420. TITY 'max|U(t)|' POSY CENT Tdess2;
  421. finsi ;
  422.  
  423.  
  424. ************* tracé des plans de phase *************
  425. * pour F=0.5
  426. F05 = 0.5;
  427. TRES1T = TAB1 . 'RESULTATS' . 1 . 'RESULTATS_TEMPORELS';
  428.  
  429. evt05 = vide 'EVOLUTIO' ;
  430. evxy05 = vide 'EVOLUTIO' ;
  431. coul05 = mots 'BLEU' 'VERT' 'ROUG'; ic05 = 0;
  432. Upas05 = lect 5 29 55;
  433. REPE BPAS ((DIME TRES1T)/2 - 1);
  434. ipas = &BPAS;
  435. Upas1 = TAB1 . 'TEMPS' . ipas;
  436. Upas2 = TAB1 . 'TEMPS' . (ipas + 1);
  437. dU1 = Upas1 - F05;
  438. dU2 = Upas2 - F05;
  439. si (dU1 * dU2 > 0); iter BPAS; finsi;
  440. * comme u(t)=\sumt f_i(t) U_i,
  441. * l'interpolation lineaire entre 2 pas U^1 et U^2 est equivalente
  442. * a l'interpolation entre u^1(t) et u^2(t)
  443. yp1 = extr TRES1T . ipas 'ORDO';
  444. zp1 = extr TRES1T . (-1*ipas) 'ORDO';
  445. yp2 = extr TRES1T . (ipas + 1) 'ORDO';
  446. zp2 = extr TRES1T . (-1*(ipas + 1)) 'ORDO';
  447. dU12 = Upas2 - Upas1;
  448. yp = ((dU2 / dU12) * yp1) - ((dU1 / dU12) * yp2);
  449. zp = ((dU2 / dU12) * zp1) - ((dU1 / dU12) * zp2);
  450. evyz = EVOL 'VIOL' 'MANU' 'u' (yp) 'du/dt' (zp);
  451. *trop de tracés! cha1pas = chai tit1 '- F=' FORMAT '(F5.3)' Upas;
  452. *trop de tracés! dess (evyz) 'CARR' 'TITRE' cha1pas;
  453. *trop de tracés! ==> on trace seulement pour F=0.5 :
  454. * si (exis Upas05 ipas);
  455. ic05 = ic05 + 1;
  456. evt05 = evt05 et (TRES1T . ipas coul (extr coul05 ic05));
  457. evxy05= evxy05 et (evyz coul (extr coul05 ic05));
  458. * finsi;
  459. FIN BPAS;
  460.  
  461. si GRAPH;
  462. dess (evt05) 'TITRE' (chai tit1 '- F=0.5');
  463. dess (evxy05) 'CARR' 'TITRE' (chai tit1 '- F=0.5');
  464. finsi;
  465.  
  466.  
  467.  
  468. si SAVE;
  469. *** recup des resultats d'integration temporelle de Matlab ***
  470. * ~/dynamique/oscillateurs/oscil_vanderpol/oscil4_2016_f05.m
  471.  
  472. opti 'ACQU' 'vdp05_1a';
  473. na = 3*2050;
  474. ACQU vdp_1a*'LISTREEL' na;
  475. * on ne recupere que les 512 derniers point = la derniere periode
  476. ideb_1a = (na - (3*512));
  477. t_1a = extr vdp_1a (lect (ideb_1a + 1) pas 3 (na - 2));
  478. u_1a = extr vdp_1a (lect (ideb_1a + 2) pas 3 (na - 1));
  479. v_1a = extr vdp_1a (lect (ideb_1a + 3) pas 3 (na - 0));
  480. t_1a = (t_1a - (extr t_1a 1)) / (2.*pi);
  481. opti 'ACQU' 'vdp05_1b';
  482. ACQU vdp_1b*'LISTREEL' na;
  483. t_1b = extr vdp_1b (lect (ideb_1a + 1) pas 3 (na - 2));
  484. u_1b = extr vdp_1b (lect (ideb_1a + 2) pas 3 (na - 1));
  485. v_1b = extr vdp_1b (lect (ideb_1a + 3) pas 3 (na - 0));
  486.  
  487. opti 'ACQU' 'vdp05_2a';
  488. na = 3*10000;
  489. ACQU vdp_2a*'LISTREEL' na;
  490. t_2a = extr vdp_2a (lect 1 pas 3 29998);
  491. u_2a = extr vdp_2a (lect 2 pas 3 29999);
  492. v_2a = extr vdp_2a (lect 3 pas 3 30000);
  493.  
  494. opti 'ACQU' 'vdp05_2b';
  495. ACQU vdp_2b*'LISTREEL' na;
  496. t_2b = extr vdp_2b (lect 1 pas 3 29998);
  497. u_2b = extr vdp_2b (lect 2 pas 3 29999);
  498. v_2b = extr vdp_2b (lect 3 pas 3 30000);
  499.  
  500. * petite reduction pour n'a'
  501. evt_1a = evol 'ROSE' manu 't' t_1a 'u' u_1a;
  502. evxy_1a = evol 'ROSE' manu 'u' u_1a 'du/dt' v_1a;
  503. evt_1b = evol 'VIOL' manu 't' t_1b 'u' u_1b;
  504. evxy_1b = evol 'VIOL' manu 'u' u_1b 'du/dt' v_1b;
  505.  
  506. evt_2a = evol 'ROSE' manu 't' t_2a 'u' u_2a;
  507. evxy_2a = evol 'ROSE' manu 'u' u_2a 'du/dt' v_2a;
  508. evt_2b = evol 'VIOL' manu 't' t_2b 'u' u_2b;
  509. evxy_2b = evol 'VIOL' manu 'u' u_2b 'du/dt' v_2b;
  510.  
  511. * dess (evt_1a et evt_1b);
  512. * dess (evt_2a et evt_2b);
  513. trevb = tabl;
  514. trevb . 1 = mot 'TIRL';
  515. trevb . 2 = mot 'TIRR';
  516. dess (evxy_2a et evxy_2b) 'CARR' trevb;
  517. dess (evxy_1a et evxy_1b) 'CARR' trevb;
  518.  
  519. *** comparaison ***
  520. Tcomp =tabl;
  521. Tcomp . 4 = trevb . 1;
  522. Tcomp . 5 = trevb . 2;
  523. dess (evxy05 et evxy_2a et evxy_2b) 'CARR' Tcomp
  524. 'TITRE' (chai tit1 '- F=0.5');
  525. dess (evxy05 et evxy_1a et evxy_1b) 'CARR' Tcomp
  526. 'TITRE' (chai tit1 '- F=0.5');
  527.  
  528. finsi;
  529.  
  530.  
  531. * autres tracé de la stabilite : lambda_R en fonction de F
  532. mycoul = mots 'VIOL' 'BLEU' 'AZUR' 'TURQ' 'VERT' 'OLIV'
  533. 'JAUN' 'ORAN' 'ROUG' 'ROSE';
  534. ncoul = dime mycoul;
  535. Fprog2 = enle Fprog (dime Fprog);
  536. evlrtot = VIDE 'EVOLUTIO' ;
  537. evlitot = VIDE 'EVOLUTIO' ;
  538. nl = dime TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL';
  539. il = 0;
  540. icoul = 0; coul1 = mot 'BLEU';
  541. repe bl nl; il = il + 1; icoul = icoul + 1;
  542. si(icoul > ncoul); icoul = 1; finsi; coul1 = extr mycoul icoul;
  543. lr1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL' . il;
  544. evlrtot = evlrtot
  545. et (evol coul1 MANU 'F' Fprog2 '\l_{R}' lr1);
  546. li1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_IMAG' . il;
  547. evlitot = evlitot
  548. et (evol coul1 MANU 'F' Fprog2 '\l_{I}' li1);
  549. fin bl;
  550. Fprog3 = prog (mini Fprog2) (maxi Fprog2);
  551. evzero = evol manu 'F' Fprog3 '\l_{R}' (0.*Fprog3);
  552. Tdess3 = tabl;
  553. Tdess3 . 1 = mot 'TIRR';
  554. si GRAPH;
  555. dess (evzero et evlrtot)
  556. 'TITY' '\l_{R} (s^{-1})' Tdess3;
  557. dess (evlitot) 'TITY' '\l_{I} (Hz)' ;
  558. finsi;
  559. * on cible les 2 premieres vp
  560. evlr12 = EXTR evlrtot 'COUR' (lect 1 2);
  561. evli12 = EXTR evlitot 'COUR' (lect 1 2);
  562. si GRAPH;
  563. dess (evzero et evlr12)
  564. 'TITY' '\l_{R} (s^{-1})' Tdess3;
  565. dess (evzero et evli12)
  566. 'TITY' '\l_{I} (s^{-1})' Tdess3;
  567. finsi;
  568.  
  569. * on met en evidence les solutions a F=0.5
  570. lr05_1 = IPOL F05 (extr evlr12 'COUR' 1) 'TOUS';
  571. lr05_2 = IPOL F05 (extr evlr12 'COUR' 2) 'TOUS';
  572. n_1 = dime lr05_1;
  573. n_2 = dime lr05_2;
  574. evr05_1 = EVOL 'VIOL' 'MANU' (prog n_1*F05) lr05_1;
  575. evr05_2 = EVOL 'BLEU' 'MANU' (prog n_2*F05) lr05_2;
  576. li05_1 = IPOL F05 (extr evli12 'COUR' 1) 'TOUS';
  577. li05_2 = IPOL F05 (extr evli12 'COUR' 2) 'TOUS';
  578. evi05_1 = EVOL 'VIOL' 'MANU' (prog n_1*F05) li05_1;
  579. evi05_2 = EVOL 'BLEU' 'MANU' (prog n_2*F05) li05_2;
  580. Tdess3 . 4 = mot 'NOLI MARQ PLUS';
  581. Tdess3 . 5 = mot 'NOLI MARQ LOSA';
  582. si GRAPH;
  583. dess (evzero et evlr12 et evr05_1 et evr05_2)
  584. 'TITY' '\l_{R} (s^{-1})' Tdess3;
  585. dess (evzero et evli12 et evi05_1 et evi05_2)
  586. 'TITY' '\l_{I} (s^{-1})' Tdess3;
  587. finsi;
  588.  
  589. * stabilite : trace dynamique pour faire une animation
  590. si GRAPH;
  591. TABORB = tabl;
  592. TABORB . 'PAS' = 1;
  593. TABORB . 'QUEUE' = MOT 'INFINIE';
  594. TABORB . 'EVOL_FIXE' = evzero ;
  595. TABORB . 'CARR' = FAUX ;
  596.  
  597. ficps = chai prefix '-ORBITE-U.ps';
  598. opti 'TRAC' PSC 'FTRA' ficps;
  599. ORBITE evuyamp TABORB;
  600.  
  601. ficps = chai prefix '-ORBITE-lR.ps';
  602. opti 'TRAC' PSC 'FTRA' ficps;
  603. ORBITE evlr12 TABORB;
  604.  
  605. ficps = chai prefix '-ORBITE-lI.ps';
  606. opti 'TRAC' PSC 'FTRA' ficps;
  607. ORBITE evli12 TABORB;
  608. finsi;
  609.  
  610.  
  611.  
  612. *-------------------------- sauvegarde -----------------------------
  613.  
  614. si SAVE ;
  615. ficxdr = chai prefix '.xdr';
  616. mess 'sauvegarde dans ' ficxdr;
  617. opti sauv ficxdr;
  618. finsi;
  619.  
  620. *--------------------- test de non regression ---------------------
  621.  
  622. si TNR ;
  623.  
  624. * reference = calcul Matlab integration temporelle directe ode23t
  625. * qq points du portrait de phase stable (= le 3eme)
  626. tref_3 = prog
  627. 0.0000 3.32031E-02 6.64062E-02 9.96094E-02 0.13281
  628. 0.16602 0.19922 0.23242 0.26562 0.29883
  629. 0.33203 0.36523 0.39844 0.43164 0.46484
  630. 0.49805 0.53125 0.56445 0.59766 0.63086
  631. 0.66406 0.69727 0.73047 0.76367 0.79687
  632. 0.83008 0.86328 0.89648 0.92969 0.96289 0.99609;
  633. uref_3 = prog
  634. -1.3333 -1.0770 -0.76083 -0.35992 0.15393
  635. 0.78173 1.4267 1.9024 2.1258 2.1685
  636. 2.1176 2.0190 1.8907 1.7376 1.5582
  637. 1.3469 1.0935 0.78151 0.38640 -0.12024
  638. -0.74273 -1.3917 -1.8815 -2.1186 -2.1693
  639. -2.1222 -2.0257 -1.8989 -1.7473 -1.5695 -1.3604;
  640. si GRAPH;
  641. evref_3 = evol JAUN MANU tref_3 uref_3;
  642. tt = tabl; tt . 4 = mot 'NOLI MARQ PLUS';
  643. dess (evt05 et evref_3) tt;
  644. finsi;
  645. * Calcul de l'Aire des portraits de phase pour F=0.5
  646. * (calculé avec la totalité des points)
  647. Aref1 = 0.91267;
  648. Aref3 = 17.520;
  649.  
  650. * Calcul de l'Aire des portraits de phase pour F=0.5
  651. Axy = intg evxy05;
  652. Axy1 = mini Axy;
  653. Axy3 = maxi Axy;
  654. mess Axy1 Axy3;
  655. * valeurs obtenues lors de la creatio ndu cas test : 0.91537 17.557
  656. * et avec nhbm=3 : 0.91540 17.638
  657. si ( ((ABS (Axy1 - Aref1 / Aref1)) > 0.01)
  658. ou ((ABS (Axy3 - Aref3 / Aref3)) > 0.01));
  659. mess 'Erreur dans les portraits de phase !';
  660. mess Axy1 Axy3;
  661. ERRE 5;
  662. finsi;
  663.  
  664. * verifs ponctuelles de u(t) solution stable pour F=0.5
  665. ev_3 = extr evt05 'COUR' 3;
  666. n_3 =dime tref_3;
  667. u_3 = IPOL tref_3 ev_3;
  668. err_3 = somm (abs (u_3 - uref_3)) / n_3;
  669. si (err_3 > 0.05);
  670. mess 'Erreur dans la solution stable pour F=0.5';
  671. list ev_3;
  672. mess err_3;
  673. ERRE 5;
  674. finsi;
  675.  
  676. finsi;
  677.  
  678.  
  679. fin ;
  680.  
  681.  
  682.  
  683.  
  684.  
  685.  

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