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