Télécharger rayle2.eso

Retour à la liste

Numérotation des lignes :

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

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