1C * * * * * * * * * * * * * * * * * * * * * * * * * 2C --- DRIVER FOR RADAU5 AT VAN DER POL'S EQUATION 3C * * * * * * * * * * * * * * * * * * * * * * * * * 4C link dr1_radau5 radau5 decsol dc_decsol or 5C link dr1_radau5 radau5 lapack lapackc dc_lapack 6 IMPLICIT REAL*8 (A-H,O-Z) 7C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) 8 PARAMETER (ND=2,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) 9 DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK) 10 EXTERNAL FVPOL,JVPOL,SOLOUT 11C --- PARAMETER IN THE DIFFERENTIAL EQUATION 12 RPAR=1.0D-6 13C --- DIMENSION OF THE SYSTEM 14 N=2 15C --- COMPUTE THE JACOBIAN ANALYTICALLY 16 IJAC=1 17C --- JACOBIAN IS A FULL MATRIX 18 MLJAC=N 19C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM 20 IMAS=0 21C --- OUTPUT ROUTINE IS USED DURING INTEGRATION 22 IOUT=1 23C --- INITIAL VALUES 24 X=0.0D0 25 Y(1)=2.0D0 26 Y(2)=-0.66D0 27C --- ENDPOINT OF INTEGRATION 28 XEND=2.0D0 29C --- REQUIRED TOLERANCE 30 RTOL=1.0D-4 31 ATOL=1.0D0*RTOL 32 ITOL=0 33C --- INITIAL STEP SIZE 34 H=1.0D-6 35C --- SET DEFAULT VALUES 36 DO I=1,20 37 IWORK(I)=0 38 WORK(I)=0.D0 39 END DO 40C --- CALL OF THE SUBROUTINE RADAU5 41 CALL RADAU5(N,FVPOL,X,Y,XEND,H, 42 & RTOL,ATOL,ITOL, 43 & JVPOL,IJAC,MLJAC,MUJAC, 44 & FVPOL,IMAS,MLMAS,MUMAS, 45 & SOLOUT,IOUT, 46 & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) 47C --- PRINT FINAL SOLUTION 48 WRITE (6,99) X,Y(1),Y(2) 49 99 FORMAT(1X,'X =',F5.2,' Y =',2E18.10) 50C --- PRINT STATISTICS 51 WRITE (6,90) RTOL 52 90 FORMAT(' rtol=',D8.2) 53 WRITE (6,91) (IWORK(J),J=14,20) 54 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4,' accpt=',I4, 55 & ' rejct=',I3,' dec=',I4,' sol=',I5) 56 STOP 57 END 58C 59C 60 SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) 61C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "CONTR5" 62 IMPLICIT REAL*8 (A-H,O-Z) 63 DIMENSION Y(N),CONT(LRC) 64 COMMON /INTERN/XOUT 65 IF (NR.EQ.1) THEN 66 WRITE (6,99) X,Y(1),Y(2),NR-1 67 XOUT=0.2D0 68 ELSE 69 10 CONTINUE 70 IF (X.GE.XOUT) THEN 71C --- CONTINUOUS OUTPUT FOR RADAU5 72 WRITE (6,99) XOUT,CONTR5(1,XOUT,CONT,LRC), 73 & CONTR5(2,XOUT,CONT,LRC),NR-1 74 XOUT=XOUT+0.2D0 75 GOTO 10 76 END IF 77 END IF 78 99 FORMAT(1X,'X =',F5.2,' Y =',2E18.10,' NSTEP =',I4) 79 RETURN 80 END 81C 82C 83 SUBROUTINE FVPOL(N,X,Y,F,RPAR,IPAR) 84C --- RIGHT-HAND SIDE OF VAN DER POL'S EQUATION 85 IMPLICIT REAL*8 (A-H,O-Z) 86 DIMENSION Y(N),F(N) 87 F(1)=Y(2) 88 F(2)=((1-Y(1)**2)*Y(2)-Y(1))/RPAR 89 RETURN 90 END 91C 92C 93 SUBROUTINE JVPOL(N,X,Y,DFY,LDFY,RPAR,IPAR) 94C --- JACOBIAN OF VAN DER POL'S EQUATION 95 IMPLICIT REAL*8 (A-H,O-Z) 96 DIMENSION Y(N),DFY(LDFY,N) 97 DFY(1,1)=0.0D0 98 DFY(1,2)=1.0D0 99 DFY(2,1)=(-2.0D0*Y(1)*Y(2)-1.0D0)/RPAR 100 DFY(2,2)=(1.0D0-Y(1)**2)/RPAR 101 RETURN 102 END 103