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