C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT CERCLE PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * C LINK dr_cercle radau5m vecteurw decsol dc_decsol OR C LINK dr_cercle radau5m vecteurw lapack lapackc dc_lapack IMPLICIT DOUBLE PRECISION (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=3,LWORK=4*ND*ND+23*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK) CHARACTER*20 NOMFICY,NOMFICZ EXTERNAL FEX,JEX,SOLOUT,SOLOUTZ,MAS C --- DEFINITION OF THE RESULT FILES WRITE(6,*)'NAME OF THE RESULT FILE FOR THE DIFFERENTIAL COMPONENTS?' READ(*,*) NOMFICY WRITE(6,*)'THE RESULT FILE FOR Y IS : ',NOMFICY NFICY = 10 OPEN(NFICY,FORM='FORMATTED',STATUS='UNKNOWN',FILE=NOMFICY) WRITE(6,*)'NAME OF THE RESULT FILE FOR THE ALGEBRAIC COMPONENTS?' READ(*,*) NOMFICZ WRITE(6,*)'THE RESULT FILE FOR Z IS : ',NOMFICZ NFICZ = 11 OPEN(NFICZ,FORM='FORMATTED',STATUS='UNKNOWN',FILE=NOMFICZ) C --- PARAMETER IN THE DIFFERENTIAL EQUATION IPAR=0 C --- DIMENSION OF THE SYSTEM N=3 C --- COMPUTE THE JACOBIAN ANALYTICALLY IJAC=1 C --- JACOBIAN IS A FULL MATRIX MLJAC=N C --- DIFFERENTIAL EQUATION IS IN IMPLICIT FORM IMAS=1 MLMAS=0 MUMAS=0 C --- OUTPUT ROUTINE IS USED DURING INTEGRATION IOUT=1 C --- INITIAL VALUES X=-1.0D0 Y(1)=1.0D0 Y(2)=0.0D0 Y(3)=0.0D0 C --- ENDPOINT OF INTEGRATION XEND=11.0D0 C --- REQUIRED TOLERANCE WRITE(*,*)'RTOL?' READ(*,*)RTOL ATOL=1.0D0*RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-7 C --- SET DEFAULT VALUES DO I=1,20 IWORK(I)=0 WORK(I)=0.D0 END DO IWORK(5)=2 IWORK(6)=1 C --- CALL OF THE SUBROUTINE RADAU5 CALL RADAU5M(N,FEX,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JEX,IJAC,MLJAC,MUJAC, & MAS,IMAS,MLMAS,MUMAS, & SOLOUT,SOLOUTZ,IOUT,NFICY,NFICZ, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- CLOSE THE RESULT FILES CLOSE(NFICY) CLOSE(NFICZ) C --- PRINT FINAL SOLUTION WRITE (6,99) X,Y(1),Y(2) WRITE (6,98) X,Y(3) 99 FORMAT(1X,'X =',F5.2,' Y =',2E22.14) 98 FORMAT(1X,'X =',F5.2,' Z =',E22.14) C --- PRINT STATISTICS WRITE (6,90) RTOL 90 FORMAT(' RTOL=',D8.2) WRITE (6,91) (IWORK(J),J=14,20) 91 FORMAT(' FCN=',I5,' JAC=',I4,' STEP=',I4,' ACCPT=',I4, & ' REJCT=',I3,' DEC=',I4,' SOL=',I5) STOP END C C SUBROUTINE SOLOUT (NFICY,NFICZ,NR,XOLD,X,Y,CONT,LRC,N, & RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "CONTR5" IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUTY,XOUTZ IF (NR.EQ.1) THEN WRITE (NFICY,99) X,Y(1),Y(2),NR-1 WRITE (NFICZ,98) X,Y(3),NR-1 XOUTY=-0.8D0 ELSE 10 CONTINUE IF (X.GE.XOUTY) THEN C --- CONTINUOUS OUTPUT FOR RADAU5M WRITE (NFICY,99) XOUTY,CONTR5(1,XOUTY,CONT,LRC), & CONTR5(2,XOUTY,CONT,LRC),NR-1 WRITE (NFICZ,98) XOUTY,CONTR5(3,XOUTY,CONT,LRC), & NR-1 XOUTY=XOUTY+0.2D0 GOTO 10 END IF END IF 99 FORMAT(1X,F5.2,2E22.14,I4) 98 FORMAT(1X,F5.2,E22.14,I4) RETURN END C C SUBROUTINE SOLOUTZ (NFICY,NFICZ,NR,XOLD,X,Y,N,VALG,LRV, & RPAR,IPAR,IRTRN,ICAS,R1,R2,HTOT,H2,H3,LAST) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "VECTWY" AND "VECTWZ" IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N) DIMENSION VWT(9),WY(6),VALG(LRV) COMMON/INTERN/XOUTY,XOUTZ EXTERNAL VECTWZ,VECTWY IF (NR.EQ.1) THEN WRITE(NFICZ,98) X,Y(3),NR-1 WRITE(NFICY,99) X,Y(1),Y(2),NR-1 XOUTZ=-0.8D0 XOUTY=-0.8D0 ELSE IF (LAST) THEN XPAR=X ELSE XPAR=XOLD ENDIF 51 CONTINUE IF (XPAR.GE.XOUTZ) THEN T=(XOUTZ+HTOT-X)/HTOT CALL VECTWZ(ICAS,VWT,R1,R2,T) W1=VWT(1) W2=VWT(2) W3=VWT(3) W4=VWT(4) W5=VWT(5) W6=VWT(6) W7=VWT(7) W8=VWT(8) W9=VWT(9) YN3=W1*VALG(3)+W2*VALG(12)+W3*VALG(21)+ & W4*VALG(6)+W5*VALG(15)+W6*VALG(24)+ & W7*VALG(9)+W8*VALG(18)+W9*VALG(27) WRITE(NFICZ,98) XOUTZ,YN3,NR-1 XOUTZ=XOUTZ+0.2D0 GOTO 51 ENDIF 52 CONTINUE IF (X.GE.XOUTY) THEN IF ((XOLD-H2).GE.XOUTY) THEN HTOTY=HTOT-H3 T=(XOUTY-XOLD+HTOTY)/HTOTY R=(HTOT-H3-H2)/HTOTY CALL VECTWY(WY,R,T) WY1=WY(1) WY2=WY(2) WY3=WY(3) WY4=WY(4) WY5=WY(5) WY6=WY(6) YN1=WY1*VALG(1)+WY2*VALG(10)+WY3*VALG(19) YN1=YN1+WY4*VALG(4)+WY5*VALG(13)+WY6*VALG(22) YN2=WY1*VALG(2)+WY2*VALG(11)+WY3*VALG(20) YN2=YN2+WY4*VALG(5)+WY5*VALG(14)+WY6*VALG(23) ELSE HTOTY=H2+H3 T=(XOUTY-X+HTOTY)/HTOTY R=H2/HTOTY CALL VECTWY(WY,R,T) WY1=WY(1) WY2=WY(2) WY3=WY(3) WY4=WY(4) WY5=WY(5) WY6=WY(6) YN1=WY1*VALG(4)+WY2*VALG(13)+WY3*VALG(22) YN1=YN1+WY4*VALG(7)+WY5*VALG(16)+WY6*VALG(25) YN2=WY1*VALG(5)+WY2*VALG(14)+WY3*VALG(23) YN2=YN2+WY4*VALG(8)+WY5*VALG(17)+WY6*VALG(26) ENDIF WRITE(NFICY,99) XOUTY,YN1,YN2,NR-1 XOUTY=XOUTY+0.2D0 GOTO 52 ENDIF ENDIF RETURN 99 FORMAT(1X,F5.2,2E22.14,I4) 98 FORMAT(1X,F5.2,E22.14,I4) END C C DOUBLE PRECISION FUNCTION PHI(X) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PI = 4.0D0*DATAN(1.0D0) EPS0=-1.0D0 EPS1=1.0D0 EPS2=4.0D0 EPS3=6.0D0 EPS4=9.0D0 EPS5=11.0D0 PHI=0.0D0 IF ((X.LT.EPS1).AND.(X.GT.EPS0))THEN X2=X*X AUX=1.0D0/(X2-1.0D0) PHI=0.5D0*PI*DEXP(X2*AUX) ENDIF IF ((X.LT.EPS3).AND.(X.GT.EPS2)) THEN XX=X-5.0D0 X2=XX*XX AUX=1.0D0/(X2-1.0D0) PHI=0.5D0*PI*DEXP(X2*AUX) ENDIF IF ((X.LT.EPS5).AND.(X.GT.EPS4))THEN XX=X-10.0D0 X2=XX*XX AUX=1.0D0/(X2-1.0D0) PHI=0.5D0*PI*DEXP(X2*AUX) ENDIF RETURN END DOUBLE PRECISION FUNCTION PHID(X) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PI = 4.0D0*DATAN(1.0D0) EPS0=-1.0D0 EPS1=1.0D0 EPS2=4.0D0 EPS3=6.0D0 EPS4=9.0D0 EPS5=11.0D0 PHID=0.0D0 IF ((X.LT.EPS1).AND.(X.GT.EPS0))THEN X2=X*X AUX=1.0D0/(X2-1.0D0) EAUX=DEXP(X2*AUX) PHID=-PI*X*AUX*AUX*EAUX ENDIF IF ((X.LT.EPS3).AND.(X.GT.EPS2)) THEN XX=X-5.0D0 X2=XX*XX AUX=1.0D0/(X2-1.0D0) EAUX=DEXP(X2*AUX) PHID=-PI*XX*AUX*AUX*EAUX ENDIF IF ((X.LT.EPS5).AND.(X.GT.EPS4))THEN XX=X-10.0D0 X2=XX*XX AUX=1.0D0/(X2-1.0D0) EAUX=DEXP(X2*AUX) PHID=-PI*XX*AUX*AUX*EAUX ENDIF RETURN END C C SUBROUTINE FEX(N,X,Y,F,RPAR,IPAR) C --- RIGHT-HAND SIDE OF CERCLE PROBLEM IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),F(N) EXTERNAL PHID DPHI= PHID(X) F(1)=-DPHI*Y(2)+Y(3)*Y(1) F(2)= DPHI*Y(1)+Y(3)*Y(2) F(3)= Y(1)*Y(1)+Y(2)*Y(2)-1.0D0 RETURN END C C SUBROUTINE JEX(N,X,Y,DFY,LDFY,RPAR,IPAR) C --- JACOBIAN OF CERCLE PROBLEM IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) EXTERNAL PHID DPHI= PHID(X) DFY(1,1)= Y(3) DFY(1,2)=-DPHI DFY(1,3)= Y(1) DFY(2,1)= DPHI DFY(2,2)= Y(3) DFY(2,3)= Y(2) DFY(3,1)= 2*Y(1) DFY(3,2)= 2*Y(2) DFY(3,3)= 0.0D0 RETURN END C C SUBROUTINE MAS(N, AM, LMAS,RPAR,IPAR) C --- MASS MATRIX IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION AM(LMAS,N) AM(1,1)=1.0D0 AM(1,2)=1.0D0 AM(1,3)=0.0D0 RETURN END