C RAYLE2 SOURCE CHAT 09/10/09 21:22:21 6519 C RAYLEI SOURCE AM 08/02/14 21:28:07 6050 SUBROUTINE RAYLE2 (IPA,IPB,IV,IW,XRVP,XIVP,IPLMOX,IPLMOY,DELTA & ,IPRX, IPIX, ICVG) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H, O-Z) ***************************************************************************** * * * RAYLEI * * ------ * ***************************************************************************** * * Fonction: * Trouve les valeurs propres d'une matrice de Rayleigh. * * création : Benoit Prabel mars 2009 * * ***************************************************************************** * * PROCEDURE D APPEL: CALL RAYLEI (IPA, IPB, IV, IW, XRVP, XIVP,IPLMOX, * IPLMOY,DELTA, IPRX, IPIX) * -INC CCOPTIO -INC SMLMOTS -INC SMCHPOI -INC SMRIGID -INC CCHAMP REAL*8 A01(1:2,1:2),B01(1:2,1:2) REAL*8 etaR(1:2),etaI(1:2) REAL*8 XRVP, XIVP, XRES, alpha0, alpha1, alpha2, DELTA ICVG = 0 *-------Projection du probleme sur la base de vecteurs (V, W) call MUCPRI (IV, IPA, IPAV) call MUCPRI (IW, IPA, IPAW) CALL MUCPRI (IV, IPB, IPBV) CALL MUCPRI (IW, IPB, IPBW) call XTY1 (IV, IPAV, IPLMOX, IPLMOY, A01(1,1) ) call XTY1 (IV, IPAW, IPLMOX, IPLMOY, A01(1,2) ) call XTY1 (IW, IPAV, IPLMOX, IPLMOY, A01(2,1) ) call XTY1 (IW, IPAW, IPLMOX, IPLMOY, A01(2,2) ) call XTY1 (IV, IPBV, IPLMOX, IPLMOY, B01(1,1) ) call XTY1 (IV, IPBW, IPLMOX, IPLMOY, B01(1,2) ) call XTY1 (IW, IPBV, IPLMOX, IPLMOY, B01(2,1) ) call XTY1 (IW, IPBW, IPLMOX, IPLMOY, B01(2,2) ) * write(*,*) A01(1,1),A01(1,2),'| -lambda |',B01(1,1),B01(1,2) * write(*,*) A01(2,1),A01(2,2),'| |',B01(2,1),B01(2,2) *-------%il faut annuler le determinant de la matrice |A01 - lambda B01| * % det(A01 - lambda B01) = alpha0 + alpha1 lambda + alpha2 lambda^2 alpha0 = (A01(1,1)*A01(2,2)) - (A01(2,1)*A01(1,2)) alpha1 = (A01(2,1)*B01(1,2)) + (B01(2,1)*A01(1,2)) & - (A01(1,1)*B01(2,2)) - (B01(1,1)*A01(2,2)) alpha2 = (B01(1,1)*B01(2,2)) - (B01(2,1)*B01(1,2)) DELTA = (alpha1**2) - (4.*alpha0*alpha2) * write(*,*) 'Rayleigh trouve Delta=',DELTA *-------% selon DELTA, on en deduit le type de val propres: * %- 1 val propre Reelle double if( (abs(DELTA)) .lt. 1.D-10) then goto 200 else * %- 2 val propres Reelles simples => il faut continuer a iterer if( DELTA .gt. 0.) then goto 300 * %- 2 val propres Complexes conjuguées => on resout elseif( DELTA .lt. 0.) then goto 100 else * %normalement, on n'arrive pas la... return endif endif 100 CONTINUE *-------%on se place dans le cas de 2 vp Complexes conjuguées XRVP +/- i*XIVP ICVG=1 XRVP = -1.*alpha1 / (2.*alpha2) XIVP = ((abs(DELTA))**0.5) / (2.*alpha2) * %reste a calculer eta (=vecteur solution du pb projeté), * %puis les vect propres associés * % on fixe la 1ere valeur à 1 + 0*i etaR(1)=1. etaI(1)=0. * % on crée a11R = A01(1,1) - XRVP*B01(1,1) a11I = - XIVP*B01(1,1) a12R = A01(1,2) - XRVP*B01(1,2) a12I = - XIVP*B01(1,2) xdenR = (a12R**2) + (a12I**2) etaR(2) = -1.*(a11R*a12R + a11I*a12I) / xdenR etaI(2) = -1.*(a11I*a12R - a11R*a12I) / xdenR * % calcul du vect propre call ADCHPO (IV, IW, IPRX, etaR(1), etaR(2)) call ADCHPO (IV, IW, IPIX, etaI(1), etaI(2)) * % normalisation Complexe => nouvelle routine NORM1C * appel a MOTS1 copié de itinv1.eso pour exclure LX * CALL MOTS1 (IPLMOT,MOTCLE) CALL MOTS2 (IPLMOT,MOTCLE) CALL NORM1C (IPRX,IPIX,IPLMOT,MOTCLE,IPRX,IPIX) * % c'est fini ! goto 900 200 CONTINUE *-------%on se place dans le cas de 1 val propre Reelle double XRVP ICVG=2 XRVP = -1.*alpha1 / (2.*alpha2) XIVP = 0.D0 * % on met les 2 vecteurs orthonormée Réels dans IPRX et IPIX !!! * % => il faudra donc bien testé DELTA dans le programme * % qui appelle Rayleigh pour différencier du cas complexe !!! * CALL MOTS1 (IPLMOT,MOTCLE) CALL MOTS2 (IPLMOT,MOTCLE) call NORMA1 (IV,IPLMOT,MOTCLE,IPRX) call NORMA1 (IW,IPLMOT,MOTCLE,IPIX) goto 900 * 300 CONTINUE *-------%on se place dans le cas de 2 val propres Reelles simples distinctes XRVP = ((DELTA**0.5) - alpha1) / (2.*alpha2) XRVP2 = -1.*( (DELTA**0.5) + alpha1 ) / (2.*alpha2) * write(*,*) XRVP,XRVP2 if( abs(1. - abs(XRVP/XRVP2)) .gt. 1.D-5) then * % => il faut donc continuer a iterer pour les séparer (ou les confondre) + nettement ICVG=0 XRVP = 0.D0 XIVP = 0.D0 IPRX = 0 IPIX = 0 else * %on a 2 val p reelle simple de meme module -> on resout * -> non! on prefere n'en donner qu'une... ICVG=3 * XIVP = XRVP2 if( (XRVP .lt. 0.) .and. (XRVP2 .gt. 0.)) XRVP=XRVP2 XIVP = 0.D0 * % on fixe la 1ere valeur du 1er vecteur propre à 1 + 0*i etaR(1) = 1. a11R = A01(1,1) - XRVP*B01(1,1) a12R = A01(1,2) - XRVP*B01(1,2) if( abs(a12R) .lt. 1.E-10 ) then etaR(1) = 0. etaR(2) = 1. else etaR(2) = -1.*(a11R/a12R) endif call ADCHPO (IV, IW, IPRX, etaR(1), etaR(2)) * CALL MOTS1 (IPLMOT,MOTCLE) CALL MOTS2 (IPLMOT,MOTCLE) CALL NORMA1 (IPRX,IPLMOT,MOTCLE,IPRX) * % on fixe la 1ere valeur du 2eme vecteur propre à 1 + 0*i * etaR(1) = 1. * a11R = A01(1,1) - XRVP2*B01(1,1) * a12R = A01(1,2) - XRVP2*B01(1,2) * if( abs(a12R) .lt. 1.E-10 ) then * etaR(1) = 0. * etaR(2) = 1. * else * etaR(2) = -1.*(a11R/a12R) * endif * call ADCHPO (IV, IW, IPRX2, etaR(1), etaR(2)) * CALL MOTS1 (IPLMOT,MOTCLE) * CALL NORMA1 (IPRX2,IPLMOT,MOTCLE,IPRX2) * IPIX = IPRX2 IPIX = 0 endif 900 CONTINUE END