C HBMBIF SOURCE OF166741 26/05/11 21:15:07 12538 *======================================================================= * Calcul des points de bifurcation par une methode de Newton-Raphson *======================================================================= * TYPC indique le type de calcul a realiser * = 'L' pour un point limite * = 'B' pour un point de branchement * = 'P' pour un doublement de periode * = 'N' pour une bifurcation de Neimark-Sacker *======================================================================= SUBROUTINE HBMBIF(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,MTPHI,KTEMP,MTLIAA, & MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,TYPC, & KPARNUM,KSORT) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMLREEL POINTEUR fvec.MLREEL, fvec0.MLREEL -INC TMDYNC * segment local (a nettoyer) SEGMENT MWORK INTEGER IPIV33(2*NT+2),INDEIG(2*NT),I4(3*NT+2),I5(2*NT+1) REAL*8 AVVP(NT),AVVP2(NT+1),V1(NT) REAL*8 EIGRE(NT),EIGIM(NT),QIN(NT),OMEGIN,KAPPA,KAPPA2 REAL*8 EIGRE2(NT+1),EIGIM2(NT+1) REAL*8 VR(NT,NT),VL(NT,NT),WORK(8*NT) REAL*8 VR2(NT+1,NT+1),VL2(NT+1,NT+1),WORK2(8*(NT+1)) REAL*8 JAC1(NT,NT),JAC2(NT,NT) REAL*8 RxwPhi(NT),RXXPHI(NT,NT),RwwPhi REAL*8 vA1(NT), vA2(NT),vA3(NT),vB1(NT),vB2(NT),vB3(NT) REAL*8 RX2(NT,NT),SS,EI(NT) REAL*8 JBP(2*NT+2,2*NT+2),solBP(2*NT+2) REAL*8 dWB, dXB(NT),dPHIR(NT),dPHII(NT),WORK2w(16*NT) REAL*8 res2(NT),Rco(NT,NT),PHIR(NT),PHII(NT),DPsiPsi(NT,NT) REAL*8 MREG(NT+1,NT+1),DEL1(NT,NT),DEL2(NT),Mi,Ci,JH(2*NT,2*NT) REAL*8 EXRE2w(2*NT),EXIM2w(2*NT),VL2w(2*NT,2*NT),VR2w(2*NT,4*NT) REAL*8 dRxDdk1(NT),dRxDdk2(NT) REAL*8 Rxmk2D(NT,NT),D10(NT,NT),D1w(NT,NT),D2(NT,NT),VQ(NT) REAL*8 RxxPHIR(NT,NT),RxxPHII(NT,NT),RxwPhiR(NT),RxwPhiI(NT) REAL*8 JNS(3*NT+2,3*NT+2),solNS(3*NT+2),resV1(NT),resV2(NT) REAL*8 JPD(3*NT+1,3*NT+1),solPD(3*NT+1),D1wP(NT,NT) ENDSEGMENT CHARACTER*1 TYPC LOGICAL CHECK,CANAL INTEGER MAXITS,MAXVCP,INFO REAL*8 ZERO,ONE,TWO,XEPS PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0, XEPS=1.D-6) C Tableaux locaux : REAL*8 AA(3,3),sol(3),arrRes(3),arrRes2(4),arrRes3(5) C Fonctions BLAS/LAPACK REAL*8 DDOT, DNRM2 EXTERNAL DDOT, DNRM2 *----------------------------------------------------------------------- PARAMETER(MAXITS=40) * TOLF=1.e-4 * TOLMIN=1.e-6 * STPMX=100. * MTQ = KTQ MTKAM = KTKAM MTEMP = KTEMP PSORT = KSORT PARNUM = KPARNUM CANAL = JANAL SEGINI,MWORK JG = NT SEGINI,fvec,fvec0 nmc = NT SEGINI,MdZw c* DO J = 1, nmd c* DO I = 1, nmd c* MdZW.MATRC(I,J)=0.D0 c* ENDDO c* ENDDO * sauvegarde du point de la courbe de reponse OMEGIN = OMEG DO I = 1,NT QIN(I) = Q1(I) ENDDO * CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec0,fref0, & .false.) CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.) c ---------------------------------------------------------------- c !!! Localisation NeimarkSacker et PeriodDoubling a debugger !!! c if(TYPC.eq.'N'.or.TYPC.eq.'P') then c WRITE(IOIMP,*) 'Localisation bifurcation ',TYPC,'non prevue !' c CHECK = .true. c return c endif c ---------------------------------------------------------------- c message IF (IIMPI.GE.2) WRITE(IOIMP,*) 'Localisation bifurcation ',TYPC IF (IIMPI.GE.3) WRITE(IOIMP,990) 0,dnrm2(NT,fvec0.prog,1),fref0 990 FORMAT('+ hbmbif : iter ',I2,' -> |R|=',E10.3,' ref=',E10.3) *======================================================================= * Initialisation *======================================================================= * La Jacobienne augmentee a la structure de base suivante: * JJ = [ Rx O Rw ] * [ (RxPhi)x Rx RxwPhi] * [ 0' 2Phi' 0 ] * qui est a adapter selon la bifurcation traitee. * * Il faut initialiser le vecteur propre Phi. * Calcul de la Jacobienne Rx = Z(OMEGA) - dFNLdX DO J=1,NT DO I=1,NT RX(I,J) = ZZ(I,J)-JAC(I,J) ENDDO ENDDO CALL COPYMAT(NT,Rx,Rco) *----------------------------------------------------------------------- * Bifurcation statique: le vecteur propre annule la Jacobienne *----------------------------------------------------------------------- IF (TYPC.EQ.'L') THEN CALL DGEEV('N','V',NT,Rco,NT,EIGRE,EIGIM,VL,1,VR,NT, & WORK,8*NT,INFO) * On prend le vecteur associe a la partie reelle la plus petite. MINVP=1 XMINVP=ABS(EIGRE(1)) DO I = 1,NT AVVP(I) = ABS(EIGRE(I)) IF(AVVP(I).LT.XMINVP) THEN MINVP=I XMINVP=AVVP(I) ENDIF ENDDO c MINVP = MINLOC(AVVP,NT) * Les vecteurs propres sont purement reels dans ce cas. DO I = 1,NT PHIR(I) = VR(I,MINVP) PHII(I) = ZERO ENDDO * Pas de nouvelle frequence pour les cycles bifurques KAPPA = ZERO MAXVCP=1 XMAXVCP=ABS(PHIR(1)) DO I = 1,NT AVVP(I) = ABS(PHIR(I)) IF(AVVP(I).GT.XMAXVCP) THEN MAXVCP=I XMAXVCP=AVVP(I) ENDIF ENDDO c MAXVCP = MAXLOC(AVVP,NT) * Elements pour la resolution par blocs SS = ZERO DO I=1,NT SS = SS + abs(RX(I,I)) ENDDO SS = SS/NT * Penalisation de la Jacobienne (evite mauvais conditionnement) CALL COPYMAT(NT,RX,RX2) DO I=1,NT EI(I) = ZERO ENDDO RX2(MAXVCP,MAXVCP) = RX2(MAXVCP,MAXVCP) + SS EI(MAXVCP) = ONE ENDIF *----------------------------------------------------------------------- * Branch point ? *----------------------------------------------------------------------- IF (TYPC.EQ.'B') THEN CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.) IF (BAL.EQ.1) THEN DO I=1,NT Rw(I) = Rw(I) - TWO*OMEG*FEXA(I) ENDDO END IF * Regularisation de la Jacobienne CALL HBMDVEC(NT,NHBM,NDDL,Q1,ONE,V1) DO I = 1,NT MREG(I,NT+1) = Rw(I) MREG(NT+1,I) = Rw(I) ENDDO znormV1 = ONE / dnrm2(NT,V1,1) DO I = 1,NT V1(I) = V1(I)*znormV1 ENDDO DO J = 1,NT DO I = 1,NT RX2(I,J) = RX(I,J) + V1(I)*V1(J) MREG(I,J) = RX2(I,J) ENDDO ENDDO MREG(NT+1,NT+1) = ZERO CALL DGEEV('N','V',NT+1,MREG,NT+1,EIGRE2,EIGIM2,VL2,1,VR2,NT+1, & WORK2,8*(NT+1),INFO) * On prend le vecteur associe a partie reelle la plus petite. MINVP=1 XMINVP=ABS(EIGRE2(1)) DO I = 1,NT+1 AVVP2(I) = ABS(EIGRE2(I)) IF(AVVP2(I).LT.XMINVP) THEN MINVP=I XMINVP=AVVP2(I) ENDIF ENDDO c MINVP = MINLOC(AVVP2,NT+1) * Les vecteurs propres sont purement reels dans ce cas. DO I = 1,NT PHIR(I) = VR2(I,MINVP) PHII(I) = ZERO ENDDO * Pas de nouvelle frequence pour les cycles bifurques KAPPA = ZERO ENDIF *----------------------------------------------------------------------- * Bifurcations dynamiques (Period doubling ou Neimark Sacker) *----------------------------------------------------------------------- * Construction de la matrice de Hill, JJ IF ((TYPC.EQ.'P').OR.(TYPC.EQ.'N')) THEN DO I = 1,NT DO J = 1,NT D10(I,J)=ZERO D1w(I,J)=ZERO D2(I,J) =ZERO ENDDO ENDDO DO I=1,2*NHBM+1 DO J=1,NDDL DEL2(NDDL*(I-1)+J) = ONE/XM(J,1) D2(NDDL*(I-1)+J,NDDL*(I-1)+J) = XM(J,1) ENDDO ENDDO DO I=1,NDDL D10(I,I) = XASM(I,1) ENDDO DO J = 2,2*NHBM,2 DO I=1,NDDL D10(NDDL*(1+(J-2))+I,NDDL*(1+(J-2))+I) = XASM(I,1) D1w(NDDL*(1+(J-2))+I,NDDL*(1+(J-1))+I) = J*XM(I,1) D1w(NDDL*(1+(J-1))+I,NDDL*(1+(J-2))+I) = -J*XM(I,1) D10(NDDL*(1+(J-1))+I,NDDL*(1+(J-1))+I) = XASM(I,1) ENDDO ENDDO * DEL1 = DEL10 + OMEG*DEL1w CALL COLIMA(NT,NT,D10,D1w,ONE,OMEG,DEL1) * JJ = [DEL2\DEL1 DEL2\Rx] * [ -I 0 ] * DO I=1,NT DO J=1,NT * JJ_11 = DEL2\DEL1 JH(I,J) = -DEL1(I,J)*DEL2(J) * JJ_12 = DEL2\Rx JH(I,NT+J) = -RX(I,J)*DEL2(J) * JJ_21 = [I] IF (I.EQ.J) THEN JH(NT+I,J) = ONE ELSE JH(NT+I,J) = ZERO ENDIF * JJ_22 = [0] JH(NT+I,NT+J) = ZERO ENDDO ENDDO CALL DGEEV('N','V',2*NT,JH,2*NT,EXRE2w,EXIM2w,VL2w,1,VR2w,2*NT, & WORK2,16*NT,INFO) CALL HBMORDO(2*NT,NDDL,OMEG,EXIM2w,INDEIG) * WRITE(*,*) ' REAL ', ' IMAG ', ' INDX ' * DO I = 1,2*NT * WRITE(*,*) EXRE2w(INDEIG(I)),EXIM2w(I),INDEIG(I) * ENDDO * WRITE(*,*) 'INDXS:',(INDEIG(IOU),IOU=1,2*NDDL) * * On choisit le vecteur associe a la valeur propre avec partie * imaginaire positive du mode instable MINVP=1 DO I=1,2*NDDL IF ((EXRE2w(INDEIG(I)).GT.ZERO).AND.(EXIM2w(I).GT.ZERO)) THEN MINVP=INDEIG(I) INDK =I * WRITE(*,*) 'Mode instable:',MINVP ENDIF ENDDO DO I=1,NT PHIR(I) = VR2w(I,MINVP) PHII(I) = VR2w(I,MINVP+1) ENDDO * WRITE(*,*) 'PHI_R ', ' PHI_I' * DO IOU = 1,NT * WRITE(*,*) PHIR(IOU),PHII(IOU) * ENDDO * STOP MAXVCP=1 XMAXVCP=ABS(PHII(1)) DO I = 1,NT AVVP(I) = ABS(PHII(I)) IF (AVVP(I).GT.XMAXVCP) THEN MAXVCP=I XMAXVCP=AVVP(I) ENDIF ENDDO c MAXVCP = MAXLOC(AVVP,NT) DO I = 1,NT VQ(I) = ZERO ENDDO VQ(MAXVCP) = ONE * --------------------- IF (TYPC.EQ.'P') THEN * Nouvelle frequence fixee a OMEGA/2 KAPPA = OMEG/TWO ENDIF * --------------------- IF (TYPC.EQ.'N') THEN * Nouvelle frequence incommensurable KAPPA KAPPA = ABS(EXIM2w(INDK)) * write(*,*) 'KAPPA^0=',KAPPA ENDIF * STOP ENDIF *======================================================================= * Iterations *======================================================================= DO its=1,MAXITS *----------------------------------------------------------------------- * CAS L: POINT LIMITE *----------------------------------------------------------------------- IF (TYPC.EQ.'L') THEN *** Calcul des derivees * RxxPhi EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS) * Perturbation de la Jacobienne: +EPS DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC1) * Perturbation de la Jacobienne: -EPS DO JJ = 1,NT Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC2) * retour a la valeur initiale DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ) ENDDO * Differences centrees r_z = ONE / (TWO*EPS1) DO I = 1,NT DO J = 1,NT RxxPHI(I,J) = (JAC2(I,J)-JAC1(I,J))*r_z ENDDO ENDDO * RxwPhi CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,MdZw) DO I=1,NT r_z = ZERO DO J = 1,NT r_z = r_z + MdZw.MATRC(I,J)*PHIR(J) ENDDO RxwPhi(I) = r_z ENDDO * Rw CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.) IF (BAL.EQ.1) THEN r_z = TWO*OMEG DO I=1,NT Rw(I) = Rw(I) - r_z * FEXA(I) ENDDO END IF * Resolution par blocs * A1 DO I = 1,NT vA1(I) = -fvec0.PROG(I) ENDDO CALL COPYMAT(NT,RX2,Rco) CALL DGESV(NT,1,Rco,NT,IPIV,vA1,NT,INFO) * A2 DO I = 1,NT vA2(I) = -Rw(I) ENDDO CALL COPYMAT(NT,RX2,Rco) CALL DGESV(NT,1,Rco,NT,IPIV,vA2,NT,INFO) * A3 DO I = 1,NT vA3(I) = EI(I) ENDDO CALL COPYMAT(NT,RX2,Rco) CALL DGESV(NT,1,Rco,NT,IPIV,vA3,NT,INFO) * B1 DO I = 1,NT r_z = ZERO DO J = 1,NT r_z = r_z + RxxPHI(I,J)*vA1(J) + RX(I,J)*PHIR(J) ENDDO vB1(I) = -r_z ENDDO CALL COPYMAT(NT,RX2,Rco) CALL DGESV(NT,1,Rco,NT,IPIV,vB1,NT,INFO) * B2 DO I = 1,NT r_z = ZERO DO J = 1,NT r_z = r_z + RxxPHI(I,J)*vA2(J) ENDDO vB2(I) = -(r_z + RxwPhi(I)) ENDDO CALL COPYMAT(NT,RX2,Rco) CALL DGESV(NT,1,Rco,NT,IPIV,vB2,NT,INFO) * B3 DO I = 1,NT r_z = ZERO DO J = 1,NT r_z = r_z + RxxPHI(I,J)*vA3(J) ENDDO vB3(I) = -r_z ENDDO CALL COPYMAT(NT,RX2,Rco) CALL DGESV(NT,1,Rco,NT,IPIV,vB3,NT,INFO) * * Construction de la matrice du probleme AA(1,1) = vA2(MAXVCP) AA(1,2) = SS*vA3(MAXVCP)-ONE AA(1,3) = ZERO AA(2,1) = vB2(MAXVCP) AA(2,2) = SS*vB3(MAXVCP) AA(2,3) = SS*vA3(MAXVCP)-ONE AA(3,1) = DDOT(NT,PHIR,1,vB2,1) AA(3,2) = SS*(DDOT(NT,PHIR,1,vB3,1)) AA(3,3) = SS*(DDOT(NT,PHIR,1,vA3,1)) * Construction du cote droit sol(1) = -vA1(MAXVCP) sol(2) = -vB1(MAXVCP) sol(3) = 0.5D0*(ONE-DDOT(NT,PHIR,1,PHIR,1)) & -DDOT(NT,PHIR,1,vB1,1) * Resolution CALL DGESV(3,1,AA,3,IPIV,sol,3,INFO) * Corrections dWB = sol(1) beta1 = sol(2) beta2 = sol(3) DO I = 1,NT dXB(I) = vA1(I) + dWB*vA2(I) +SS*beta1*vA3(I) dPHIR(I)= vB1(I)+ dWB*vB2(I) +SS*beta1*vB3(I) +SS*beta2*vA3(I) Q1(I) = Q1(I) + dXB(I) PHIR(I) = PHIR(I) + dPHIR(I) ENDDO OMEG = OMEG + dWB ENDIF *----------------------------------------------------------------------- * CAS B: POINT DE BRANCHEMENT *----------------------------------------------------------------------- IF (TYPC.EQ.'B') THEN *** Calcul des derivees * RxxPhi EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS) * Perturbation de la Jacobienne: +EPS DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC1) * Perturbation de la Jacobienne: -EPS DO JJ = 1,NT Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC2) DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ) ENDDO c DPsiPsi CALL HBMDPP(NT,NDDL,Q1,PHIR,DPsiPsi) * Differences centrees DO I = 1,NT DO J = 1,NT RxxPHI(I,J)=((JAC2(I,J)-JAC1(I,J))/(TWO*EPS1)) + DPsiPsi(I,J) ENDDO ENDDO * RxwPhi CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,MdZw) DO I=1,NT r_z = ZERO DO J = 1,NT r_z = r_z + MdZw.MATRC(I,J)*PHIR(J) ENDDO RxwPhi(I) = r_z ENDDO * Rw CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.) IF (BAL.EQ.1) THEN DO I=1,NT Rw(I) = Rw(I) - TWO*OMEG*FEXA(I) ENDDO END IF * Rww*Phi CALL HBMRWW(NT,NDDL,Q1,PHIR,XM,RwwPhi) IF (BAL.EQ.1) THEN DO I=1,NT Rw(I) = Rw(I) - TWO*OMEG*FEXA(I) RwwPhi = RwwPhi - TWO*FEXA(I)*PHIR(I) ENDDO END IF * Construction de la matrice du probleme et du cote droit solBP(2*NT+1) = ZERO DO I = 1,NT solBP(NT+I) = ZERO DO J = 1,NT JBP(I,J) = RX(I,J) JBP(NT+I,J) = RxxPHI(I,J) JBP(I,NT+J) = ZERO JBP(NT+I,NT+J) = RX2(I,J) solBP(NT+I) = solBP(NT+I) + RX2(I,J)*PHIR(J) ENDDO JBP(2*NT+1,NT+I) = RxwPhi(I) JBP(2*NT+1,I) = Rw(I) JBP(2*NT+2,I) = ZERO JBP(2*NT+2,NT+I) = TWO*PHIR(I) JBP(I,2*NT+1) = Rw(I) JBP(NT+I,2*NT+1) = RxwPhi(I) JBP(I,2*NT+2) = PHIR(I) JBP(NT+I,2*NT+2) = ZERO solBP(I) = -fvec0.PROG(I) solBP(NT+I) = -solBP(NT+I) ENDDO solBP(2*NT+1) = -(DDOT(NT,Rw,1,PHIR,1)) solBP(2*NT+2) = -(DDOT(NT,PHIR,1,PHIR,1)-ONE) JBP(2*NT+1,2*NT+1) = RwwPhi JBP(2*NT+1,2*NT+2) = ZERO JBP(2*NT+2,2*NT+1) = ZERO JBP(2*NT+2,2*NT+2) = ZERO * * Resolution CALL DGESV(2*NT+2,1,JBP,2*NT+2,IPIV33,solBP,2*NT+2,INFO) * Corrections dWB = solBP(2*NT+1) DO I = 1,NT dXB(I) = solBP(I) dPHIR(I)= solBP(NT+I) Q1(I) = Q1(I) + dXB(I) PHIR(I) = PHIR(I) + dPHIR(I) ENDDO OMEG = OMEG + dWB ENDIF *----------------------------------------------------------------------- * CAS P: DOUBLEMENT DE PERIODE *----------------------------------------------------------------------- IF (TYPC.EQ.'P') THEN *** Calcul des derivees * RxxPhiR EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS) * Perturbation de la Jacobienne: +EPS1 DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC1) * Perturbation de la Jacobienne: -EPS1 DO JJ = 1,NT Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC2) DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ) ENDDO * Differences centrees DO I = 1,NT DO J = 1,NT RxxPHIR(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS1) ENDDO ENDDO * RxxPhiI EPS2 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHII,1) + XEPS) * Perturbation de la Jacobienne: +EPS2 DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC1) * Perturbation de la Jacobienne: -EPS2 DO JJ = 1,NT Q1(JJ) = Q1(JJ) - TWO*EPS2*PHII(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC2) DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ) ENDDO * Differences centrees DO I = 1,NT DO J = 1,NT RxxPHII(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS2) ENDDO ENDDO * Rw CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.) IF (BAL.EQ.1) THEN DO I=1,NT Rw(I) = Rw(I) - TWO*OMEG*FEXA(I) ENDDO END IF * Calcul de: (Rxw-(OMEG/2)*D2)*PhiR-(0.5*D10+OMEG*D1w)*PhiI et * (Rxw-(OMEG/2)*D2)*PhiI+(0.5*D10+OMEG*D1w)*PhiR CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,MdZw) CALL COLIMA(NT,NT,DEL1,D10,ONE,-0.5D0,D1wP) DO I=1,NT RxwPhiR(I) = ZERO RxwPhiI(I) = ZERO DO J = 1,NT RxwPhiR(I)=RxwPhiR(I)+(MdZw.MATRC(I,J)-KAPPA*D2(I,J))*PHIR(J) & - D1wP(I,J)*PHII(J) RxwPhiI(I)=RxwPhiI(I)+(MdZw.MATRC(I,J)-KAPPA*D2(I,J))*PHII(J) & + D1wP(I,J)*PHIR(J) ENDDO ENDDO * Rxmk2D = Rx-(KAPPA**2)*D2 CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D) * Construction de la matrice du probleme et du cote droit DO I = 1,NT solPD(NT+I) = ZERO solPD(2*NT+I) = ZERO DO J = 1,NT JPD(I,J) = RX(I,J) JPD(NT+I,J) = RxxPHIR(I,J) JPD(2*NT+I,J) = RxxPHII(I,J) JPD(I,NT+J) = ZERO JPD(I,2*NT+J) = ZERO JPD(NT+I,NT+J) = Rxmk2D(I,J) JPD(NT+I,2*NT+J) = -KAPPA*DEL1(I,J) JPD(2*NT+I,NT+J) = KAPPA*DEL1(I,J) JPD(2*NT+I,2*NT+J) = Rxmk2D(I,J) solPD(NT+I)=solPD(NT+I)+Rxmk2D(I,J)*PHIR(J) & -KAPPA*DEL1(I,J)*PHII(J) solPD(2*NT+I)=solPD(2*NT+I)+Rxmk2D(I,J)*PHII(J) & +KAPPA*DEL1(I,J)*PHIR(J) ENDDO JPD( I,3*NT+1) = Rw(I) JPD( NT+I,3*NT+1) = RxwPhiR(I) JPD(2*NT+I,3*NT+1) = RxwPhiI(I) JPD(3*NT+1,I) = ZERO JPD(3*NT+1,2*NT+1) = VQ(I) JPD(3*NT+1,NT+1) = ZERO solPD(I) = -fvec0.PROG(I) solPD(NT+I) = -solNS(NT+I) solPD(2*NT+I) = -solNS(2*NT+I) ENDDO solPD(3*NT+1) = -(DDOT(NT,VQ,1,PHII,1)-0.5D0) JNS(3*NT+1,3*NT+1) = ZERO * Resolution CALL DGESV(3*NT+1,1,JPD,3*NT+1,I5,solPD,3*NT+1,INFO) * Corrections dWB = solPD(3*NT+1) DO I = 1,NT dXB(I) = solPD(I) dPHIR(I)= solPD(NT+I) dPHII(I)= solPD(2*NT+I) Q1(I) = Q1(I) + dXB(I) PHIR(I) = PHIR(I) + dPHIR(I) PHII(I) = PHII(I) + dPHII(I) ENDDO OMEG = OMEG + dWB KAPPA= OMEG/TWO END IF *----------------------------------------------------------------------- * CAS N: BIFURCATION DE NEIMARK-SACKER *----------------------------------------------------------------------- IF (TYPC.EQ.'N') THEN *** Calcul des derivees * RxxPhiR EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS) * Perturbation de la Jacobienne: +EPS1 DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC1) * Perturbation de la Jacobienne: -EPS1 DO JJ = 1,NT Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC2) DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ) ENDDO * Differences centrees DO I = 1,NT DO J = 1,NT RxxPHIR(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS1) ENDDO ENDDO * RxxPhiI EPS2 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHII,1) + XEPS) * Perturbation de la Jacobienne: +EPS2 DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC1) * Perturbation de la Jacobienne: -EPS2 DO JJ = 1,NT Q1(JJ) = Q1(JJ) - TWO*EPS2*PHII(JJ) ENDDO CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref, & .false.) CALL COPYMAT(NT,JAC,JAC2) DO JJ = 1,NT Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ) ENDDO * Differences centrees DO I = 1,NT DO J = 1,NT RxxPHII(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS2) ENDDO ENDDO * RxwPhiR-KAPPA*D1w*PhiI, RxwPhiI+KAPPA*D1w*PhiR CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,MdZw) DO I=1,NT RxwPhiR(I) = ZERO RxwPhiI(I) = ZERO dRxDdk1(I) = ZERO dRxDdk2(I) = ZERO DO J = 1,NT RxwPhiR(I)=RxwPhiR(I)+MdZw.MATRC(I,J)*PHIR(J) & -KAPPA*D1w(I,J)*PHII(J) RxwPhiI(I)=RxwPhiI(I)+MdZw.MATRC(I,J)*PHII(J) & +KAPPA*D1w(I,J)*PHIR(J) dRxDdk1(I)=dRxDdk1(I)-TWO*KAPPA*D2(I,J)*PHIR(J) & -DEL1(I,J)*PHII(J) dRxDdk2(I)=dRxDdk2(I)-TWO*KAPPA*D2(I,J)*PHII(J) & +DEL1(I,J)*PHIR(J) ENDDO ENDDO * Rxmk2D = Rx-(KAPPA**2)*D2 CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D) * Rw CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.) IF (BAL.EQ.1) THEN DO I=1,NT Rw(I) = Rw(I) - TWO*OMEG*FEXA(I) ENDDO END IF * Construction de la matrice du probleme et du cote droit DO I = 1,NT solNS(NT+I) = ZERO solNS(2*NT+I) = ZERO DO J = 1,NT JNS(I,J) = RX(I,J) JNS(NT+I,J) = RxxPHIR(I,J) JNS(2*NT+I,J) = RxxPHII(I,J) JNS(I,NT+J) = ZERO JNS(I,2*NT+J) = ZERO JNS(NT+I,NT+J) = Rxmk2D(I,J) JNS(NT+I,2*NT+J) = -KAPPA*DEL1(I,J) JNS(2*NT+I,NT+J) = KAPPA*DEL1(I,J) JNS(2*NT+I,2*NT+J) = Rxmk2D(I,J) solNS(NT+I)=solNS(NT+I)+Rxmk2D(I,J)*PHIR(J) & -KAPPA*DEL1(I,J)*PHII(J) solNS(2*NT+I)=solNS(2*NT+I)+Rxmk2D(I,J)*PHII(J) & +KAPPA*DEL1(I,J)*PHIR(J) ENDDO JNS( I,3*NT+1) = ZERO JNS( NT+I,3*NT+1) = dRxDdk1(I) JNS(2*NT+I,3*NT+1) = dRxDdk2(I) JNS( I,3*NT+2) = Rw(I) JNS( NT+I,3*NT+2) = RxwPhiR(I) JNS(2*NT+I,3*NT+2) = RxwPhiI(I) JNS(3*NT+1,I) = ZERO JNS(3*NT+1,NT+I) = VQ(I) JNS(3*NT+1,2*NT+I) = ZERO JNS(3*NT+2,I) = ZERO JNS(3*NT+2,NT+I) = ZERO JNS(3*NT+2,2*NT+I) = VQ(I) solNS(I) = -fvec0.PROG(I) solNS(NT+I) = -solNS(NT+I) solNS(2*NT+I) = -solNS(2*NT+I) ENDDO solNS(3*NT+1) = -(DDOT(NT,VQ,1,PHIR,1)-ONE) solNS(3*NT+2) = -(DDOT(NT,VQ,1,PHII,1)) JNS(3*NT+1,3*NT+1) = ZERO JNS(3*NT+1,3*NT+2) = ZERO JNS(3*NT+2,3*NT+1) = ZERO JNS(3*NT+2,3*NT+2) = ZERO * Resolution CALL DGESV(3*NT+2,1,JNS,3*NT+2,I4,solNS,3*NT+2,INFO) DO IOU=1,3*NT+2 ENDDO * Corrections dWB = solNS(3*NT+2) dK = solNS(3*NT+1) DO I = 1,NT dXB(I) = solNS(I) dPHIR(I)= solNS(NT+I) dPHII(I)= solNS(2*NT+I) Q1(I) = Q1(I) + dXB(I) PHIR(I) = PHIR(I) + dPHIR(I) PHII(I) = PHII(I) + dPHII(I) ENDDO OMEG = OMEG + dWB KAPPA= KAPPA+ dK * write(*,*) 'w +',dWB,'=',OMEG * write(*,*) 'kappa +',dK,'=',KAPPA END IF *======================================================================= * Evaluation du residu *======================================================================= * a) R CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA, & MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec0,fref, & .false.) *----------------------------------------------------------------------- IF (TYPC.EQ.'L') THEN *----------------------------------------------------------------------- * b) Rxs*Phi CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.) DO I=1,NT res2(I) = ZERO DO J=1,NT RX(I,J) = ZZ(I,J) - JAC(I,J) res2(I) = res2(I) + RX(I,J)*PHIR(J) ENDDO ENDDO * c) norm(Phi)= 1. res3 = DDOT(NT,PHIR,1,PHIR,1)-ONE * Y = ( R ; RQ * \phi ; \phi*\phi-1 ) arrRes(1) = dnrm2(NT,fvec0.prog(1),1) arrRes(2) = dnrm2(NT,res2,1) arrRes(3) = abs(res3) res = MAXVAL(arrRes,1) IF (IIMPI.GE.3) & WRITE(IOIMP,991) its,(arrRes(iou),iou=1,3),res 991 FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';', & E10.3,';','} -> res=',E10.3) ENDIF *----------------------------------------------------------------------- IF (TYPC.EQ.'B') THEN *----------------------------------------------------------------------- * b) Rxs*Phi CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.) CALL HBMDVEC(NT,NHBM,NDDL,Q1,ONE,V1) znormV1 = ONE / dnrm2(NT,V1,1) DO I = 1,NT V1(I) = V1(I)*znormV1 ENDDO DO I = 1,NT r_z = ZERO DO J = 1,NT RX(I,J) = ZZ(I,J) - JAC(I,J) RX2(I,J) = RX(I,J) + V1(I)*V1(J) r_z = r_z + RX2(I,J)*PHIR(J) ENDDO res2(I) = r_z ENDDO * c) Rw'*Phi CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.) IF (BAL.EQ.1) THEN DO I=1,NT Rw(I) = Rw(I) - TWO*OMEG*FEXA(I) ENDDO END IF res3 = DDOT(NT,Rw,1,PHIR,1) * d) norm(Phi) = 1 res4 = DDOT(NT,PHIR,1,PHIR,1)-ONE * Y = ( R ; RQ * \phi ; Rw\phi ; \phi*\phi-1 ) arrRes2(1) = dnrm2(NT,fvec0.prog(1),1) arrRes2(2) = dnrm2(NT,res2,1) arrRes2(3) = abs(res3) arrRes2(4) = abs(res4) res = MAXVAL(arrRes2,1) IF (IIMPI.GE.3) & WRITE(IOIMP,992) its,(arrRes2(iou),iou=1,4),res 992 FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';', & E10.3,';',E10.3,'} -> res=',E10.3) ENDIF *----------------------------------------------------------------------- IF (TYPC.EQ.'P') THEN *----------------------------------------------------------------------- CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.) DO I = 1,NT DO J = 1,NT RX(I,J) = ZZ(I,J) - JAC(I,J) ENDDO ENDDO CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D) CALL COLIMA(NT,NT,D10,D1w,ONE,OMEG,DEL1) * b) (Rx-((OMEG/2)**2)*D2)*PhiR - KAPPA*D1*PhiI * c) (Rx-((OMEG/2)**2)*D2)*PhiI + KAPPA*D1*PhiR DO I = 1,NT r_z1 = ZERO r_z2 = ZERO DO J = 1,NT r_z1 = r_z1 + Rxmk2D(I,J)*PHIR(J) & - KAPPA*DEL1(I,J)*PHII(J) r_z2 = r_z2 + Rxmk2D(I,J)*PHII(J) & + KAPPA*DEL1(I,J)*PHIR(J) ENDDO resV1(I) = r_z1 resV2(I) = r_z2 ENDDO * d) (VQ'*PHIR) - 1 res4 = DDOT(NT,VQ,1,PHIR,1)-ONE * Y = ( R ; [D0-k^2*D2]*\phi_R - k*D1*\phi_I ; * [D0-k^2*D2]*\phi_I + k*D1*\phi_R ; Q*\phi_I-1 ) arrRes2(1) = dnrm2(NT,fvec0.prog(1),1) arrRes2(2) = dnrm2(NT,resV1,1) arrRes2(3) = dnrm2(NT,resV2,1) arrRes2(4) = abs(res4) res = MAXVAL(arrRes2,1) IF (IIMPI.GE.3) & WRITE(IOIMP,993) its,(arrRes2(iou),iou=1,4),res 993 FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';', & E10.3,';',E10.3,'} -> res=',E10.3) ENDIF *----------------------------------------------------------------------- IF (TYPC.EQ.'N') THEN *----------------------------------------------------------------------- CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.) DO J = 1,NT DO I = 1,NT RX(I,J) = ZZ(I,J) - JAC(I,J) ENDDO ENDDO CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D) CALL COLIMA(NT,NT,D10,D1w,ONE,OMEG,DEL1) * b) (Rx-(KAPPA**2)*D2)*PhiR - KAPPA*D1*PhiI * c) (Rx-(KAPPA**2)*D2)*PhiI + KAPPA*D1*PhiR DO I = 1,NT r_z1 = ZERO r_z2 = ZERO DO J = 1,NT r_z1 = r_z1 + Rxmk2D(I,J)*PHIR(J) & - KAPPA*DEL1(I,J)*PHII(J) r_z2 = r_z2 + Rxmk2D(I,J)*PHII(J) & + KAPPA*DEL1(I,J)*PHIR(J) ENDDO resV1(I) = r_z1 resV2(I) = r_z2 ENDDO * d) (VQ'*PHIR) - 1 res4 = DDOT(NT,VQ,1,PHIR,1)-ONE * e) VQ'*PHII res5 = DDOT(NT,VQ,1,PHII,1) * Y = ( R ; [D0-k^2*D2]*\phi_R - k*D1*\phi_I ; * [D0-k^2*D2]*\phi_I + k*D1*\phi_R ; * vQ*\phi_R-1 ; vQ*\phi_I) arrRes3(1) = dnrm2(NT,fvec0.prog(1),1) arrRes3(2) = dnrm2(NT,resV1,1) arrRes3(3) = dnrm2(NT,resV2,1) arrRes3(4) = abs(res4) arrRes3(5) = abs(res5) res = MAXVAL(arrRes3,1) IF (IIMPI.GE.3) & WRITE(IOIMP,994) its,(arrRes3(iou),iou=1,5),res c & WRITE(IOIMP,994) its,KAPPA,(arrRes3(iou),iou=1,4),res 994 FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';', & E10.3,';',E10.3,';',E10.3,'} -> res=',E10.3) ENDIF *----------------------------------------------------------------------- c IF (IIMPI.GE.3) WRITE(IOIMP,*) 'Iter ',its,' -> res =',res IF (res.LE.TOLMIN) THEN CHECK=.false. CBIF = CBIF + 1 TYPBIF(CBIF) = TYPC DO I=1,NT QBIFU(I,CBIF) = Q1(I) QPSIR(I,CBIF) = PHIR(I) QPSII(I,CBIF) = PHII(I) Q1(NT) = QIN(I) ENDDO WBIFU(CBIF) = OMEG WBIF2(CBIF) = KAPPA OMEG = OMEGIN c RETURN GOTO 666 ELSE CHECK = .true. ENDIF ENDDO DO I = 1,NT Q1(I) = QIN(I) ENDDO OMEG = OMEGIN c message IF (CHECK) THEN c Pas de convergence apres %i1 iterations. L'execution continue INTERR(1)=MAXITS CALL ERREUR(151) c RETURN GOTO 666 ELSEIF (IIMPI.GE.2) THEN WRITE(IOIMP,*) '--> convergence apres ',its,'iterations' ENDIF * MENAGE 666 CONTINUE SEGSUP,MWORK SEGSUP,fvec,fvec0 RETURN END