Télécharger rayle2.eso

Retour à la liste

Numérotation des lignes :

  1. C RAYLE2 SOURCE CB215821 17/07/21 21:15:30 9513
  2. C RAYLEI SOURCE AM 08/02/14 21:28:07 6050
  3. SUBROUTINE RAYLE2 (IPA,IPB,IV,IW,XRVP,XIVP,IPLMOX,IPLMOY,DELTA
  4. & ,IPRX, IPIX, ICVG)
  5. IMPLICIT INTEGER(I-N)
  6. IMPLICIT REAL*8 (A-H, O-Z)
  7.  
  8. *****************************************************************************
  9. * *
  10. * RAYLEI *
  11. * ------ *
  12. *****************************************************************************
  13. *
  14. * Fonction:
  15. * Trouve les valeurs propres d'une matrice de Rayleigh.
  16. *
  17. * création : Benoit Prabel mars 2009
  18. *
  19. *
  20. *****************************************************************************
  21. *
  22. * PROCEDURE D APPEL: CALL RAYLEI (IPA, IPB, IV, IW, XRVP, XIVP,IPLMOX,
  23. * IPLMOY,DELTA, IPRX, IPIX)
  24. *
  25. -INC CCOPTIO
  26. -INC SMLMOTS
  27. -INC SMCHPOI
  28. -INC SMRIGID
  29. -INC CCHAMP
  30.  
  31. REAL*8 A01(1:2,1:2),B01(1:2,1:2)
  32. REAL*8 etaR(1:2),etaI(1:2)
  33. REAL*8 XRVP, XIVP, XRES, alpha0, alpha1, alpha2, DELTA
  34. CHARACTER*4 MOTCLE
  35.  
  36. ICVG = 0
  37.  
  38. *-------Projection du probleme sur la base de vecteurs (V, W)
  39. call MUCPRI (IV, IPA, IPAV)
  40. call MUCPRI (IW, IPA, IPAW)
  41. CALL MUCPRI (IV, IPB, IPBV)
  42. CALL MUCPRI (IW, IPB, IPBW)
  43. call XTY1 (IV, IPAV, IPLMOX, IPLMOY, A01(1,1) )
  44. call XTY1 (IV, IPAW, IPLMOX, IPLMOY, A01(1,2) )
  45. call XTY1 (IW, IPAV, IPLMOX, IPLMOY, A01(2,1) )
  46. call XTY1 (IW, IPAW, IPLMOX, IPLMOY, A01(2,2) )
  47. call XTY1 (IV, IPBV, IPLMOX, IPLMOY, B01(1,1) )
  48. call XTY1 (IV, IPBW, IPLMOX, IPLMOY, B01(1,2) )
  49. call XTY1 (IW, IPBV, IPLMOX, IPLMOY, B01(2,1) )
  50. call XTY1 (IW, IPBW, IPLMOX, IPLMOY, B01(2,2) )
  51. * write(*,*) A01(1,1),A01(1,2),'| -lambda |',B01(1,1),B01(1,2)
  52. * write(*,*) A01(2,1),A01(2,2),'| |',B01(2,1),B01(2,2)
  53.  
  54.  
  55. *-------%il faut annuler le determinant de la matrice |A01 - lambda B01|
  56. * % det(A01 - lambda B01) = alpha0 + alpha1 lambda + alpha2 lambda^2
  57. alpha0 = (A01(1,1)*A01(2,2)) - (A01(2,1)*A01(1,2))
  58. alpha1 = (A01(2,1)*B01(1,2)) + (B01(2,1)*A01(1,2))
  59. & - (A01(1,1)*B01(2,2)) - (B01(1,1)*A01(2,2))
  60. alpha2 = (B01(1,1)*B01(2,2)) - (B01(2,1)*B01(1,2))
  61. DELTA = (alpha1**2) - (4.*alpha0*alpha2)
  62. * write(*,*) 'Rayleigh trouve Delta=',DELTA
  63.  
  64. *-------% selon DELTA, on en deduit le type de val propres:
  65. * %- 1 val propre Reelle double
  66. if( (abs(DELTA)) .lt. 1.D-10) then
  67. goto 200
  68. else
  69. * %- 2 val propres Reelles simples => il faut continuer a iterer
  70. if( DELTA .gt. 0.) then
  71. goto 300
  72. * %- 2 val propres Complexes conjuguées => on resout
  73. elseif( DELTA .lt. 0.) then
  74. goto 100
  75. else
  76. * %normalement, on n'arrive pas la...
  77. return
  78. endif
  79. endif
  80.  
  81. 100 CONTINUE
  82. *-------%on se place dans le cas de 2 vp Complexes conjuguées XRVP +/- i*XIVP
  83. ICVG=1
  84. XRVP = -1.*alpha1 / (2.*alpha2)
  85. XIVP = ((abs(DELTA))**0.5) / (2.*alpha2)
  86. * %reste a calculer eta (=vecteur solution du pb projeté),
  87. * %puis les vect propres associés
  88. * % on fixe la 1ere valeur à 1 + 0*i
  89. etaR(1)=1.
  90. etaI(1)=0.
  91. * % on crée
  92. a11R = A01(1,1) - XRVP*B01(1,1)
  93. a11I = - XIVP*B01(1,1)
  94. a12R = A01(1,2) - XRVP*B01(1,2)
  95. a12I = - XIVP*B01(1,2)
  96. xdenR = (a12R**2) + (a12I**2)
  97. etaR(2) = -1.*(a11R*a12R + a11I*a12I) / xdenR
  98. etaI(2) = -1.*(a11I*a12R - a11R*a12I) / xdenR
  99. * % calcul du vect propre
  100. call ADCHPO (IV, IW, IPRX, etaR(1), etaR(2))
  101. call ADCHPO (IV, IW, IPIX, etaI(1), etaI(2))
  102. * % normalisation Complexe => nouvelle routine NORM1C
  103. * appel a MOTS1 copié de itinv1.eso pour exclure LX
  104. * CALL MOTS1 (IPLMOT,MOTCLE)
  105. CALL MOTS2 (IPLMOT,MOTCLE)
  106. CALL NORM1C (IPRX,IPIX,IPLMOT,MOTCLE,IPRX,IPIX)
  107. * % c'est fini !
  108. goto 900
  109.  
  110. 200 CONTINUE
  111. *-------%on se place dans le cas de 1 val propre Reelle double XRVP
  112. ICVG=2
  113. XRVP = -1.*alpha1 / (2.*alpha2)
  114. XIVP = 0.D0
  115. * % on met les 2 vecteurs orthonormée Réels dans IPRX et IPIX !!!
  116. * % => il faudra donc bien testé DELTA dans le programme
  117. * % qui appelle Rayleigh pour différencier du cas complexe !!!
  118. * CALL MOTS1 (IPLMOT,MOTCLE)
  119. CALL MOTS2 (IPLMOT,MOTCLE)
  120. call NORMA1 (IV,IPLMOT,MOTCLE,IPRX)
  121. call NORMA1 (IW,IPLMOT,MOTCLE,IPIX)
  122. goto 900
  123. *
  124. 300 CONTINUE
  125. *-------%on se place dans le cas de 2 val propres Reelles simples distinctes
  126. XRVP = ((DELTA**0.5) - alpha1) / (2.*alpha2)
  127. XRVP2 = -1.*( (DELTA**0.5) + alpha1 ) / (2.*alpha2)
  128. * write(*,*) XRVP,XRVP2
  129. if( abs(1. - abs(XRVP/XRVP2)) .gt. 1.D-5) then
  130. * % => il faut donc continuer a iterer pour les séparer (ou les confondre) + nettement
  131. ICVG=0
  132. XRVP = 0.D0
  133. XIVP = 0.D0
  134. IPRX = 0
  135. IPIX = 0
  136. else
  137. * %on a 2 val p reelle simple de meme module -> on resout
  138. * -> non! on prefere n'en donner qu'une...
  139. ICVG=3
  140. * XIVP = XRVP2
  141. if( (XRVP .lt. 0.) .and. (XRVP2 .gt. 0.)) XRVP=XRVP2
  142. XIVP = 0.D0
  143. * % on fixe la 1ere valeur du 1er vecteur propre à 1 + 0*i
  144. etaR(1) = 1.
  145. a11R = A01(1,1) - XRVP*B01(1,1)
  146. a12R = A01(1,2) - XRVP*B01(1,2)
  147. if( abs(a12R) .lt. 1.E-10 ) then
  148. etaR(1) = 0.
  149. etaR(2) = 1.
  150. else
  151. etaR(2) = -1.*(a11R/a12R)
  152. endif
  153. call ADCHPO (IV, IW, IPRX, etaR(1), etaR(2))
  154. * CALL MOTS1 (IPLMOT,MOTCLE)
  155. CALL MOTS2 (IPLMOT,MOTCLE)
  156. CALL NORMA1 (IPRX,IPLMOT,MOTCLE,IPRX)
  157. * % on fixe la 1ere valeur du 2eme vecteur propre à 1 + 0*i
  158. * etaR(1) = 1.
  159. * a11R = A01(1,1) - XRVP2*B01(1,1)
  160. * a12R = A01(1,2) - XRVP2*B01(1,2)
  161. * if( abs(a12R) .lt. 1.E-10 ) then
  162. * etaR(1) = 0.
  163. * etaR(2) = 1.
  164. * else
  165. * etaR(2) = -1.*(a11R/a12R)
  166. * endif
  167. * call ADCHPO (IV, IW, IPRX2, etaR(1), etaR(2))
  168. * CALL MOTS1 (IPLMOT,MOTCLE)
  169. * CALL NORMA1 (IPRX2,IPLMOT,MOTCLE,IPRX2)
  170. * IPIX = IPRX2
  171. IPIX = 0
  172. endif
  173.  
  174. 900 CONTINUE
  175.  
  176. END
  177.  
  178.  
  179.  
  180.  
  181.  
  182.  

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