C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT PENDULUM PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * C LINK dr_pendulum radau5m vecteurw decsol dc_decsol OR C LINK dr_pendulum radau5m vecteurw lapack lapackc dc_lapack IMPLICIT DOUBLE PRECISION (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=6,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=6 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=0.0D0 Y(1)=1.0D0 Y(2)=0.0D0 Y(3)=0.0D0 Y(4)=0.0D0 Y(5)=0.0D0 Y(6)=0.0D0 C --- ENDPOINT OF INTEGRATION XEND=10.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)=4 IWORK(6)=2 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,99) X,Y(3),Y(4) WRITE (6,98) X,Y(5),Y(6) 99 FORMAT(1X,'X =',F5.2,' Y =',2E22.14) 98 FORMAT(1X,'X =',F5.2,' Z =',2E22.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 (NFICY,99) X,Y(3),Y(4),NR-1 WRITE (NFICZ,99) X,Y(5),Y(6),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 (NFICY,99) XOUTY,CONTR5(3,XOUTY,CONT,LRC), & CONTR5(4,XOUTY,CONT,LRC),NR-1 WRITE (NFICZ,99) XOUTY,CONTR5(5,XOUTY,CONT,LRC), & CONTR5(6,XOUTY,CONT,LRC),NR-1 XOUTY=XOUTY+0.2D0 GOTO 10 END IF END IF 99 FORMAT(1X,F5.2,2E22.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,99) X,Y(5),Y(6),NR-1 WRITE(NFICY,99) X,Y(1),Y(2),NR-1 WRITE(NFICY,99) X,Y(3),Y(4),NR-1 XOUTZ=0.2D0 XOUTY=0.2D0 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) YN5= W1*VALG(5)+W2*VALG(23)+W3*VALG(41) YN5=YN5+ W4*VALG(11)+W5*VALG(29)+W6*VALG(47) YN5=YN5+ W7*VALG(17)+W8*VALG(35)+W9*VALG(53) YN6= W1*VALG(6)+W2*VALG(24)+W3*VALG(42) YN6=YN6+W4*VALG(12)+W5*VALG(30)+W6*VALG(48) YN6=YN6+W7*VALG(18)+W8*VALG(36)+W9*VALG(54) WRITE(NFICZ,99) XOUTZ,YN5,YN6,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(19)+WY3*VALG(37) YN1=YN1+WY4*VALG(7)+WY5*VALG(25)+WY6*VALG(43) YN2=WY1*VALG(2)+WY2*VALG(20)+WY3*VALG(38) YN2=YN2+WY4*VALG(8)+WY5*VALG(26)+WY6*VALG(44) YN3=WY1*VALG(3)+WY2*VALG(21)+WY3*VALG(39) YN3=YN3+WY4*VALG(9)+WY5*VALG(27)+WY6*VALG(45) YN4=WY1*VALG(4)+WY2*VALG(22)+WY3*VALG(40) YN4=YN4+WY4*VALG(10)+WY5*VALG(28)+WY6*VALG(46) 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(7)+WY2*VALG(25)+WY3*VALG(43) YN1=YN1+WY4*VALG(13)+WY5*VALG(31)+WY6*VALG(49) YN2=WY1*VALG(8)+WY2*VALG(26)+WY3*VALG(44) YN2=YN2+WY4*VALG(14)+WY5*VALG(32)+WY6*VALG(50) YN3=WY1*VALG(9)+WY2*VALG(27)+WY3*VALG(45) YN3=YN3+WY4*VALG(15)+WY5*VALG(33)+WY6*VALG(51) YN4=WY1*VALG(10)+WY2*VALG(28)+WY3*VALG(46) YN4=YN4+WY4*VALG(16)+WY5*VALG(34)+WY6*VALG(52) ENDIF WRITE(NFICY,99) XOUTY,YN1,YN2,NR-1 WRITE(NFICY,99) XOUTY,YN3,YN4,NR-1 XOUTY=XOUTY+0.2D0 GOTO 52 ENDIF ENDIF RETURN 99 FORMAT(1X,F5.2,2E22.14,I4) END C C SUBROUTINE FEX(N,X,Y,F,RPAR,IPAR) C --- RIGHT-HAND SIDE OF PENDULUM PROBLEM IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),F(N) F(1) = Y(3)-Y(1)*Y(6) F(2) = Y(4)-Y(2)*Y(6) F(3) =-Y(1)*Y(5) F(4) =-Y(2)*Y(5)-1.0D0 F(5) = Y(1)*Y(1)+Y(2)*Y(2)-1.0D0 F(6) = Y(1)*Y(3)+Y(2)*Y(4) RETURN END C C SUBROUTINE JEX(N,X,Y,DFY,LDFY,RPAR,IPAR) C --- JACOBIAN OF PENDULUM PROBLEM IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) DO I=1,LDFY DO J=1,N DFY(I,J)=0.0D0 END DO END DO DFY(1,1)=-Y(6) DFY(1,3)= 1.0D0 DFY(1,6)=-Y(1) DFY(2,2)=-Y(6) DFY(2,4)= 1.0D0 DFY(2,6)=-Y(2) DFY(3,1)=-Y(5) DFY(3,5)=-Y(1) DFY(4,2)=-Y(5) DFY(4,5)=-Y(2) DFY(5,1)=2*Y(1) DFY(5,2)=2*Y(2) DFY(6,1)= Y(3) DFY(6,2)= Y(4) DFY(6,3)= Y(1) DFY(6,4)= Y(2) 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)=1.0D0 AM(1,4)=1.0D0 AM(1,5)=0.0D0 AM(1,6)=0.0D0 RETURN END