Télécharger corman.procedur

Retour à la liste

Numérotation des lignes :

  1. * CORMAN PROCEDUR CHAT 13/02/05 21:15:00 7691
  2.  
  3.  
  4.  
  5. 'DEBPROC' corman rit*'RIGIDITE' MO*'MMODEL' MA*'MCHAML'
  6. ORDRE*'ENTIER' U0*'CHPOINT' S0*'MCHAML' ffi*CHPOINT
  7. ETAB12*'TABLE';
  8. * quan,d on arrive ici on se trouve sur la config en debut
  9. * de pas. ffi est le résidu et s0 est la contrainte en début de pas
  10. *
  11. *
  12. * attention dans cette version on n'utilise pas les approximants de Padé.
  13. *
  14. sgdlamda='SIGN' 1 ;
  15. *** mess ' entree dans corman ' ;
  16. IDYN=ETAB12.'DYNAMIQUE';
  17. 'SI' IDYN;
  18. MMA= ETAB12. 'MASSE';
  19. QUADTT2=4.D0 / ETAB12.'DT' / ETAB12.'DT';
  20. MMAPLU=MMA * QUADTT2;
  21. * mess ' quadtt2 ' quadtt2;
  22. 'FINSI';
  23. N=ORDRE ;
  24. DEP1=U0;
  25. ncontmax=250;
  26. isol1=faux;
  27. prec = 1.D-8;
  28. ttdep=table; ttdep.0 = U0; ttlam=table ; ttlam.0=0.;
  29. ttsig=table; ttsig.0=s0;NA=0 ;
  30. depreel0 = U0 'ENLE' 'LX' ;
  31. depini=U0 'ENLE' 'LX' ;
  32. sigreel0 = S0 ;
  33. ffs0= 'BSIGMA' S0 MO MA;
  34. lamreel0 = 0.d0;
  35. IOUT=1;
  36. ritblo='EXTR' rit 'RIGI' 'MULT';
  37. * on veut lambda=1.d0
  38. llld = mots 'UX' 'UY' 'UZ' 'RX' 'RY' 'RZ' ;
  39. IMAS=FAUX;
  40. IDIMV= 'VALE' 'DIME';
  41. IMODD= 'VALE' 'MODE';
  42. 'SI' ( 'EXISTE' MO 'ELEM' 'TRI3' 'QUA4' 'TRI6' 'QUA8'
  43. 'CUB8' 'TET4' 'PRI6' 'PYR5' 'CU20' 'TE10' 'PR15' 'PY13' ) ;
  44. IMAS=VRAI;MOMAS= 'EXTR' MO 'ELEM''TRI3' 'QUA4' 'TRI6' 'QUA8'
  45. 'CUB8' 'TET4' 'PRI6' 'PYR5' 'CU20' 'TE10' 'PR15' 'PY13';
  46. 'FINSI';
  47. ICOQ=FAUX;
  48. 'SI' ( 'EXISTE' MO 'ELEM' 'COQ2' 'DKT' );ICOQ=VRAI;
  49. MOCOQ= 'EXTR' MO 'ELEM' 'COQ2' 'DKT'; 'FINSI';
  50. 'SI' ('EXIS' MO ELEM 'POUT'); IOUT=0;'QUITTER' CORMAN; 'FINSI';
  51. ff = ffi ;
  52. lam0 = 0. ;
  53. geoini = form;
  54. sigcour=s0;
  55. *m0
  56. deltatot= u0*0.;
  57. 'REPE' IT ;
  58. 'SI' ('EXIS' MO ELEM 'POUT'); IOUT=0;'QUITTER' IT; 'FINSI';
  59. 'SI' ('EXIS' MO ELEM 'COQ2'); IOUT=0;'QUITTER' IT; 'FINSI';
  60. *'SI' (&it > 1) ; mess ' continuation ' ( &it - 1) lam0;'FINSI';
  61. 'SI' (&it > ncontmax);MESS ' Pas de solution apres ' ncontmax
  62. 'continuations'; IOUT=0;
  63. 'QUITTER' IT; 'FINSI';
  64. **************************Initialisation**************************
  65.  
  66. *
  67. * quelques initialisation de table servant aux stockages intermediaire
  68. *
  69. tadep=table; talam=table;tasig=table;taksi = table;
  70. tagrad = table; tagraf=table;
  71. tabase = table;tabalpha = table; tabvald = table; tabpoly = table;
  72. *
  73. ************************** calcul du premier ordre******************
  74. *
  75. deltau0= u0*0.;
  76. kspp = 'KSIGMA' sigreel0 mo ma;
  77. rit= 'RIGI' MO MA;
  78. rit=rit 'ET' ritblo;
  79. ritot = rit 'ET' kspp;
  80. 'SI' IDYN;
  81. Ritot= ritot 'ET' MMAPLU;
  82. 'FINSI';
  83. DE1= 'RESOU' ritot ff 'NOUNIL';
  84. psc = 'XTX' de1 ; psc2 = psc + 1.D0 ** 0.5D0;
  85. lambda1 = 1.d0 / psc2 ;
  86. dep1 = de1 * lambda1;
  87. s1 = 'SIGMA' 'LINE' dep1 mo ma;
  88. *
  89. * calcul d'une norme de convergence
  90. *
  91. * remplissage des tables pour ordre 1
  92. tadep . 1 = dep1;
  93. talam . 1 = lambda1 ;
  94. tasig . 1 = s1;
  95. taksi . 1 = 'KSIGMA' s1 mo ma;
  96. tagrad . 1 = 'GRAD' dep1 mo ma;
  97. 'SI' ICOQ ;tagraf . 1 = 'GRAF' dep1 mo ma ;'FINSI';
  98. tabalpha. 1 = table;
  99. tabalpha . 1 . 1 = ('XTX' tadep. 1)**0.5D0 ;
  100. tabase. 1 = (tadep. 1)/(tabalpha. 1 . 1);
  101. p=1;
  102. *mess ' ordre ' p ' lambda ' lambda1;
  103. *temps place;
  104. *
  105. ******************** calcul des ordres superieurs*****************
  106. *
  107. *mess ' debut de boucle';
  108. 'REPETER' bordre (N - 1);
  109. p = &bordre + 1;
  110. *
  111. *calcul des ksig* (si) * u(p - i + 1);
  112. nfois = p - 1 ;
  113. 'REPETER' ks nfois;
  114. fksi = taksi . &ks * tadep .( p - &ks);
  115. 'SI' ( &ks 'EGA' 1);
  116. fks = -1.D0 * fksi;
  117. 'SINON';
  118. *menagement de l'espace memoire
  119. XXX = fks;
  120. fks = fks - fksi;
  121. 'FINSI';
  122. 'FIN' ks;
  123. fks = fks ;
  124. *
  125. * calcul de bsigma ( D *0.5( gradui * graduj))
  126. *
  127. 'REPETER' bs nfois;
  128. gradi = tagrad . &bs;
  129. gradj = tagrad . ( p - &bs );
  130. 'SI' ICOQ ;
  131. grafi = tagraf . &bs;
  132. grafj = tagraf . ( p - &bs );
  133. 'FINSI';
  134.  
  135. * les cartes suivantes sont pour les elements coques
  136. 'SI' ICOQ ;
  137. gradi = 'REDU' gradi MOCOQ; gradj = 'REDU' gradj MOCOQ;
  138. grafi = 'REDU' grafi MOCOQ; grafj = 'REDU' grafj MOCOQ;
  139.  
  140. Uxxi='EXCO' gradi 'UX,X' 'SCAL' chan TYPE 'SCAL';
  141. Uxyi='EXCO' gradi 'UX,Y' 'SCAL' chan TYPE 'SCAL' ;
  142. Uxzi='EXCO' gradi 'UX,Z' 'SCAL' chan TYPE 'SCAL' ;
  143. Uyxi='EXCO' gradi 'UY,X' 'SCAL' chan TYPE 'SCAL' ;
  144. Uyyi='EXCO' gradi 'UY,Y' 'SCAL' chan TYPE 'SCAL' ;
  145. Uyzi='EXCO' gradi 'UY,Z' 'SCAL' chan TYPE 'SCAL' ;
  146. Uzxi='EXCO' gradi 'UZ,X' 'SCAL' chan TYPE 'SCAL' ;
  147. Uzyi='EXCO' gradi 'UZ,Y' 'SCAL' chan TYPE 'SCAL' ;
  148.  
  149. Uxxj='EXCO' gradj 'UX,X' 'SCAL' chan TYPE 'SCAL' ;
  150. Uxyj='EXCO' gradj 'UX,Y' 'SCAL' chan TYPE 'SCAL' ;
  151. Uxzj='EXCO' gradj 'UX,Z' 'SCAL' chan TYPE 'SCAL' ;
  152. Uyxj='EXCO' gradj 'UY,X' 'SCAL' chan TYPE 'SCAL' ;
  153. Uyyj='EXCO' gradj 'UY,Y' 'SCAL' chan TYPE 'SCAL' ;
  154. Uyzj='EXCO' gradj 'UY,Z' 'SCAL' chan TYPE 'SCAL' ;
  155. Uzxj='EXCO' gradj 'UZ,X' 'SCAL' chan TYPE 'SCAL' ;
  156. Uzyj='EXCO' gradj 'UZ,Y' 'SCAL' chan TYPE 'SCAL' ;
  157. Bxxi='EXCO' grafi 'BX,X' 'SCAL' chan TYPE 'SCAL' ;
  158. Bxyi='EXCO' grafi 'BX,Y' 'SCAL' chan TYPE 'SCAL' ;
  159. Byxi='EXCO' grafi 'BY,X' 'SCAL' chan TYPE 'SCAL' ;
  160. Byyi='EXCO' grafi 'BY,Y' 'SCAL' chan TYPE 'SCAL' ;
  161. Bxxj='EXCO' grafj 'BX,X' 'SCAL' chan TYPE 'SCAL' ;
  162. Bxyj='EXCO' grafj 'BX,Y' 'SCAL' chan TYPE 'SCAL' ;
  163. Byxj='EXCO' grafj 'BY,X' 'SCAL' chan TYPE 'SCAL';
  164. * Nouvelle version
  165. epss=( 0.5 * ( Uxxi*Uxxj + (Uyxi*Uyxj) + (Uzxi*Uzxj ) ) )
  166. 'NOMC' 'EPSS';
  167. eptt=( 0.5 * ( Uxyi*Uxyj + (Uyyi*Uyyj) + (Uzyi*Uzyj ) ) )
  168. 'NOMC' 'EPTT';
  169. gast=(0.5 * (Uxyi*Uxxj + (Uyyi*Uyxj) + (Uzyi*Uzxj) + (Uxyj*Uxxi) +
  170. (Uyyj*Uyxi) + (Uzxi*Uzyj) ) )'NOMC' 'GAST';
  171. rtss=(epss *0.) 'NOMC' 'RTSS';
  172. rttt=rtss 'NOMC' 'RTTT';
  173. rtst=rtss 'NOMC' 'RTST';
  174.  
  175. ep2 = epss + eptt + gast + rtss + rttt + rtst;
  176.  
  177. 'FINSI';
  178.  
  179. * les cartes suivantes sont pour les elements massifs
  180.  
  181. 'SI' IMAS ;
  182. gradi = 'REDU' gradi MOMAS;gradj = 'REDU' gradj MOMAS;
  183. gru2 =(( 'TAGR' gradi) * gradj mo) * 0.5D0 ;
  184. 'SI' ('NEG' imodd 'AXIS');
  185. epxyz= 'EXCO' gru2 ('MOTS' 'UX,X' 'UY,Y' 'UZ,Z')
  186. ('MOTS' 'EPXX' 'EPYY' 'EPZZ');
  187. gaxyz= ('EXCO' gru2 ('MOTS' 'UX,Y' 'UX,Z' 'UY,Z')
  188. ('MOTS' 'GAXY' 'GAXZ' 'GAYZ'))*2.d0;
  189. ep2 = epxyz + gaxyz;
  190. 'SINON';
  191. epxx = 'EXCO' gru2 UR,R 'EPRR';
  192. epyy = 'EXCO' gru2 UZ,Z 'EPZZ';
  193. epzz = 'EXCO' gru2 UT,T 'EPTT';
  194. gaxy = ('EXCO' gru2 UR,Z 'GARZ')*2.D0 ;
  195. gaxz = ('EXCO' gru2 UR,T 'GART')*2.D0 ;
  196. gayz = ('EXCO' gru2 UT,Z 'GATZ')*2.D0 ;
  197. ep2 = epxx + epyy + epzz + gaxy + gaxz + gayz;
  198. 'FINSI';
  199. 'FINSI';
  200.  
  201. * fin du traitement spécifique
  202. ep2 = 'CHANG' ep2 'TYPE' 'DEFORMATIONS';
  203. 'SI' ( 'EGA' &bs 1);
  204. ep2tot = ep2 ;
  205. 'SINON';
  206. XXX = ep2;
  207. ep2tot = ep2tot + ep2;
  208. 'FINSI';
  209. 'FIN' bs;
  210. sig = 'ELAS' ep2tot mo ma ;
  211. fbs = ('BSIGMA' sig mo ma) * -1.D0 ;
  212. u2 = 'RESOU' ritot ( fks + fbs) 'NOUNIL';
  213. lamdai = (u2 'XTY' tadep.1 llld llld )* -1.D0 * lambda1;
  214.  
  215. *
  216. *
  217. * remplissage resultats ordre p
  218. *
  219. talam . p = lamdai ;
  220. tadep . p = lamdai / lambda1 * dep1 + u2;
  221. tasig . p = tadep . p 'SIGMA' 'LINE' mo ma + sig;
  222. taksi . p = 'KSIGMA' tasig . p mo ma;
  223. tagrad . p = tadep . p 'GRAD' mo ma;
  224. 'SI' ICOQ;tagraf . p = tadep . p 'GRAF' mo ma;'FINSI';
  225. * petit calcul pour voir l'evolution du rayon de convergence
  226. * norU1 = ('XTX' tadep. 1)**0.5D0 ;
  227. * norUn = ('XTX' tadep. P)**0.5D0 ;
  228. * aserie = ( prec * norU1 / norUn ) ** ( 1.D0 / ( p - 1 ) );
  229. * mess ' ordre ' p ' lambda ' lamdai 'raymax' aserie;
  230. *
  231. **Orthogonalisation de la base pour le calcul des approximants de Pade**
  232. *
  233. *** tabalpha. p = table;
  234. *calcul somme de (Up.Wq).Wq
  235. *** 'REPE' ort (p - 1 );
  236. *** alphapq = (tadep. p 'XTY' tabase. &ort llld llld );
  237. *** 'SI' ('EGA' &ort 1);
  238. *** vec = tadep. p - (alphapq*tabase. &ort) ;
  239. *** 'SINON' ;
  240. *** vec = vec - (alphapq*tabase. &ort);
  241. *** 'FINSI';
  242. *remplissage de la table contenant les alphapq issu de la projection
  243. *** tabalpha. p . &ort = alphapq;
  244. *** 'FIN' ort ;
  245.  
  246. *calcul des vecteurs de la nouvelle base
  247. *** Norvec = ('XTX' vec)**0.5D0 ;
  248. *** base = vec / Norvec ;
  249. *** tabalpha. p . p = norvec ;
  250. *remplissage de la table des vecteurs de la base orthonormee
  251. *** tabase. p = base;
  252. *temps place;
  253. 'FIN' bordre;
  254.  
  255.  
  256. **********************************************************************
  257. *************calcul du rayon de convergence maximal 'aserie'**********
  258. **********************************************************************
  259.  
  260. norU1 = ('XTX' tadep. 1)**0.5D0 ;
  261. norUn = ('XTX' tadep. N)**0.5D0;
  262. aserie = ( prec * norU1 / norUn ) ** ( 1.D0 / ( N - 1 ) );
  263. *
  264. **menage de la place memoire
  265. *
  266. 'REPE' titi ( 'DIME' taksi );
  267. detr taksi. &titi;
  268. 'FIN' titi;
  269. 'OUBLI' taksi;
  270. 'REPE' titi ( 'DIME' tagrad );
  271. 'DETR' tagrad. &titi;
  272. 'SI' ICOQ; 'DETR' tagraf. &titi;'FINSI';
  273. 'FIN' titi;
  274. 'OUBLI' tagrad; 'OUBLI' tagraf; 'OUBLI' gradi; 'OUBLI' gradj;
  275.  
  276. *
  277. ****Calcul de dn (coeff de Dn) en resolvant un systeme triangulaire****
  278. *
  279.  
  280. ***'REPE' dn (N - 1);
  281. *** tabvald. &dn = table;
  282. *** 'REPE' equa &dn;
  283. *** addi = tabalpha. (&dn + 1) . ( &dn + 1 - &equa );
  284. *** 'SI' ( &equa 'NEG' 1 ) ;
  285. *** 'REPE' ilig (&equa - 1);
  286. *** ope = tabalpha. ( &dn + 1 - &ilig ) . (&dn - &equa + 1) ;
  287. *** mult = ope * tabvald. &dn. &ilig ;
  288. *** addi = mult + addi ;
  289. *** 'FIN' ilig;
  290. *** 'FINSI';
  291. *** den = tabalpha. (&dn + 1 - &equa). (&dn + 1 - &equa);
  292. *** tabvald. &dn . &equa = (-1. * addi)/ den;
  293. *** 'FIN' equa ;
  294. ***'FIN' dn;
  295.  
  296. ******************************************************************
  297. *****Calcul du critere pour les approximants de Pade 'apade'******
  298. ****** Recherche de la valeur apade ***************
  299. ******************************************************************
  300. adep = aserie*3.;
  301. acou= adep;
  302.  
  303. ainf = 0.D0;
  304. recher = faux;
  305. epsipade =prec;
  306. puis2=1;
  307. **'REPE' cherchea ;
  308. **************creation d'une table puissance de a *************
  309. tapuisa = table;
  310.  
  311. 'REPE' puis (N - 1);
  312. puisacou = acou**&puis;
  313. tapuisa. &puis = puisacou;
  314. 'FIN' puis;
  315.  
  316. ********
  317. **********Calcul des polynomes Dn-1(a) pour chaque ordre***************
  318. ********
  319. *** tabDn = table;
  320. **la 1er valeur de cette table vaut 1**
  321. *** tabDn. 0 = 1.D0;
  322. *** adpoly= 1 ;
  323. *** 'REPE' polyDn (N - 1);
  324. *** adpoly = (tapuisa. &polyDn*tabvald.(N - 1). &polyDn) + adpoly;
  325. *** tabdn. &polyDn = adpoly;
  326. *** 'FIN' polyDn;
  327. **
  328. ******************Expression des approximants de Pade***************
  329. **
  330. ************Calcul de Pn(U(a)) - Pn-1 (U(a))*************
  331.  
  332. **Expression de Pn(U(a))
  333. *** PnUa = U0;
  334. * mess ' avant boucle boupa'; temps place;
  335. *** 'REPE' boupa (N - 1);
  336. *** som=tapuisa.&boupa/tabDn.(N - 1)*tadep. &boupa*
  337. *** tabDn.(N - 1 - &boupa);
  338. *** PnUa = som + PnUa ;
  339. *** 'FIN' boupa;
  340. * mess ' apres boucle pnua'; temps place;
  341. **Expression de Pn-1 (U(amp))
  342. *** Pmoins1 = U0;
  343. *** 'REPE' boupm (N - 2);
  344. *** sam=tapuisa. &boupm/tabDn.(N - 2)*tadep.&boupm*
  345. *** tabDn.(N - &boupm - 2);
  346. *** Pmoins1 = sam + Pmoins1 ;
  347. *** 'FIN' boupm;
  348. *****Calcul de (Pn(a) - Pn-1 (a)) / ( Pn(a) - U0)***
  349. *calcul norme de (Pn(a) - Pn-1 (a) )
  350. *** nt = PnUa - Pmoins1 ;
  351. *** nornt = (XTX nt)**0.5D0;
  352. *calcul norme de ( Pn(a) - U0)
  353. *** np = PnUa - U0;
  354. *** nornp = (XTX np)**0.5D0;
  355. *** crit = nornt / nornp ;
  356. **recherche de apade par dichotomie
  357. *** 'SI' ( 'EGA' recher faux);
  358. *** 'SI' (crit < epsipade);
  359. *** ainf = acou;
  360. *** acou = acou + adep;
  361. *** 'ITER' cherchea;
  362. *** 'SINON';
  363. *** asup = acou;
  364. *** acou = ( ainf + asup )/2;
  365. *** recher = vrai;
  366. *** compteur = 1;
  367. *** 'FINSI';
  368. *** 'SINON';
  369. *** compteur = compteur + 1;
  370. *** 'SI' (crit '>' epsipade);
  371. *** asup = acou;
  372. *** acou = (ainf + asup ) /2.D0;
  373. *** 'SINON';
  374. *** ainf = acou;
  375. *** acou = (asup + ainf)/2.D0;
  376. *** 'FINSI';
  377. *** 'SI' ( compteur '>' 16) ;
  378. *** apade = ainf;
  379. *** 'QUITTER' cherchea;
  380. *** 'SINON';
  381. *** 'ITERER' cherchea;
  382. *** 'FINSI';
  383. *** 'FINSI';
  384. ***'FIN' cherchea;
  385. *mess ' aserie ' aserie ' apade' apade ' compteur ' compteur;
  386. **
  387. ********************Resolution de Dn-1 = 0*********************
  388. **
  389. **Creation d'un listreel contenant les valeurs des coeff de Dn-1 **
  390. ***lll= prog 1.;
  391. ***'REPE' resol (N - 1) ;
  392. *** lla= prog tabvald. (N - 1) . &resol ;
  393. *** lll = lll et lla;
  394. ***
  395. ***'FIN' resol;
  396.  
  397. ****Resolution ******
  398. ***solpet = 'RACP' 'INTE' 0.D0 apade lll ;
  399. ***Recherche de la plus petite racine reelle positive***
  400. ***solpet = 1.e10;
  401. ***'SI'( ('DIME' soluu) > 0) ;
  402. *** 'REPE' ppsol ('DIME' soluu) ;
  403. *** aDn = 'EXTR' soluu &ppsol ;
  404. *** 'SI' ( aDn '>' 0.D0);
  405. *** 'SI' (aDn '<' solpet) ;
  406. *** solpet = aDn;
  407. *** 'FINSI';
  408. *** 'QUIT' ppsol;
  409. *** 'FINSI';
  410. *** 'FIN' ppsol ;
  411. ***'FINSI';
  412. ***
  413. ***
  414. *mess ' solpet ' solpet;
  415. ***'SI' ( 'NEG' solpet 'VIDE');
  416. *** 'SI' ( solpet '<' apade);
  417. *** apade =solpet ;
  418. ***'FINSI';
  419.  
  420.  
  421. *mess 'aDn ' solpet ' apade' apade 'aserie' aserie ;
  422.  
  423. ********************************************************************
  424. ********Recherche de 'a' pour que lambda vaille 1 pour les 2 cas****
  425. ********************************************************************
  426.  
  427.  
  428. ***xpetpa = apade;lamxpet=0.;
  429. ifini=faux;
  430. * mess 'Recherche de la valeur de a pour que lambda vaille 1: xpet' ;
  431. prec1= prec * 0.1;
  432. ifiser=faux;ifipad=faux;
  433. *'SI'( 'NON' PADE);
  434.  
  435. ************Pour le developpement en serie**************************
  436.  
  437. * recherche pour savoir si on croise lambda=1 Pour cela on evalue d'abord la
  438. * valeur obtenu pour aserie et si dlam + lam0 >1 on recherche les
  439. *
  440. * zeros du polynome donnant lambda
  441. * talam contient les coeff du polynoem sauf talam .0
  442. xnew = lamreel0;
  443. 'REPE' valr N;
  444. apuisn = aserie ** &valr;
  445. xnew=apuisn*talam.&valr + xnew;
  446. 'FIN' valr;
  447. * mess ' xnew ' xnew;
  448. 'SI' ('>' xnew 1.) ;
  449. lll = prog -1.D0 ;
  450. 'REPE' bobo N ;
  451. lla = prog talam . &bobo;
  452. lll= lll et lla;
  453. 'FIN' bobo;
  454. * 'LIST' lll;;
  455. * xsol= racp lll ;
  456. xpet= racp 'INTE' 0.D0 aserie lll;
  457. * mess ' c est bon xpet ' xpet ;
  458. ifiser=vrai;
  459. * list xraci;
  460. * recherche de la plus petite racine relle positive
  461. * 'SI' ( ('DIME' xsol ) > 0);
  462. * 'REPE' bobo ( 'DIME' xsol);
  463. * ab = 'EXTR' xsol &bobo ;
  464. * mess ' racine ' ab;
  465. * 'SI' ( ab '>' 0.D0 );
  466. * 'SI' ( ab '<' autil) ;
  467. * xpet = ab;
  468. * list ( xpet - xraci);
  469. * ifiser=vrai;
  470. * 'FINSI';
  471. ** 'QUIT' bobo;
  472. * 'FINSI';
  473. * 'FIN' bobo;
  474. * 'FINSI';
  475. 'SINON';
  476. xpet = aserie;
  477. * mess ' il faut recommencer xpet ' xpet;
  478. 'FINSI';
  479. *'SINON' ;
  480. *******************Pour les approximants de Pade********************
  481. * on ne fait ce travail que si le dev en serie n'arrive pas a la solution
  482. *** 'SI' ('NON' ifiser);
  483. *** ttpuis = table;
  484. *** solution = faux;
  485. *** xpetinf = 0.;
  486. *** lamrech= 1. ;
  487. *** 'REPE' solua ;
  488. *** 'REPE' pui (N - 1);
  489. *** puisxpet = xpetpa**&pui;
  490. *** ttpuis. &pui = puisxpet;
  491. *** 'FIN' pui;
  492. **********Calcul des polynomes Dn-1(asol)**********
  493. *** Dnxpet = table;
  494. *** Dnxpet. 0 = 1.;
  495. *** solpoly= 1 ;
  496.  
  497. *** 'REPE' polysol (N - 1);
  498. *** solpoly = (ttpuis. &polysol*tabvald.(N - 1). &polysol) + solpoly;
  499. *** Dnxpet. &polysol = solpoly;
  500. *** 'FIN' polysol;
  501.  
  502. *********Expression de Pn(lambda(asol))*************
  503. *** lamxpet = lam0;
  504.  
  505. *** 'REPE' sola (N - 1);
  506. *** sim=ttpuis.&sola/Dnxpet.(N - 1)*talam. &sola*
  507. *** Dnxpet.(N - 1 - &sola);
  508. *** lamxpet = sim + lamxpet ;
  509. *** 'FIN' sola;
  510. * 'MESS' ' solua ' &solua ' lamxpet' lamxpet;
  511. ********Resolution de lamda egale a 1********
  512. *** 'SI' ( 'EGA' solution faux);
  513. *** 'SI' (lamxpet '<' lamrech);
  514. *** ifpad=faux;
  515. *** 'QUITTER'solua;
  516. *** 'SINON';
  517. *** ifipad=vrai;
  518. ****** xpetsup = xpetpa;
  519. *** xpetpa = ( xpetinf + xpetsup )/2;
  520. *** solution = vrai;
  521. *** tour = 1;
  522. *** 'FINSI';
  523. *** 'SINON';
  524. *** tour = tour + 1;
  525. *** 'SI' (lamxpet > lamrech);
  526. *** xpetsup = xpetpa;
  527. *** xpetpa = (xpetinf + xpetsup ) / 2;
  528. *** 'SINON';
  529. *** xpetinf = xpetpa;
  530. *** xpetpa = (xpetsup + xpetinf)/2;
  531. *** 'FINSI';
  532. *** 'SI' ( (lamxpet '>' (lamrech - prec1) ) 'ET'
  533. *** (lamxpet '<' (lamrech + prec1) )) ;
  534. *** xpetpa = xpetinf;
  535. *** 'QUITTER' solua;
  536. *** 'SINON';
  537. *** 'ITERER' solua;
  538. *** 'FINSI';
  539. *** 'FINSI';
  540. * mess ' nombre de tour' tour;
  541. *** 'FIN' solua;
  542. *** 'FINSI';
  543.  
  544. *'FINSI';
  545. * mess 'serie ' xpet 'apade' lamxpet;
  546. *
  547. *
  548. * choix de la methode
  549. *
  550. PADE=FAUX;
  551. 'SI' ifiser ;
  552. *** pade=faux;
  553. ifini=vrai;
  554. *** 'SINON';
  555. * 'SI' ifipad;
  556. * ifini=vrai;
  557. * xpet=xpetpa;
  558. * 'SINON';
  559. * 'SI' ( xnew '>' lamxpet );
  560. * pade=faux;
  561. * 'SINON';
  562. * xpet=xpetpa;
  563. * 'FINSI';
  564. * 'FINSI';
  565. *** pade=faux;
  566. 'FINSI';
  567.  
  568. **creation de table pour les differentes valeurs de a****
  569. valeura = table;
  570. dray = xpet ;
  571. valeura. 1 = dray;
  572.  
  573. ********************************************************************
  574. ******************Calcul de lambda et de U *************************
  575. ************** par les approximants de Pade *********************
  576. *********** et par le developpement en series ******************
  577. *******************************************************************
  578.  
  579. *********** par le developpement en series ******************
  580.  
  581. ***'SI' ( 'NON' pade);
  582.  
  583. 'REPE' valree N;
  584. apuisn = valeura.1 ** &valree;
  585. depreel0=apuisn*tadep.&valree + depreel0;
  586. deltau0=apuisn*tadep.&valree + deltau0;
  587. sigreel0=apuisn*tasig.&valree + sigreel0;
  588. lamreel0=apuisn*talam.&valree + lamreel0;
  589. 'FIN' valree;
  590. ************** par les approximants de Pade *********************
  591. ***'SINON';
  592. *** tabamax = table;
  593. *** ttDn = table;
  594. *** tabPn = table;
  595. *** tabamax. 1 = table;
  596. *** 'REPE' amax (N - 1);
  597. *** puisamax = valeura. 1**&amax;
  598. *** tabamax. &amax = puisamax;
  599. *** 'FIN' amax;
  600. *** ttDn. 0 = 1.;
  601. *** adpolyn= 1 ;
  602.  
  603. **Expression des Dn-1 *********
  604. *** 'REPE' Dnmax (N - 1);
  605. *** eta = tabamax. &Dnmax * tabvald. (N - 1). &Dnmax;
  606. *** adpolyn = eta + adpolyn;
  607. *** ttDn. &Dnmax = adpolyn;
  608. *** 'FIN' Dnmax;
  609. * mess ' apres boucle dnmax'; temps place;
  610. *** 'REPE' padeexp (N - 1);
  611. **Expression de Pn(U(a))********
  612. *** exp = ttDn. (N - 1 - &padeexp)/ttDn. (N - 1);
  613. *** PnUa=exp*tabamax. &padeexp*tadep. &padeexp;
  614. *** depreel0 = PnUa + depreel0 ;
  615. *** deltau0=PnUa + deltau0;
  616. **Expression de lambda(U(a))****
  617. *** Pnlama=exp*tabamax. &padeexp*talam. &padeexp;
  618. *** lamreel0 = Pnlama + lamreel0 ;
  619. **Expression de sigma(U(a))****
  620. *** Pnsiga = exp*tabamax. &padeexp*tasig. &padeexp;
  621. *** sigreel0 = Pnsiga + sigreel0;
  622. *** 'FIN' padeexp;
  623. * mess ' apres padeexp' ; temps place;
  624. ***'FINSI';
  625.  
  626. ttdep. 1 = depreel0;
  627.  
  628. ttsig. 1 = sigreel0 'PICA' ( depreel0 - U0 ) mo ;
  629. tsetse = ttsig. 1 ;
  630. ttlam. 1 = lamreel0;
  631. *
  632. *test pour savoir si on repart dans le bon sens
  633. *
  634. dlamda = ttlam . 1 - ttlam . (1 - 1) ;
  635. dlam = dlamda /dray ;
  636. sgd = signe( dlam) ; sgd= sgdlamda;
  637. 'SI' (sgdlamda 'NEG' sgd ) ;
  638. MESS ' On doit rediriger le pas ' ;
  639. 'REPE' bouval1 1;
  640. valeura. &bouval1 = dray * &bouval1 * (-1) ;
  641. TAB . 'A' . ( &bouval1 + NA) = valeura. &bouval1 +a0 ;
  642. 'FIN' bouval1;
  643. ********** par le developpement en series ******************
  644. * list valeura ;
  645. *** 'SI' ( 'NON' pade);
  646. depreel0 = depreel0 'ENLE' LX;
  647. 'REPE' valree1 N ;
  648. apuisn = valeura.1 ** &valree1;
  649. depreel0=apuisn*tadep.&valree1 + depreel0 ;
  650. sigreel0=apuisn*tasig.&valree1 + sigreel0 ;
  651. lamreel0=apuisn*talam.&valree1 + lamreel0 ;
  652. 'FIN' valree1;
  653.  
  654. ************** par les approximants de Pade *********************
  655. *** 'SINON';
  656. *** tabamax = table;
  657. *** ttDn = table;
  658. *** tabPn = table;
  659. *** tabamax. 1 = table;
  660. *** 'REPE' am1 (N - 1);
  661. *** puisamax = valeura. 1**&am1;
  662. *** tabamax. &am1 = puisamax;
  663. *** 'FIN' am1;
  664. *** ttDn. 0 = 1.;
  665. *** adpolyn= 1 ;
  666. **Expression des Dn-1 *********
  667. *** 'REPE' Dn1 (N - 1);
  668. *** eta = tabamax. &Dn1 * tabvald. (N - 1). &Dn1;
  669. *** adpolyn = eta + adpolyn;
  670. *** ttDn. &Dn1 = adpolyn;
  671. *** 'FIN' Dn1;
  672. *** 'REPE' pade1 (N - 1);
  673. **Expression de Pn(U(a))********
  674. *** exp = ttDn. (N - 1 - &pade1)/ttDn. (N - 1);
  675. *** PnUa=exp*tabamax. &pade1*tadep. &pade1;
  676. *** depreel0 = PnUa + depreel0 ;
  677. **Expression de lambda(U(a))****
  678. *** Pnlama=exp*tabamax. &pade1*talam. &pade1;
  679. *** lamreel0 = Pnlama + lamreel0 ;
  680. **Expression de sigma(U(a))****
  681. *** Pnsiga = exp*tabamax. &pade1*tasig. &pade1;
  682. *** sigreel0 = Pnsiga + sigreel0;
  683. *** 'FIN' pade1;
  684. *** 'FINSI';
  685. dlamda = ttlam . 1 - ttlam . (1 - 1) ;
  686. dlam = dlamda /dray * (-1 ) ;
  687. sgd = signe( dlam) ;
  688. *
  689. * MESS ' valeur de dlam aprés redirection ' dlam ;
  690. * MESS ' LE SIGNE après redirection ' SGd ;
  691. *
  692. ttdep . 1 = depreel0 ;
  693. ttsig . 1 = sigreel0 PICA ( depreel0 - U0 ) mo ;
  694. tsetse=ttsig . 1 ;
  695. ttlam . 1 = lamreel0;
  696. * U0= depreel0 ;
  697. 'FINSI';
  698. sigreel0=sigreel0 'PICA' deltau0 mo ;
  699. deltatot=deltau0 + deltatot ;
  700. sigpiok= ('SIGMA' deltatot mo ma) + S0;
  701. sigcau= sigpiok 'PICA' deltatot mo;
  702. sigreel0= sigcau;
  703. ttdep.0 = depreel0;
  704. ttsig.0 = S0;
  705. ttlam.0= lamreel0;
  706. 'SI' ('NON' ifini) ;
  707. *mess ' lamrell0 ' lamreel0;
  708. ff= ff *( 1.d0 - lamreel0);
  709. depreel0= deltatot 'ENLE' 'LX';
  710. lamreel0=0.D0;
  711. 'FORM' geoini;
  712. 'FORM' deltatot;
  713. U0= deltatot;
  714. 'FINSI';
  715. *a0 = tab . a . (1 + na) ;
  716. dlamda = (ttlam . 1 - ttlam . 0 );
  717. dlamb = dlamda /dray ;
  718. sgdlamda = signe( dlamb) ;
  719. *tab.'SIGNE'=sgdlamda;
  720. 'FORM' geoini;
  721. 'FORM' deltatot;
  722. 'SI' ifini;
  723. 'MESS' ' sortie de la MAN nombre de segments ' &it IOUT;
  724. form geoini;
  725. 'QUITTER' IT;
  726. 'FINSI';
  727. 'FIN' IT;
  728. 'FINPROC' deltatot IOUT;
  729.  
  730.  
  731.  

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