$$$$ CORMAN * CORMAN PROCEDUR CHAT 13/02/05 21:15:00 7691 'DEBPROC' corman rit*'RIGIDITE' MO*'MMODEL' MA*'MCHAML' ORDRE*'ENTIER' U0*'CHPOINT' S0*'MCHAML' ffi*CHPOINT ETAB12*'TABLE'; * quan,d on arrive ici on se trouve sur la config en debut * de pas. ffi est le résidu et s0 est la contrainte en début de pas * * * attention dans cette version on n'utilise pas les approximants de Padé. * sgdlamda='SIGN' 1 ; *** mess ' entree dans corman ' ; IDYN=ETAB12.'DYNAMIQUE'; 'SI' IDYN; MMA= ETAB12. 'MASSE'; QUADTT2=4.D0 / ETAB12.'DT' / ETAB12.'DT'; MMAPLU=MMA * QUADTT2; * mess ' quadtt2 ' quadtt2; 'FINSI'; N=ORDRE ; DEP1=U0; ncontmax=250; isol1=faux; prec = 1.D-8; ttdep=table; ttdep.0 = U0; ttlam=table ; ttlam.0=0.; ttsig=table; ttsig.0=s0;NA=0 ; depreel0 = U0 'ENLE' 'LX' ; depini=U0 'ENLE' 'LX' ; sigreel0 = S0 ; ffs0= 'BSIGMA' S0 MO MA; lamreel0 = 0.d0; IOUT=1; ritblo='EXTR' rit 'RIGI' 'MULT'; * on veut lambda=1.d0 llld = mots 'UX' 'UY' 'UZ' 'RX' 'RY' 'RZ' ; IMAS=FAUX; IDIMV= 'VALE' 'DIME'; IMODD= 'VALE' 'MODE'; 'SI' ( 'EXISTE' MO 'ELEM' 'TRI3' 'QUA4' 'TRI6' 'QUA8' 'CUB8' 'TET4' 'PRI6' 'PYR5' 'CU20' 'TE10' 'PR15' 'PY13' ) ; IMAS=VRAI;MOMAS= 'EXTR' MO 'ELEM''TRI3' 'QUA4' 'TRI6' 'QUA8' 'CUB8' 'TET4' 'PRI6' 'PYR5' 'CU20' 'TE10' 'PR15' 'PY13'; 'FINSI'; ICOQ=FAUX; 'SI' ( 'EXISTE' MO 'ELEM' 'COQ2' 'DKT' );ICOQ=VRAI; MOCOQ= 'EXTR' MO 'ELEM' 'COQ2' 'DKT'; 'FINSI'; 'SI' ('EXIS' MO ELEM 'POUT'); IOUT=0;'QUITTER' CORMAN; 'FINSI'; ff = ffi ; lam0 = 0. ; geoini = form; sigcour=s0; *m0 deltatot= u0*0.; 'REPE' IT ; 'SI' ('EXIS' MO ELEM 'POUT'); IOUT=0;'QUITTER' IT; 'FINSI'; 'SI' ('EXIS' MO ELEM 'COQ2'); IOUT=0;'QUITTER' IT; 'FINSI'; *'SI' (&it > 1) ; mess ' continuation ' ( &it - 1) lam0;'FINSI'; 'SI' (&it > ncontmax);MESS ' Pas de solution apres ' ncontmax 'continuations'; IOUT=0; 'QUITTER' IT; 'FINSI'; **************************Initialisation************************** * * quelques initialisation de table servant aux stockages intermediaire * tadep=table; talam=table;tasig=table;taksi = table; tagrad = table; tagraf=table; tabase = table;tabalpha = table; tabvald = table; tabpoly = table; * ************************** calcul du premier ordre****************** * deltau0= u0*0.; kspp = 'KSIGMA' sigreel0 mo ma; rit= 'RIGI' MO MA; rit=rit 'ET' ritblo; ritot = rit 'ET' kspp; 'SI' IDYN; Ritot= ritot 'ET' MMAPLU; 'FINSI'; DE1= 'RESOU' ritot ff 'NOUNIL'; psc = 'XTX' de1 ; psc2 = psc + 1.D0 ** 0.5D0; lambda1 = 1.d0 / psc2 ; dep1 = de1 * lambda1; s1 = 'SIGMA' 'LINE' dep1 mo ma; * * calcul d'une norme de convergence * * remplissage des tables pour ordre 1 tadep . 1 = dep1; talam . 1 = lambda1 ; tasig . 1 = s1; taksi . 1 = 'KSIGMA' s1 mo ma; tagrad . 1 = 'GRAD' dep1 mo ma; 'SI' ICOQ ;tagraf . 1 = 'GRAF' dep1 mo ma ;'FINSI'; tabalpha. 1 = table; tabalpha . 1 . 1 = ('XTX' tadep. 1)**0.5D0 ; tabase. 1 = (tadep. 1)/(tabalpha. 1 . 1); p=1; *mess ' ordre ' p ' lambda ' lambda1; *temps place; * ******************** calcul des ordres superieurs***************** * *mess ' debut de boucle'; 'REPETER' bordre (N - 1); p = &bordre + 1; * *calcul des ksig* (si) * u(p - i + 1); nfois = p - 1 ; 'REPETER' ks nfois; fksi = taksi . &ks * tadep .( p - &ks); 'SI' ( &ks 'EGA' 1); fks = -1.D0 * fksi; 'SINON'; *menagement de l'espace memoire XXX = fks; fks = fks - fksi; 'FINSI'; 'FIN' ks; fks = fks ; * * calcul de bsigma ( D *0.5( gradui * graduj)) * 'REPETER' bs nfois; gradi = tagrad . &bs; gradj = tagrad . ( p - &bs ); 'SI' ICOQ ; grafi = tagraf . &bs; grafj = tagraf . ( p - &bs ); 'FINSI'; * les cartes suivantes sont pour les elements coques 'SI' ICOQ ; gradi = 'REDU' gradi MOCOQ; gradj = 'REDU' gradj MOCOQ; grafi = 'REDU' grafi MOCOQ; grafj = 'REDU' grafj MOCOQ; Uxxi='EXCO' gradi 'UX,X' 'SCAL' chan TYPE 'SCAL'; Uxyi='EXCO' gradi 'UX,Y' 'SCAL' chan TYPE 'SCAL' ; Uxzi='EXCO' gradi 'UX,Z' 'SCAL' chan TYPE 'SCAL' ; Uyxi='EXCO' gradi 'UY,X' 'SCAL' chan TYPE 'SCAL' ; Uyyi='EXCO' gradi 'UY,Y' 'SCAL' chan TYPE 'SCAL' ; Uyzi='EXCO' gradi 'UY,Z' 'SCAL' chan TYPE 'SCAL' ; Uzxi='EXCO' gradi 'UZ,X' 'SCAL' chan TYPE 'SCAL' ; Uzyi='EXCO' gradi 'UZ,Y' 'SCAL' chan TYPE 'SCAL' ; Uxxj='EXCO' gradj 'UX,X' 'SCAL' chan TYPE 'SCAL' ; Uxyj='EXCO' gradj 'UX,Y' 'SCAL' chan TYPE 'SCAL' ; Uxzj='EXCO' gradj 'UX,Z' 'SCAL' chan TYPE 'SCAL' ; Uyxj='EXCO' gradj 'UY,X' 'SCAL' chan TYPE 'SCAL' ; Uyyj='EXCO' gradj 'UY,Y' 'SCAL' chan TYPE 'SCAL' ; Uyzj='EXCO' gradj 'UY,Z' 'SCAL' chan TYPE 'SCAL' ; Uzxj='EXCO' gradj 'UZ,X' 'SCAL' chan TYPE 'SCAL' ; Uzyj='EXCO' gradj 'UZ,Y' 'SCAL' chan TYPE 'SCAL' ; Bxxi='EXCO' grafi 'BX,X' 'SCAL' chan TYPE 'SCAL' ; Bxyi='EXCO' grafi 'BX,Y' 'SCAL' chan TYPE 'SCAL' ; Byxi='EXCO' grafi 'BY,X' 'SCAL' chan TYPE 'SCAL' ; Byyi='EXCO' grafi 'BY,Y' 'SCAL' chan TYPE 'SCAL' ; Bxxj='EXCO' grafj 'BX,X' 'SCAL' chan TYPE 'SCAL' ; Bxyj='EXCO' grafj 'BX,Y' 'SCAL' chan TYPE 'SCAL' ; Byxj='EXCO' grafj 'BY,X' 'SCAL' chan TYPE 'SCAL'; * Nouvelle version epss=( 0.5 * ( Uxxi*Uxxj + (Uyxi*Uyxj) + (Uzxi*Uzxj ) ) ) 'NOMC' 'EPSS'; eptt=( 0.5 * ( Uxyi*Uxyj + (Uyyi*Uyyj) + (Uzyi*Uzyj ) ) ) 'NOMC' 'EPTT'; gast=(0.5 * (Uxyi*Uxxj + (Uyyi*Uyxj) + (Uzyi*Uzxj) + (Uxyj*Uxxi) + (Uyyj*Uyxi) + (Uzxi*Uzyj) ) )'NOMC' 'GAST'; rtss=(epss *0.) 'NOMC' 'RTSS'; rttt=rtss 'NOMC' 'RTTT'; rtst=rtss 'NOMC' 'RTST'; ep2 = epss + eptt + gast + rtss + rttt + rtst; 'FINSI'; * les cartes suivantes sont pour les elements massifs 'SI' IMAS ; gradi = 'REDU' gradi MOMAS;gradj = 'REDU' gradj MOMAS; gru2 =(( 'TAGR' gradi) * gradj mo) * 0.5D0 ; 'SI' ('NEG' imodd 'AXIS'); epxyz= 'EXCO' gru2 ('MOTS' 'UX,X' 'UY,Y' 'UZ,Z') ('MOTS' 'EPXX' 'EPYY' 'EPZZ'); gaxyz= ('EXCO' gru2 ('MOTS' 'UX,Y' 'UX,Z' 'UY,Z') ('MOTS' 'GAXY' 'GAXZ' 'GAYZ'))*2.d0; ep2 = epxyz + gaxyz; 'SINON'; epxx = 'EXCO' gru2 UR,R 'EPRR'; epyy = 'EXCO' gru2 UZ,Z 'EPZZ'; epzz = 'EXCO' gru2 UT,T 'EPTT'; gaxy = ('EXCO' gru2 UR,Z 'GARZ')*2.D0 ; gaxz = ('EXCO' gru2 UR,T 'GART')*2.D0 ; gayz = ('EXCO' gru2 UT,Z 'GATZ')*2.D0 ; ep2 = epxx + epyy + epzz + gaxy + gaxz + gayz; 'FINSI'; 'FINSI'; * fin du traitement spécifique ep2 = 'CHANG' ep2 'TYPE' 'DEFORMATIONS'; 'SI' ( 'EGA' &bs 1); ep2tot = ep2 ; 'SINON'; XXX = ep2; ep2tot = ep2tot + ep2; 'FINSI'; 'FIN' bs; sig = 'ELAS' ep2tot mo ma ; fbs = ('BSIGMA' sig mo ma) * -1.D0 ; u2 = 'RESOU' ritot ( fks + fbs) 'NOUNIL'; lamdai = (u2 'XTY' tadep.1 llld llld )* -1.D0 * lambda1; * * * remplissage resultats ordre p * talam . p = lamdai ; tadep . p = lamdai / lambda1 * dep1 + u2; tasig . p = tadep . p 'SIGMA' 'LINE' mo ma + sig; taksi . p = 'KSIGMA' tasig . p mo ma; tagrad . p = tadep . p 'GRAD' mo ma; 'SI' ICOQ;tagraf . p = tadep . p 'GRAF' mo ma;'FINSI'; * petit calcul pour voir l'evolution du rayon de convergence * norU1 = ('XTX' tadep. 1)**0.5D0 ; * norUn = ('XTX' tadep. P)**0.5D0 ; * aserie = ( prec * norU1 / norUn ) ** ( 1.D0 / ( p - 1 ) ); * mess ' ordre ' p ' lambda ' lamdai 'raymax' aserie; * **Orthogonalisation de la base pour le calcul des approximants de Pade** * *** tabalpha. p = table; *calcul somme de (Up.Wq).Wq *** 'REPE' ort (p - 1 ); *** alphapq = (tadep. p 'XTY' tabase. &ort llld llld ); *** 'SI' ('EGA' &ort 1); *** vec = tadep. p - (alphapq*tabase. &ort) ; *** 'SINON' ; *** vec = vec - (alphapq*tabase. &ort); *** 'FINSI'; *remplissage de la table contenant les alphapq issu de la projection *** tabalpha. p . &ort = alphapq; *** 'FIN' ort ; *calcul des vecteurs de la nouvelle base *** Norvec = ('XTX' vec)**0.5D0 ; *** base = vec / Norvec ; *** tabalpha. p . p = norvec ; *remplissage de la table des vecteurs de la base orthonormee *** tabase. p = base; *temps place; 'FIN' bordre; ********************************************************************** *************calcul du rayon de convergence maximal 'aserie'********** ********************************************************************** norU1 = ('XTX' tadep. 1)**0.5D0 ; norUn = ('XTX' tadep. N)**0.5D0; aserie = ( prec * norU1 / norUn ) ** ( 1.D0 / ( N - 1 ) ); * **menage de la place memoire * 'REPE' titi ( 'DIME' taksi ); detr taksi. &titi; 'FIN' titi; 'OUBLI' taksi; 'REPE' titi ( 'DIME' tagrad ); 'DETR' tagrad. &titi; 'SI' ICOQ; 'DETR' tagraf. &titi;'FINSI'; 'FIN' titi; 'OUBLI' tagrad; 'OUBLI' tagraf; 'OUBLI' gradi; 'OUBLI' gradj; * ****Calcul de dn (coeff de Dn) en resolvant un systeme triangulaire**** * ***'REPE' dn (N - 1); *** tabvald. &dn = table; *** 'REPE' equa &dn; *** addi = tabalpha. (&dn + 1) . ( &dn + 1 - &equa ); *** 'SI' ( &equa 'NEG' 1 ) ; *** 'REPE' ilig (&equa - 1); *** ope = tabalpha. ( &dn + 1 - &ilig ) . (&dn - &equa + 1) ; *** mult = ope * tabvald. &dn. &ilig ; *** addi = mult + addi ; *** 'FIN' ilig; *** 'FINSI'; *** den = tabalpha. (&dn + 1 - &equa). (&dn + 1 - &equa); *** tabvald. &dn . &equa = (-1. * addi)/ den; *** 'FIN' equa ; ***'FIN' dn; ****************************************************************** *****Calcul du critere pour les approximants de Pade 'apade'****** ****** Recherche de la valeur apade *************** ****************************************************************** adep = aserie*3.; acou= adep; ainf = 0.D0; recher = faux; epsipade =prec; puis2=1; **'REPE' cherchea ; **************creation d'une table puissance de a ************* tapuisa = table; 'REPE' puis (N - 1); puisacou = acou**&puis; tapuisa. &puis = puisacou; 'FIN' puis; ******** **********Calcul des polynomes Dn-1(a) pour chaque ordre*************** ******** *** tabDn = table; **la 1er valeur de cette table vaut 1** *** tabDn. 0 = 1.D0; *** adpoly= 1 ; *** 'REPE' polyDn (N - 1); *** adpoly = (tapuisa. &polyDn*tabvald.(N - 1). &polyDn) + adpoly; *** tabdn. &polyDn = adpoly; *** 'FIN' polyDn; ** ******************Expression des approximants de Pade*************** ** ************Calcul de Pn(U(a)) - Pn-1 (U(a))************* **Expression de Pn(U(a)) *** PnUa = U0; * mess ' avant boucle boupa'; temps place; *** 'REPE' boupa (N - 1); *** som=tapuisa.&boupa/tabDn.(N - 1)*tadep. &boupa* *** tabDn.(N - 1 - &boupa); *** PnUa = som + PnUa ; *** 'FIN' boupa; * mess ' apres boucle pnua'; temps place; **Expression de Pn-1 (U(amp)) *** Pmoins1 = U0; *** 'REPE' boupm (N - 2); *** sam=tapuisa. &boupm/tabDn.(N - 2)*tadep.&boupm* *** tabDn.(N - &boupm - 2); *** Pmoins1 = sam + Pmoins1 ; *** 'FIN' boupm; *****Calcul de (Pn(a) - Pn-1 (a)) / ( Pn(a) - U0)*** *calcul norme de (Pn(a) - Pn-1 (a) ) *** nt = PnUa - Pmoins1 ; *** nornt = (XTX nt)**0.5D0; *calcul norme de ( Pn(a) - U0) *** np = PnUa - U0; *** nornp = (XTX np)**0.5D0; *** crit = nornt / nornp ; **recherche de apade par dichotomie *** 'SI' ( 'EGA' recher faux); *** 'SI' (crit < epsipade); *** ainf = acou; *** acou = acou + adep; *** 'ITER' cherchea; *** 'SINON'; *** asup = acou; *** acou = ( ainf + asup )/2; *** recher = vrai; *** compteur = 1; *** 'FINSI'; *** 'SINON'; *** compteur = compteur + 1; *** 'SI' (crit '>' epsipade); *** asup = acou; *** acou = (ainf + asup ) /2.D0; *** 'SINON'; *** ainf = acou; *** acou = (asup + ainf)/2.D0; *** 'FINSI'; *** 'SI' ( compteur '>' 16) ; *** apade = ainf; *** 'QUITTER' cherchea; *** 'SINON'; *** 'ITERER' cherchea; *** 'FINSI'; *** 'FINSI'; ***'FIN' cherchea; *mess ' aserie ' aserie ' apade' apade ' compteur ' compteur; ** ********************Resolution de Dn-1 = 0********************* ** **Creation d'un listreel contenant les valeurs des coeff de Dn-1 ** ***lll= prog 1.; ***'REPE' resol (N - 1) ; *** lla= prog tabvald. (N - 1) . &resol ; *** lll = lll et lla; *** ***'FIN' resol; ****Resolution ****** ***solpet = 'RACP' 'INTE' 0.D0 apade lll ; ***Recherche de la plus petite racine reelle positive*** ***solpet = 1.e10; ***'SI'( ('DIME' soluu) > 0) ; *** 'REPE' ppsol ('DIME' soluu) ; *** aDn = 'EXTR' soluu &ppsol ; *** 'SI' ( aDn '>' 0.D0); *** 'SI' (aDn '<' solpet) ; *** solpet = aDn; *** 'FINSI'; *** 'QUIT' ppsol; *** 'FINSI'; *** 'FIN' ppsol ; ***'FINSI'; *** *** *mess ' solpet ' solpet; ***'SI' ( 'NEG' solpet 'VIDE'); *** 'SI' ( solpet '<' apade); *** apade =solpet ; ***'FINSI'; *mess 'aDn ' solpet ' apade' apade 'aserie' aserie ; ******************************************************************** ********Recherche de 'a' pour que lambda vaille 1 pour les 2 cas**** ******************************************************************** ***xpetpa = apade;lamxpet=0.; ifini=faux; * mess 'Recherche de la valeur de a pour que lambda vaille 1: xpet' ; prec1= prec * 0.1; ifiser=faux;ifipad=faux; *'SI'( 'NON' PADE); ************Pour le developpement en serie************************** * recherche pour savoir si on croise lambda=1 Pour cela on evalue d'abord la * valeur obtenu pour aserie et si dlam + lam0 >1 on recherche les * * zeros du polynome donnant lambda * talam contient les coeff du polynoem sauf talam .0 xnew = lamreel0; 'REPE' valr N; apuisn = aserie ** &valr; xnew=apuisn*talam.&valr + xnew; 'FIN' valr; * mess ' xnew ' xnew; 'SI' ('>' xnew 1.) ; lll = prog -1.D0 ; 'REPE' bobo N ; lla = prog talam . &bobo; lll= lll et lla; 'FIN' bobo; * 'LIST' lll;; * xsol= racp lll ; xpet= racp 'INTE' 0.D0 aserie lll; * mess ' c est bon xpet ' xpet ; ifiser=vrai; * list xraci; * recherche de la plus petite racine relle positive * 'SI' ( ('DIME' xsol ) > 0); * 'REPE' bobo ( 'DIME' xsol); * ab = 'EXTR' xsol &bobo ; * mess ' racine ' ab; * 'SI' ( ab '>' 0.D0 ); * 'SI' ( ab '<' autil) ; * xpet = ab; * list ( xpet - xraci); * ifiser=vrai; * 'FINSI'; ** 'QUIT' bobo; * 'FINSI'; * 'FIN' bobo; * 'FINSI'; 'SINON'; xpet = aserie; * mess ' il faut recommencer xpet ' xpet; 'FINSI'; *'SINON' ; *******************Pour les approximants de Pade******************** * on ne fait ce travail que si le dev en serie n'arrive pas a la solution *** 'SI' ('NON' ifiser); *** ttpuis = table; *** solution = faux; *** xpetinf = 0.; *** lamrech= 1. ; *** 'REPE' solua ; *** 'REPE' pui (N - 1); *** puisxpet = xpetpa**&pui; *** ttpuis. &pui = puisxpet; *** 'FIN' pui; **********Calcul des polynomes Dn-1(asol)********** *** Dnxpet = table; *** Dnxpet. 0 = 1.; *** solpoly= 1 ; *** 'REPE' polysol (N - 1); *** solpoly = (ttpuis. &polysol*tabvald.(N - 1). &polysol) + solpoly; *** Dnxpet. &polysol = solpoly; *** 'FIN' polysol; *********Expression de Pn(lambda(asol))************* *** lamxpet = lam0; *** 'REPE' sola (N - 1); *** sim=ttpuis.&sola/Dnxpet.(N - 1)*talam. &sola* *** Dnxpet.(N - 1 - &sola); *** lamxpet = sim + lamxpet ; *** 'FIN' sola; * 'MESS' ' solua ' &solua ' lamxpet' lamxpet; ********Resolution de lamda egale a 1******** *** 'SI' ( 'EGA' solution faux); *** 'SI' (lamxpet '<' lamrech); *** ifpad=faux; *** 'QUITTER'solua; *** 'SINON'; *** ifipad=vrai; ****** xpetsup = xpetpa; *** xpetpa = ( xpetinf + xpetsup )/2; *** solution = vrai; *** tour = 1; *** 'FINSI'; *** 'SINON'; *** tour = tour + 1; *** 'SI' (lamxpet > lamrech); *** xpetsup = xpetpa; *** xpetpa = (xpetinf + xpetsup ) / 2; *** 'SINON'; *** xpetinf = xpetpa; *** xpetpa = (xpetsup + xpetinf)/2; *** 'FINSI'; *** 'SI' ( (lamxpet '>' (lamrech - prec1) ) 'ET' *** (lamxpet '<' (lamrech + prec1) )) ; *** xpetpa = xpetinf; *** 'QUITTER' solua; *** 'SINON'; *** 'ITERER' solua; *** 'FINSI'; *** 'FINSI'; * mess ' nombre de tour' tour; *** 'FIN' solua; *** 'FINSI'; *'FINSI'; * mess 'serie ' xpet 'apade' lamxpet; * * * choix de la methode * PADE=FAUX; 'SI' ifiser ; *** pade=faux; ifini=vrai; *** 'SINON'; * 'SI' ifipad; * ifini=vrai; * xpet=xpetpa; * 'SINON'; * 'SI' ( xnew '>' lamxpet ); * pade=faux; * 'SINON'; * xpet=xpetpa; * 'FINSI'; * 'FINSI'; *** pade=faux; 'FINSI'; **creation de table pour les differentes valeurs de a**** valeura = table; dray = xpet ; valeura. 1 = dray; ******************************************************************** ******************Calcul de lambda et de U ************************* ************** par les approximants de Pade ********************* *********** et par le developpement en series ****************** ******************************************************************* *********** par le developpement en series ****************** ***'SI' ( 'NON' pade); 'REPE' valree N; apuisn = valeura.1 ** &valree; depreel0=apuisn*tadep.&valree + depreel0; deltau0=apuisn*tadep.&valree + deltau0; sigreel0=apuisn*tasig.&valree + sigreel0; lamreel0=apuisn*talam.&valree + lamreel0; 'FIN' valree; ************** par les approximants de Pade ********************* ***'SINON'; *** tabamax = table; *** ttDn = table; *** tabPn = table; *** tabamax. 1 = table; *** 'REPE' amax (N - 1); *** puisamax = valeura. 1**&amax; *** tabamax. &amax = puisamax; *** 'FIN' amax; *** ttDn. 0 = 1.; *** adpolyn= 1 ; **Expression des Dn-1 ********* *** 'REPE' Dnmax (N - 1); *** eta = tabamax. &Dnmax * tabvald. (N - 1). &Dnmax; *** adpolyn = eta + adpolyn; *** ttDn. &Dnmax = adpolyn; *** 'FIN' Dnmax; * mess ' apres boucle dnmax'; temps place; *** 'REPE' padeexp (N - 1); **Expression de Pn(U(a))******** *** exp = ttDn. (N - 1 - &padeexp)/ttDn. (N - 1); *** PnUa=exp*tabamax. &padeexp*tadep. &padeexp; *** depreel0 = PnUa + depreel0 ; *** deltau0=PnUa + deltau0; **Expression de lambda(U(a))**** *** Pnlama=exp*tabamax. &padeexp*talam. &padeexp; *** lamreel0 = Pnlama + lamreel0 ; **Expression de sigma(U(a))**** *** Pnsiga = exp*tabamax. &padeexp*tasig. &padeexp; *** sigreel0 = Pnsiga + sigreel0; *** 'FIN' padeexp; * mess ' apres padeexp' ; temps place; ***'FINSI'; ttdep. 1 = depreel0; ttsig. 1 = sigreel0 'PICA' ( depreel0 - U0 ) mo ; tsetse = ttsig. 1 ; ttlam. 1 = lamreel0; * *test pour savoir si on repart dans le bon sens * dlamda = ttlam . 1 - ttlam . (1 - 1) ; dlam = dlamda /dray ; sgd = signe( dlam) ; sgd= sgdlamda; 'SI' (sgdlamda 'NEG' sgd ) ; MESS ' On doit rediriger le pas ' ; 'REPE' bouval1 1; valeura. &bouval1 = dray * &bouval1 * (-1) ; TAB . 'A' . ( &bouval1 + NA) = valeura. &bouval1 +a0 ; 'FIN' bouval1; ********** par le developpement en series ****************** * list valeura ; *** 'SI' ( 'NON' pade); depreel0 = depreel0 'ENLE' LX; 'REPE' valree1 N ; apuisn = valeura.1 ** &valree1; depreel0=apuisn*tadep.&valree1 + depreel0 ; sigreel0=apuisn*tasig.&valree1 + sigreel0 ; lamreel0=apuisn*talam.&valree1 + lamreel0 ; 'FIN' valree1; ************** par les approximants de Pade ********************* *** 'SINON'; *** tabamax = table; *** ttDn = table; *** tabPn = table; *** tabamax. 1 = table; *** 'REPE' am1 (N - 1); *** puisamax = valeura. 1**&am1; *** tabamax. &am1 = puisamax; *** 'FIN' am1; *** ttDn. 0 = 1.; *** adpolyn= 1 ; **Expression des Dn-1 ********* *** 'REPE' Dn1 (N - 1); *** eta = tabamax. &Dn1 * tabvald. (N - 1). &Dn1; *** adpolyn = eta + adpolyn; *** ttDn. &Dn1 = adpolyn; *** 'FIN' Dn1; *** 'REPE' pade1 (N - 1); **Expression de Pn(U(a))******** *** exp = ttDn. (N - 1 - &pade1)/ttDn. (N - 1); *** PnUa=exp*tabamax. &pade1*tadep. &pade1; *** depreel0 = PnUa + depreel0 ; **Expression de lambda(U(a))**** *** Pnlama=exp*tabamax. &pade1*talam. &pade1; *** lamreel0 = Pnlama + lamreel0 ; **Expression de sigma(U(a))**** *** Pnsiga = exp*tabamax. &pade1*tasig. &pade1; *** sigreel0 = Pnsiga + sigreel0; *** 'FIN' pade1; *** 'FINSI'; dlamda = ttlam . 1 - ttlam . (1 - 1) ; dlam = dlamda /dray * (-1 ) ; sgd = signe( dlam) ; * * MESS ' valeur de dlam aprés redirection ' dlam ; * MESS ' LE SIGNE après redirection ' SGd ; * ttdep . 1 = depreel0 ; ttsig . 1 = sigreel0 PICA ( depreel0 - U0 ) mo ; tsetse=ttsig . 1 ; ttlam . 1 = lamreel0; * U0= depreel0 ; 'FINSI'; sigreel0=sigreel0 'PICA' deltau0 mo ; deltatot=deltau0 + deltatot ; sigpiok= ('SIGMA' deltatot mo ma) + S0; sigcau= sigpiok 'PICA' deltatot mo; sigreel0= sigcau; ttdep.0 = depreel0; ttsig.0 = S0; ttlam.0= lamreel0; 'SI' ('NON' ifini) ; *mess ' lamrell0 ' lamreel0; ff= ff *( 1.d0 - lamreel0); depreel0= deltatot 'ENLE' 'LX'; lamreel0=0.D0; 'FORM' geoini; 'FORM' deltatot; U0= deltatot; 'FINSI'; *a0 = tab . a . (1 + na) ; dlamda = (ttlam . 1 - ttlam . 0 ); dlamb = dlamda /dray ; sgdlamda = signe( dlamb) ; *tab.'SIGNE'=sgdlamda; 'FORM' geoini; 'FORM' deltatot; 'SI' ifini; 'MESS' ' sortie de la MAN nombre de segments ' &it IOUT; form geoini; 'QUITTER' IT; 'FINSI'; 'FIN' IT; 'FINPROC' deltatot IOUT;