C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT VDPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=2,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) + ,true(nd),error(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) EXTERNAL FVANDER,JVANDER,SOLOUT c ------ FILE DE DONNEES ---------- write(6,3500) 3500 format(1x,'radau5 on van der pol') RPAR=1.D-6 C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM N=2 C --- COMPUTE THE JACOBIAN ANALYTICALLY IJAC=1 C --- JACOBIAN IS A FULL MATRIX MLJAC=N C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 C --- INITIAL VALUES X=0.0D0 Y(1)=2.0D0 Y(2)=0.0D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=1.0D0 it1=mclock() DO 20 I=1,11 C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) CALL RADAU5(N,FVANDER,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JVANDER,IJAC,MLJAC,MUJAC, & FVANDER,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c WRITE (8,*) Y(1) c WRITE (8,*) Y(2) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO XEND=XEND+1.D0 if(i.eq.1) then true(1)=-0.1863646254808130E+01 true(2)= 0.7535430865435460E+00 else if(i.eq.2) then true(1)= 0.1706167732170456E+01 true(2)= -0.8928097010248257E+00 else if(i.eq.3) then true(1)= -0.1510606936744095E+01 true(2)= 0.1178380000730945E+01 else if(i.eq.4) then true(1)= 0.1194414677695905E+01 true(2)= -0.2799585996540082E+01 else if(i.eq.5) then true(1)= 0.1890428596416747E+01 true(2)= -0.7345118680166940E+00 else if(i.eq.6) then true(1)=-0.1737716306805883E+01 true(2)= 0.8604008653025923E+00 else if(i.eq.7) then true(1)= 0.1551614645548223E+01 true(2)= -0.1102382892533321E+01 else if(i.eq.8) then true(1)=-0.1278631984330405E+01 true(2)= 0.2013890883009155E+01 else if(i.eq.9) then true(1)= -0.1916552949489830E+01 true(2)= 0.7169573003463228E+00 else if(i.eq.10) then true(1)= 0.1768163792391936E+01 true(2)=-0.8315276407898496E+00 else if(i.eq.11) then true(1)=-0.1590150544829062E+01 true(2)= 0.1040279389212485E+01 endif sum=0.0D+0 do 270 k=1,2 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 continue it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) WRITE(6,*)TARRAY(1),summ,sumh/11.0d+0 summ=summ*atol sumh=sumh*atol/(11.0D+0) write(6,*) summ,sumh WRITE (6,*)(ISTAT(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) RETURN END SUBROUTINE FVANDER(N,X,Y,F,EPS,IPAR) C --- RIGHT-HAND SIDE OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),F(N) F(1)=Y(2) PROD=1.D0-Y(1)*Y(1) F(2)=(PROD*Y(2)-Y(1))/EPS RETURN END C SUBROUTINE JVANDER(N,X,Y,DFY,LDFY,EPS,IPAR) C --- JACOBIAN OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) DFY(1,1)=0.D0 DFY(1,2)=1.D0 DFY(2,1)=(-2.0D0*Y(1)*Y(2)-1.0D0)/EPS DFY(2,2)=(1.0D0-Y(1)**2)/EPS RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT ROBERTSON PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=3,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20),error(nd) + ,true(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) EXTERNAL FROBER,JROBER,SOLOUT c ------ FILE DE DONNEES ---------- OPEN(8,FILE='res_rad5') REWIND 8 write(8,3050) 3050 format(1x,'these are the robertson results by radau') C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+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 EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 C --- INITIAL VALUES X=0.0D0 Y(1)=1.0D0 Y(2)=0.0D0 Y(3)=0.0D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=1.0D-6*RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=1.0D0 c CALL DTIME(TARRAY) it1=mclock() DO 20 I=1,12 C --- CALL OF THE SUBROUTINE RADAU5 CALL RADAU5(N,FROBER,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JROBER,IJAC,MLJAC,MUJAC, & FROBER,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c WRITE (8,*) Y(1) c WRITE (8,*) Y(2) c WRITE (8,*) Y(3) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO XEND=XEND*10.D0 if(i.eq.1) then true(1)= 0.9664597373330035E+00 true(2)=0.3074626578578675E-04 true(3)=0.3350951640121071E-01 else if(i.eq.2) then true(1)=0.8413699238414729E+00 true(2)=0.1623390937990473E-04 true(3)=0.1586138422491472E+00 else if(i.eq.3) then true(1)=0.6172348823960878E+00 true(2)=0.6153591274639123E-05 true(3)=0.3827589640126376E+00 else if(i.eq.4) then true(1)=0.3368745306607069E+00 true(2)=0.2013702318261393E-05 true(3)=0.6631234556369748E+00 else if(i.eq.5) then true(1)=0.1073004285378040E+00 true(2)=0.4800166972571660E-06 true(3)=0.8926990914454987E+00 else if(i.eq.6) then true(1)=0.1786592114209946E-01 true(2)=0.7274751468436319E-07 true(3)=0.9821340061103859E+00 else if(i.eq.7) then true(1)= 0.2031483924973415E-02 true(2)=0.8142277783356159E-08 true(3)=0.9979685079327488E+00 else if(i.eq.8) then true(1)=0.2076093439016395E-03 true(2)=0.8306077485067610E-09 true(3)=0.9997923898254906E+00 else if(i.eq.9) then true(1)=0.2082417512179460E-04 true(2)=0.8329841429908955E-10 true(3)=0.9999791757415798E+00 else if(i.eq.10) then true(1)=0.2083229471647004E-05 true(2)=0.8332935037760723E-11 true(3)=0.9999979167621954E+00 else if(i.eq.11) then true(1)=0.2083328471883087E-06 true(2)=0.8333315602809495E-12 true(3)=0.9999997916663195E+00 else if(i.eq.12) then true(1)=0.2083340149701284E-07 true(2)=0.8333360770334744E-13 true(3)=0.9999999791665152E+00 end if sum=0.0d+0 do 270 k=1,3 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() write(6,*) y(1) write(6,*) y(2) write(6,*) y(3) tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh/12.0D+0 summ=summ*atol sumh=sumh*atol/12.0D+0 write(6,*) SUMM,SUMH WRITE (8,*)(ISTAT(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT WRITE (6,99) X,Y(1),Y(2),Y(3),NR-1 99 FORMAT(1X,'X =',F5.2,' Y =',3E18.10,' NSTEP =',I4) RETURN END SUBROUTINE FROBER(N,X,Y,F,RPAR,IPAR) C --- RIGHT-HAND SIDE OF ROBERTSON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),F(N) F(1)=-0.04D0*Y(1)+1.D4*Y(2)*Y(3) F(3)=3.D7*Y(2)*Y(2) F(2)=-F(1)-F(3) RETURN END C SUBROUTINE JROBER(N,X,Y,DFY,LDFY,RPAR,IPAR) C --- JACOBIAN OF ROBERTSON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) PROD1=1.0D4*Y(2) PROD2=1.0D4*Y(3) PROD3=6.0D7*Y(2) DFY(1,1)=-0.04D0 DFY(1,2)=PROD2 DFY(1,3)=PROD1 DFY(2,1)=0.04D0 DFY(2,2)=-PROD2-PROD3 DFY(2,3)=-PROD1 DFY(3,1)=0.D0 DFY(3,2)=PROD3 DFY(3,3)=0.D0 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 C * * * * * * * * * * * * * * * * * * * * * * * * * compile equation compile //venus/user/hairer/hailib/time cfeh driver_radau5 ~book2/appendix/radau5 ~book2/appendix/decsol equation time IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=8,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) dimension true(nd),error(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) EXTERNAL FHIRES,JHIRES,SOLOUT c ------ FILE DE DONNEES ---------- OPEN(8,FILE='res_rad5') write(8,3050) 3050 format(1x,'results on Hires problem') REWIND 8 C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP sumh=0.0D+0 summ=0.0D+0 C --- DIMENSION OF THE SYSTEM N=8 C --- COMPUTE THE JACOBIAN IJAC=1 C --- JACOBIAN IS A FULL MATRIX MLJAC=N C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 C --- INITIAL VALUES X=0.0D0 Y (1) = 1.D0 Y (2) = 0.D0 Y (3) = 0.D0 Y (4) = 0.D0 Y (5) = 0.D0 Y (6) = 0.D0 Y (7) = 0.D0 Y (8) = 0.0057D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL*1.D-4 ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=321.8122D0 c CALL DTIME(TARRAY) it1=mclock() DO 20 I=1,2 C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) CALL RADAU5(N,FHIRES,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JHIRES,IJAC,MLJAC,MUJAC, & FHIRES,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c DO 15 K=1,8 c 15 WRITE (8,*)Y(K) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO XEND=XEND+100.D0 c CALL DTIME(TARRAY) if(i.eq.1) then true(1)=0.000737131257332567 true(2)=0.000144248572631618 true(3)=0.000058887297409676 true(4)=0.001175651343283149 true(5)=0.002386356198831330 true(6)=0.006238968252742796 true(7)=0.002849998395185769 true(8)=0.002850001604814231 else if(i.eq.2) then true(1)=0.000670305503581864 true(2)=0.000130996846986347 true(3)=0.000046862231597733 true(4)=0.001044668020551705 true(5)=0.000594883830951485 true(6)=0.001399628833942774 true(7)=0.001014492757718480 true(8)=0.004685507242281520 end if sum=0.0D+0 do 270 k=1,8 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() do 15 k=1,8 write(6,*) y(k) 15 continue tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh/2.0D+0 summ=summ*atol sumh=sumh*atol/2.0D+0 write(6,*) summ,sumh WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- WRITE (6,*)(ISTAT(J),J=14,20) WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) RETURN END C SUBROUTINE FHIRES (N, X, Y, DY, RPAR, IPAR) INTEGER N DOUBLE PRECISION X, Y, DY DIMENSION Y (8), DY (8) DY (1) = -1.71D0*Y(1) + 0.43D0*Y(2) + 8.32D0*Y(3) + 0.0007D0 DY (2) = 1.71D0*Y(1) - 8.75D0*Y(2) DY (3) = -10.03D0*Y(3) + 0.43D0*Y(4) + 0.035D0*Y(5) DY (4) = 8.32D0*Y(2) + 1.71D0*Y(3) - 1.12D0*Y(4) DY (5) = -1.745D0*Y(5) + 0.43D0*Y(6) + 0.43D0*Y(7) DY (6) = -280.0D0*Y(6)*Y(8) + 0.69D0*Y(4) + 1.71D0*Y(5) - & 0.43D0*Y(6) + 0.69D0*Y(7) DY (7) = 280.0D0*Y(6)*Y(8) - 1.81D0*Y(7) DY (8) = -DY (7) RETURN END SUBROUTINE JHIRES(N,X,Y,DFY,LDFY,RPAR,IPAR) C -------- JACOBIAN FOR HIRES PROBLEM -------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) C------ METTRE A ZERO ------- DO 1 I=1,N DO 1 J=1,N 1 DFY(I,J)=0.D0 C DFY(1,1)= -1.71D0 DFY(1,2)= 0.43D0 DFY(1,3)= + 8.32D0 C DFY(2,1)= 1.71D0 DFY(2,2)= - 8.75D0 C DFY(3,3)= -10.03D0 DFY(3,4)= 0.43D0 DFY(3,5)= + 0.035D0 C DFY(4,2)= 8.32D0 DFY(4,3)= + 1.71D0 DFY(4,4)= - 1.12D0 C DFY(5,5)= -1.745D0 DFY(5,6)= + 0.43D0 DFY(5,7)= + 0.43D0 C DFY(6,4)= + 0.69D0 DFY(6,5)= + 1.71D0 DFY(6,6)= - 0.43D0 -280.0D0*Y(8) DFY(6,7)= + 0.69D0 DFY(6,8)= -280.0D0*Y(6) C DFY(7,6)= 280.0D0*Y(8) DFY(7,7)= - 1.81D0 DFY(7,8)= 280.0D0*Y(6) C DFY(8,6)= - 280.0D0*Y(8) DFY(8,7)= 1.81D0 DFY(8,8)= - 280.0D0*Y(6) RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT OREGO PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * compile equation compile //venus/user/hairer/hailib/time cfeh driver_radau5 ~book2/appendix/radau5 ~book2/appendix/decsol equation time IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=3,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) dimension true(nd),error(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) EXTERNAL FOREGON,JOREGON,SOLOUT c ------ FILE DE DONNEES ---------- OPEN(8,FILE='res_rad5') write(8,3050) 3050 format(1x,'radau problem on oreg') REWIND 8 C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP sumh=0.0D+0 summ=0.0D+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 EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 C --- INITIAL VALUES X=0.0D0 Y(1)=1.0D0 Y(2)=2.0D0 Y(3)=3.0D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=1.0D-6*RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=30.0D0 c CALL DTIME(TARRAY) it1=mclock() DO 20 I=1,12 C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) CALL RADAU5(N,FOREGON,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JOREGON,IJAC,MLJAC,MUJAC, & FOREGON,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c WRITE (8,*) Y(1) c WRITE (8,*) Y(2) c WRITE (8,*) Y(3) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO XEND=XEND+30.D0 c CALL DTIME(TARRAY) if(i.eq.1) then true(1)=0.1000661467180497E+01 true(2)=0.1512778937348249E+04 true(3)=0.1035854312767229E+05 else if(i.eq.2) then true(1)=0.1000874625199626E+01 true(2)=0.1144336972384497E+04 true(3)=0.8372149966624639E+02 else if(i.eq.3) then true(1)=0.1001890368438751E+01 true(2)=0.5299926232295553E+03 true(3)=0.1662279579042420E+01 else if(i.eq.4) then true(1)=0.1004118022612645E+01 true(2)=0.2438326079910346E+03 true(3)=0.1008822224048647E+01 else if(i.eq.5) then true(1)=0.1008995416634061E+01 true(2)=0.1121664388662539E+03 true(3)=0.1007783229065319E+01 else if(i.eq.6) then true(1)=0.1019763472537298E+01 true(2)=0.5159761322947535E+02 true(3)=0.1016985778956374E+01 else if(i.eq.7) then true(1)= 0.1043985088527474E+01 true(2)=0.2373442027531524E+02 true(3)=0.1037691843544522E+01 else if(i.eq.8) then true(1)=0.1100849071667922E+01 true(2)=0.1091533805469020E+02 true(3)=0.1085831969810860E+01 else if(i.eq.9) then true(1)=0.1249102130020572E+01 true(2)=0.5013945178605446E+01 true(3)=0.1208326626237875E+01 else if(i.eq.10) then true(1)=0.1779724751937019E+01 true(2)=0.2281852385542403E+01 true(3)=0.1613754023671725E+01 else if(i.eq.11) then true(1)=0.1000889326903503E+01 true(2)=0.1125438585746596E+04 true(3)=0.1641049483777168E+05 else if(i.eq.12) then true(1)=0.1000814870318523E+01 true(2)=0.1228178521549889E+04 true(3)=0.1320554942846513E+03 endif sum=0.0D+0 do 270 k=1,3 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh/12.0D+0 summ=summ*atol sumh=sumh*atol/12.0D+0 write(6,*) summ,sumh ID=0 WRITE (8,*)(ISTAT(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT WRITE (6,99) X,Y(1),Y(2),Y(3),NR-1 99 FORMAT(1X,'X =',F11.4,' Y =',3F16.5,' NSTEP =',I4) RETURN END SUBROUTINE FOREGON(N,X,Y,F,RPAR,IPAR) C --- RIGHT-HAND SIDE OF OREGON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),F(N) F(1)=77.27D0*(Y(2)+Y(1)*(1.D0-8.375D-6*Y(1)-Y(2))) F(2)=(Y(3)-(1.D0+Y(1))*Y(2))/77.27D0 F(3)=0.161D0*(Y(1)-Y(3)) RETURN END C SUBROUTINE JOREGON(N,X,Y,DFY,LDFY,RPAR,IPAR) C --- JACOBIAN OF OREGON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) DFY(1,1)=77.27D0*(1.D0-2.D0*8.375D-6*Y(1)-Y(2)) DFY(1,2)=77.27D0*(1.D0-Y(1)) DFY(1,3)=0.D0 DFY(2,1)=-Y(2)/77.27D0 DFY(2,2)=-(1.D0+Y(1))/77.27D0 DFY(2,3)=1.D0/77.27D0 DFY(3,1)=.161D0 DFY(3,2)=.0D0 DFY(3,3)=-.161D0 RETURN END C --- DRIVER FOR RADAU5 AT ROBERTSON PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=4,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) dimension true(nd),error(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) EXTERNAL FE5,JE5,SOLOUT c ------ FILE DE DONNEES ---------- write(6,3050) 3050 format(1x,'results on E5 by radau') C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP sumh=0.0D+0 summ=0.0D+0 C --- DIMENSION OF THE SYSTEM N=4 C --- COMPUTE THE JACOBIAN ANALYTICALLY IJAC=1 C --- JACOBIAN IS A FULL MATRIX MLJAC=N C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 C --- INITIAL VALUES X=0.0D0 Y(1)=1.76D-3 Y(2)=0.0D0 Y(3)=0.0D0 Y(4)=0.0D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=1.7D-24 ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=10.0D0 c CALL DTIME(TARRAY) it1=mclock() DO 20 I=1,7 C --- CALL OF THE SUBROUTINE RADAU5 CALL RADAU5(N,FE5,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JE5,IJAC,MLJAC,MUJAC, & FE5,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c WRITE (6,*) Y(1) c WRITE (6,*) Y(2) c WRITE (6,*) Y(3) c WRITE (6,*) Y(4) c write (6,*) y(2), y(3)+y(4) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO XEND=XEND*100.D0 c CALL DTIME(TARRAY) if(i.eq.1) then true(1)=1.7599259497677897058D-003 true(2)=1.3846281519376516449D-011 true(3)=7.6370038530073911180D-013 true(4)=1.3082581134075777338D-011 else if(i.eq.2) then true(1)=1.6180769999072942552D-003 true(2)=1.3822370304983735443D-010 true(3)=8.2515735006838336088D-012 true(4)=1.2997212954915352082D-010 else if(i.eq.3) then true(1)=7.4813208224292220114D-006 true(2)=2.3734781561205975019D-012 true(3)=2.2123586689581663654D-012 true(4)=1.6111948716243113653D-013 else if(i.eq.4) then true(1)=4.7150333630401632232D-010 true(2)=1.8188895860807021729D-014 true(3)=1.8188812376786725407D-014 true(4)=8.3484020296321693074D-020 else if(i.eq.5) then true(1)=3.1317148329356996037D-014 true(2)=1.4840957952870064294D-016 true(3)=1.4840957948345691466D-016 true(4)=4.5243728279782625194D-026 else if(i.eq.6) then true(1)=3.8139035189787091771D-049 true(2)=1.0192582567660293322D-020 true(3)=1.0192582567660293322D-020 true(4)=3.7844935507486221171D-065 else if(i.eq.7) then true(1)=0.0000000000000000000D-000 true(2)=8.8612334976263783420D-023 true(3)=8.8612334976263783421D-023 true(4)=0.0000000000000000000D-000 end if sum=0.0D+0 do 270 k=1,4 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh/7.0D+0 summ=summ*rtol sumh=sumh*rtol/7.0D+0 write(6,*) summ,sumh ID=0 WRITE (6,*)(ISTAT(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT WRITE (6,99) X,Y(1),Y(2),Y(3),NR-1 99 FORMAT(1X,'X =',F5.2,' Y =',3E18.10,' NSTEP =',I4) RETURN END SUBROUTINE FE5(N,X,Y,F,RPAR,IPAR) C --- RIGHT-HAND SIDE OF ROBERTSON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),F(N) PROD1=7.89D-10*Y(1) PROD2=1.1D7*Y(1)*Y(3) PROD3=1.13D9*Y(2)*Y(3) PROD4=1.13D3*Y(4) F(1)=-PROD1-PROD2 F(2)=PROD1-PROD3 F(4)=PROD2-PROD4 F(3)=F(2)-F(4) RETURN END C SUBROUTINE JE5(N,X,Y,DFY,LDFY,RPAR,IPAR) C --- JACOBIAN OF ROBERTSON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) A=7.89D-10 B=1.1D7 CM=1.13D9 C=1.13D3 DFY(1,1)=-A-B*Y(3) DFY(1,2)=0.D0 DFY(1,3)=-B*Y(1) DFY(1,4)=0.D0 DFY(2,1)=A DFY(2,2)=-CM*Y(3) DFY(2,3)=-CM*Y(2) DFY(2,4)=0.D0 DFY(3,1)=A-B*Y(3) DFY(3,2)=-CM*Y(3) DFY(3,3)=-B*Y(1)-CM*Y(2) DFY(3,4)=C DFY(4,1)=B*Y(3) DFY(4,2)=0.D0 DFY(4,3)=B*Y(1) DFY(4,4)=-C RETURN END C SUBROUTINE FTIGE(NN,T,TH,F,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION F(NN),TH(150),U(150),V(150),W(150) DIMENSION ALPHA(150),BETA(150),STH(150),CTH(150) COMMON/NNNN/N,NNCOM,NSQ,NQUATR,DELTAS C ----- CALCUL DES TH(I) ET DES SIN ET COS ------------- DO 22 I=2,N THDIFF=TH(I)-TH(I-1) STH(I)=DSIN(THDIFF) 22 CTH(I)=DCOS(THDIFF) C -------- CALCUL DU COTE DROIT DU SYSTEME LINEAIRE ----- IF(T.GT.3.14159265358979324D0)THEN C --------- CASE T GREATER PI ------------ C ---------- I=1 ------------ TERM1=(-3.D0*TH(1)+TH(2))*NQUATR V(1)=TERM1 C -------- I=2,..,N-1 ----------- DO 32 I=2,N-1 TERM1=(TH(I-1)-2.D0*TH(I)+TH(I+1))*NQUATR 32 V(I)=TERM1 C ----------- I=N ------------- TERM1=(TH(N-1)-TH(N))*NQUATR V(N)=TERM1 ELSE C --------- CASE T LESS EQUAL PI ------------ FABS=1.5D0*DSIN(T)*DSIN(T) FX=-FABS FY= FABS C ---------- I=1 ------------ TERM1=(-3.D0*TH(1)+TH(2))*NQUATR TERM2=NSQ*(FY*DCOS(TH(1))-FX*DSIN(TH(1))) V(1)=TERM1+TERM2 C -------- I=2,..,N-1 ----------- DO 34 I=2,N-1 TERM1=(TH(I-1)-2.D0*TH(I)+TH(I+1))*NQUATR TERM2=NSQ*(FY*DCOS(TH(I))-FX*DSIN(TH(I))) 34 V(I)=TERM1+TERM2 C ----------- I=N ------------- TERM1=(TH(N-1)-TH(N))*NQUATR TERM2=NSQ*(FY*DCOS(TH(N))-FX*DSIN(TH(N))) V(N)=TERM1+TERM2 END IF C -------- COMPUTE PRODUCT DV=W ------------ W(1)=STH(2)*V(2) DO 43 I=2,N-1 43 W(I)=-STH(I)*V(I-1)+STH(I+1)*V(I+1) W(N)=-STH(N)*V(N-1) C -------- TERME 3 ----------------- DO 435 I=1,N 435 W(I)=W(I)+TH(N+I)*TH(N+I) C ------- SOLVE SYSTEM CW=W --------- ALPHA(1)=1.D0 DO 44 I=2,N ALPHA(I)=2.D0 44 BETA(I-1)=-CTH(I) ALPHA(N)=3.D0 DO 45 I=N-1,1,-1 Q=BETA(I)/ALPHA(I+1) W(I)=W(I)-W(I+1)*Q 45 ALPHA(I)=ALPHA(I)-BETA(I)*Q W(1)=W(1)/ALPHA(1) DO 46 I=2,N 46 W(I)=(W(I)-BETA(I-1)*W(I-1))/ALPHA(I) C -------- COMPUTE U=CV+DW --------- U(1)=V(1)-CTH(2)*V(2)+STH(2)*W(2) DO 47 I=2,N-1 47 U(I)=2.D0*V(I)-CTH(I)*V(I-1)-CTH(I+1)*V(I+1) & -STH(I)*W(I-1)+STH(I+1)*W(I+1) U(N)=3.D0*V(N)-CTH(N)*V(N-1)-STH(N)*W(N-1) C -------- PUT DERIVATIVES IN RIGHT PLACE ------------- DO 54 I=1,N F(I)=TH(N+I) 54 F(N+I)=U(I) RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT TIGE PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=80,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) COMMON/NNNN/NCOM,NNCOM,NSQ,NQUATR,DELTAS DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK) dimension true(nd),error(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) EXTERNAL FTIGE,SOLOUT c ------ FILE DE DONNEES ---------- write(8,3050) 3050 format(1x,'radau on the beam problem') C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=5 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP sumh=0.0D+0 summ=0.0D+0 C ---------- CONSTANTS -------------- N=40 NN=2*N NCOM=N NSQ=N*N NQUATR=NSQ*NSQ NNCOM=NN AN=N DELTAS=1.D0/AN C --- SET DEFAULT VALUES DO 12 I=1,20 WORK(I)=0.D0 12 IWORK(I)=0 c WORK(3)=0.1 c WORK(4)=0.3 c WORK(5)=0.99 c WORK(6)=2.0 C --------- INITIAL VALUES ------------- T=0.D0 DO 1 I=1,NN 1 Y(I)=0.D0 TEND=5.D0 H=1.D-3 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --------- SWITCH FOR OUTPUT -------- IOUT=0 C --- SECOND ORDER OPTION IJAC=0 c IWORK(9)=N c MLJAC=N MLJAC=NN C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM IMAS=0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) CALL RADAU5(NN,FTIGE,T,Y,TEND,H, & RTOL,ATOL,ITOL, & FTIGE,IJAC,MLJAC,MUJAC, & FTIGE,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) c CALL DTIME(TARRAY) true(1)=-0.005792366591294675 true(2)=-0.016952985507199259 true(3)=-0.027691033129713322 true(4)=-0.038008156558781729 true(5)=-0.047906168597422688 true(6)=-0.057387104352737008 true(7)=-0.066453273134522699 true(8)=-0.075107305819780661 true(9)=-0.083352197654124544 true(10)=-0.091191346546446469 true(11)=-0.098628587001297248 true(12)=-0.105668220037774708 true(13)=-0.112315039540924422 true(14)=-0.1 18574355272698475 true(15)=-0.124452012875526880 true(16)=-0.129954411326390999 true(17)=-0.135088518061004200 true(18)=-0.139861881919410397 true(19)=-0.144282644101482929 true(20)=-0.148359547246256976 true(21)=-0.152101942900106414 true(22)=-0.155519797806080921 true(23)=-0.158623699341992299 true(24)=-0.161424860370167541 true(25)=-0.163935123819275499 true(26)=-0.166166967344037066 true(27)=-0.168133508177817718 true(28)=-0.169848508060189926 true(29)=-0.171326378244038509 true(30)=-0.172582184746215274 true(31)=-0.173631653797526901 true(32)=-0.174491177383960691 true(33)=-0.175177818786287100 true(34)=-0.175709317871242317 true(35)=-0.176104096022807288 true(36)=-0.176381260717507812 true(37)=-0.176560609756417469 true(38)=-0.176662635226010517 true(39)=-0.176708527080694206 true(40)=-0.176720176107510191 true(41)= 0.037473626808570053 true(42)=0.109911788012810762 true(43)=0.179836047447039129 true(44)=0.247242730557127186 true(45)=0.312129382035491301 true(46)=0.374494737701689822 true(47)=0.434338607372647125 true(48)=0.491662035432760524 true(49)=0.546467785483476383 true(50)=0.598760970245279030 true(51)=0.648549361126755851 true(52)=0.695843516905088648 true(53)=0.740657266848912124 true(54)=0.783008174791347177 true(55)=0.822917665884869456 true(56)=0.860411030561688098 true(57)=0.895517550233742218 true(58)=0.928270826293034365 true(59)=0.958708933474210358 true(60)=0.986874782150222219 true(61)=1.012816579967983789 true(62)=1.036587736684594479 true(63)=1.058246826485315033 true(64)=1.077857811432700289 true(65)=1.095490221995530989 true(66)=1.111219164319120026 true(67)=1.125125269269998022 true(68)=1.137294526582397119 true(69)=1.147818025203744592 true(70)=1.156792131966898566 true(71)=1.164318845152484938 true(72)=1.170505992580311363 true(73)=1.175467424328008220 true(74)=1.179323003206967714 true(75)=1.182198586301326345 true(76)=1.184226111211404704 true(77)=1.185543909813440450 true(78)=1.186297084230907673 true(79)=1.186637618874913665 true(80)=1.186724615129383839 sum=0.0D+0 do 270 k=1,80 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/80.0D+0) sumh=sumh+sum summ=max(summ,sum) 20 continue it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh ID=0 c WRITE(6,9921)(Y(I),I=1,NN) 9921 FORMAT(1X,F22.16) WRITE (6,*) (IWORK(I),I=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (IWORK(I),I=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- IF(TARRAY(1).GT.1000.)STOP TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER(MX=8,MY=5,MACHS1=2,MACHS2=4) PARAMETER (ND=2*MX*MY,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) dimension true(nd),error(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) EXTERNAL FPLATE,JPLATE,JPLATS,JPLATSB,SOLOUT C -------- CONSTANTES ------------------ COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT NX=MX NY=MY NACHS1=MACHS1 NACHS2=MACHS2 NXM1=NX-1 NYM1=NY-1 NDEMI=NX*NY OMEGA=1000.D0 STIFFN=100.D0 WEIGHT=200.D0 DENOM=NX+1 DELX=2.D0/DENOM USH4=1.D0/(DELX**4) FAC=STIFFN*USH4 c ------ FILE DE DONNEES ---------- write(6,3500) 3500 format(1x,'radau results for plate problem') C -------- TWO OPTIONS ------- C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP sumh=0.0D+0 summ=0.0D+0 C --- DIMENSION OF THE SYSTEM N=ND C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 C --- INITIAL VALUES X=0.0D0 DO 1 I=1,N 1 Y(I)=0.D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL*1.0D-3 ITOL=0 C --- INITIAL STEP SIZE H=1.0D-2 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- SECOND ORDER OPTION AND BANDED IJAC=1 IWORK(9)=N/2 MLJAC=2*MX MUJAC=2*MX c MLJAC=N C --- ENDPOINT OF INTEGRATION XEND=7.D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL OF THE SUBROUTINE RADAU5 CALL RADAU5(N,FPLATE,X,Y,XEND,H, & RTOL,ATOL,ITOL, c & JPLATE,IJAC,MLJAC,MUJAC, & JPLATSB,IJAC,MLJAC,MUJAC, & FPLATE,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION DO 15 K=1,N 15 WRITE (6,*)Y(K) c CALL DTIME(TARRAY) true(1)=0.000490143813851336 true(2)=0.000980081485560611 true(3)=0.001462893811482190 true(4)=0.001915822464411935 true(5)=0.002285152533727002 true(6)=0.002461353376688549 true(7)=0.002254597413097122 true(8)=0.001438312591933600 true(9)=0.000849025149228402 true(10)=0.001697885005625757 true(11)=0.002535239886068847 true(12)=0.003323989552181772 true(13)=0.003977902193560667 true(14)=0.004320231736082990 true(15)=0.004025679955083897 true(16)=0.002643206356123840 true(17)=0.000980287627702671 true(18)=0.001960162971121222 true(19)=0.002925787622964379 true(20)=0.003831644928823870 true(21)=0.004570305067454005 true(22)=0.004922706753377098 true(23)=0.004509194826194244 true(24)=0.002876625183867201 true(25)=0.000849025149228402 true(26)= 0.001697885005625757 true(27)=0.002535239886068847 true(28)=0.003323989552181772 true(29)=0.003977902193560667 true(30)=0.004320231736082990 true(31)=0.004025679955083897 true(32)=0.002643206356123840 true(33)=0.000490143813851336 true(34)=0.000980081485560611 true(35)=0.001462893811482190 true(36)=0.001915822464411935 true(37)=0.002285152533727002 true(38)=0.002461353376688549 true(39)=0.002254597413097122 true(40)=0.001438312591933600 true(41)=- 0.001177590304545409 true(42)=-0.002409005827992214 true(43)=-0.003722140831656533 true(44)=-0.005078780056048207 true(45)=-0.006302661811097914 true(46)=-0.006973399942926759 true(47)=-0.006394575120415784 true(48)=-0.003960464551310118 true(49)=-0.002040148244040460 true(50)=-0.004174829877953482 true(51)=-0.006456510337516159 true(52)=-0.008832503276738242 true(53)=-0.011029624807177369 true(54)=-0.012352389570141255 true(55)=-0.011524177328690540 true(56)=-0.007253301886026949 true(57)=-0.002355180609090818 true(58)=-0.004818011655984428 true(59)=-0.007444281663313065 true(60)=-0.010157560112096413 true(61)=-0.012605323622195828 true(62)=-0.013946799885853517 true(63)=-0.012789150240831569 true(64)=-0.007920929102620235 true(65)=-0.002040148244040460 true(66)=-0.004174829877953482 true(67)=-0.006456510337516159 true(68)=-0.008832503276738242 true(69)=-0.011029624807177369 true(70)=-0.012352389570141255 true(71)=-0.011524177328690540 true(72)=-0.007253301886026949 true(73)=-0.001177590304545409 true(74)=-0.002409005827992214 true(75)=-0.003722140831656533 true(76)=-0.005078780056048207 true(77)=-0.006302661811097914 true(78)=-0.006973399942926759 true(79)=-0.006394575120415784 true(80)=-0.003960464551310118 sum=0.0D+0 do 270 k=1,80 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/80.0D+0) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh WRITE(6,*)TARRAY(1) WRITE (6,*)(ISTAT(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT WRITE (6,991) X,NR-1 991 FORMAT(1X,'X =',F9.5,' NSTEP =',I4) C WRITE(6,992)(Y(I),I=1,N,N/10) 992 FORMAT(10X,10F14.10) RETURN END C SUBROUTINE FPLATE (N, X, Y, F, RPAR, IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), F(N) COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT C -------- LA BOUCLE ------- DO 1 I=1,NX DO 1 J=1,NY K=I+NX*(J-1) C -------- DERIVEE DEUXIEME ---- F(K)=Y(K+NDEMI) C ------ POINT CENTRAL --- UC=16.D0*Y(K) IF(I.GT.1)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K-1) END IF IF(I.LT.NX)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K+1) END IF IF(J.GT.1)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K-NX) END IF IF(J.LT.NY)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K+NX) END IF IF(I.GT.1 .AND.J.GT.1 )UC=UC+2.D0*Y(K-NX-1) IF(I.LT.NX.AND.J.GT.1 )UC=UC+2.D0*Y(K-NX+1) IF(I.GT.1 .AND.J.LT.NY)UC=UC+2.D0*Y(K+NX-1) IF(I.LT.NX.AND.J.LT.NY)UC=UC+2.D0*Y(K+NX+1) IF(I.GT.2)UC=UC+Y(K-2) IF(I.LT.NXM1)UC=UC+Y(K+2) IF(J.GT.2)UC=UC+Y(K-2*NX) IF(J.LT.NYM1)UC=UC+Y(K+2*NX) IF(J.EQ.NACHS1.OR.J.EQ.NACHS2)THEN XI=I*DELX FORCE=EXP(-5.D0*(X-XI-2.D0)**2)+EXP(-5.D0*(X-XI-5.D0)**2) ELSE FORCE=0.D0 END IF F(K+NDEMI)=-OMEGA*Y(K+NDEMI)-FAC*UC+FORCE*WEIGHT 1 CONTINUE RETURN END c SUBROUTINE JPLATE(N,X,Y,DFY,LDFY,RPAR,IPAR) C -------- JACOBIAN FOR PLATE PROBLEM -------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) c ########### c DIMENSION AVP(80,80),WR(80),WI(80),ANGL(80) C ############## COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT C------ METTRE A ZERO ------- DO 1 I=1,N DO 1 J=1,N 1 DFY(I,J)=0.D0 C -------- LA BOUCLE ------- DO 2 I=1,NX DO 2 J=1,NY K=I+NX*(J-1) C -------- DERIVEE DEUXIEME ---- DFY(K,K+NDEMI)=1.D0 C ------ POINT CENTRAL --- DFY(K+NDEMI,K)=-FAC*16.D0 IF(I.GT.1)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K-1)=FAC*8.D0 END IF IF(I.LT.NX)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K+1)=FAC*8.D0 END IF IF(J.GT.1)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K-NX)=FAC*8.D0 END IF IF(J.LT.NY)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K+NX)=FAC*8.D0 END IF IF(I.GT.1 .AND.J.GT.1 )DFY(K+NDEMI,K-NX-1)=-FAC*2.D0 IF(I.LT.NX.AND.J.GT.1 )DFY(K+NDEMI,K-NX+1)=-FAC*2.D0 IF(I.GT.1 .AND.J.LT.NY)DFY(K+NDEMI,K+NX-1)=-FAC*2.D0 IF(I.LT.NX.AND.J.LT.NY)DFY(K+NDEMI,K+NX+1)=-FAC*2.D0 IF(I.GT.2)DFY(K+NDEMI,K-2)=-FAC IF(I.LT.NXM1)DFY(K+NDEMI,K+2)=-FAC IF(J.GT.2)DFY(K+NDEMI,K-2*NX)=-FAC IF(J.LT.NYM1)DFY(K+NDEMI,K+2*NX)=-FAC DFY(K+NDEMI,K+NDEMI)= -OMEGA 2 CONTINUE C ########## CALCUL VALEURS PROPRES ###### c DO 3 I=1,N c DO 3 J=1,N c 3 AVP(I,J)=DFY(I,J) c CALL VALPRP(80,N,AVP,WR,WI,IERR) c ANGLMX =0.D0 c PI=3.1415926535D0 c DO 4 I=1,N c ANGL(I)=180.D0-ABS(180.D0*ATAN2(WI(I),WR(I))/PI) c ANGLMX=MAX(ANGLMX,ANGL(I)) c WRITE(6,993)I,WR(I),WI(I),ANGL(I) c 4 CONTINUE c 993 FORMAT(' Eigenvalues -->',i4,2f12.5,' Angle=',f6.2) c WRITE(6,*)' ANGLEMAX=',ANGLMX C RETURN END c SUBROUTINE JPLATS(N,X,Y,DFY,LDFY,RPAR,IPAR) C -------- JACOBIAN FOR PLATE PROBLEM -------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT C------ METTRE A ZERO ------- DO 1 I=1,NDEMI DO 1 J=1,N 1 DFY(I,J)=0.D0 C -------- LA BOUCLE ------- DO 2 I=1,NX DO 2 J=1,NY K=I+NX*(J-1) C ------ POINT CENTRAL --- DFY(K,K)=-FAC*16.D0 IF(I.GT.1)THEN DFY(K,K)=DFY(K,K)-FAC DFY(K,K-1)=FAC*8.D0 END IF IF(I.LT.NX)THEN DFY(K,K)=DFY(K,K)-FAC DFY(K,K+1)=FAC*8.D0 END IF IF(J.GT.1)THEN DFY(K,K)=DFY(K,K)-FAC DFY(K,K-NX)=FAC*8.D0 END IF IF(J.LT.NY)THEN DFY(K,K)=DFY(K,K)-FAC DFY(K,K+NX)=FAC*8.D0 END IF IF(I.GT.1 .AND.J.GT.1 )DFY(K,K-NX-1)=-FAC*2.D0 IF(I.LT.NX.AND.J.GT.1 )DFY(K,K-NX+1)=-FAC*2.D0 IF(I.GT.1 .AND.J.LT.NY)DFY(K,K+NX-1)=-FAC*2.D0 IF(I.LT.NX.AND.J.LT.NY)DFY(K,K+NX+1)=-FAC*2.D0 IF(I.GT.2)DFY(K,K-2)=-FAC IF(I.LT.NXM1)DFY(K,K+2)=-FAC IF(J.GT.2)DFY(K,K-2*NX)=-FAC IF(J.LT.NYM1)DFY(K,K+2*NX)=-FAC DFY(K,K+NDEMI)= -OMEGA 2 CONTINUE RETURN END c SUBROUTINE JPLATSB(N,X,Y,DFY,LDFY,RPAR,IPAR) C -------- JACOBIAN FOR PLATE PROBLEM -------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT C------ METTRE A ZERO ------- DO 1 I=1,LDFY DO 1 J=1,N 1 DFY(I,J)=0.D0 MU=2*NX+1 FAC2=FAC*2.0D0 FAC8=FAC*8.0D0 FAC16=FAC*16.0D0 C -------- LA BOUCLE ------- DO 2 I=1,NX DO 2 J=1,NY K=I+NX*(J-1) C ------ POINT CENTRAL --- DFY(MU,K)=-FAC16 IF(I.GT.1)THEN DFY(MU,K)=DFY(MU,K)-FAC DFY(MU+1,K-1)=FAC8 END IF IF(I.LT.NX)THEN DFY(MU,K)=DFY(MU,K)-FAC DFY(MU-1,K+1)=FAC8 END IF IF(J.GT.1)THEN DFY(MU,K)=DFY(MU,K)-FAC DFY(MU+NX,K-NX)=FAC8 END IF IF(J.LT.NY)THEN DFY(MU,K)=DFY(MU,K)-FAC DFY(MU-NX,K+NX)=FAC8 END IF IF(I.GT.1 .AND.J.GT.1 )DFY(MU+NX+1,K-NX-1)=-FAC2 IF(I.LT.NX.AND.J.GT.1 )DFY(MU+NX-1,K-NX+1)=-FAC2 IF(I.GT.1 .AND.J.LT.NY)DFY(MU-NX+1,K+NX-1)=-FAC2 IF(I.LT.NX.AND.J.LT.NY)DFY(MU-NX-1,K+NX+1)=-FAC2 IF(I.GT.2)DFY(MU+2,K-2)=-FAC IF(I.LT.NXM1)DFY(MU-2,K+2)=-FAC IF(J.GT.2)DFY(MU+2*NX,K-2*NX)=-FAC IF(J.LT.NYM1)DFY(MU-2*NX,K+2*NX)=-FAC DFY(MU,K+NDEMI)= -OMEGA 2 CONTINUE RETURN END SUBROUTINE XPLATE(N,X,Y,FX,RPAR,IPAR) C -------- JACOBIAN FOR PLATE PROBLEM -------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),FX(N) COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT C -------- METTRE A ZERO -------- DO 1 I=1,N 1 FX(I)=0.D0 C -------- LA BOUCLE ------- DO 4 I=1,NX DO 4 J=1,NY K=I+NX*(J-1) IF(J.EQ.NACHS1.OR.J.EQ.NACHS2)THEN XI=I*DELX FX(K+NDEMI)=WEIGHT*( & -10.D0*(X-XI-2.D0)*EXP(-5.D0*(X-XI-2.D0)**2) & -10.D0*(X-XI-5.D0)*EXP(-5.D0*(X-XI-5.D0)**2)) END IF 4 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT CUSP PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=96,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK) dimension true(nd),error(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) COMMON/NERVES/NNERV COMMON/DIFFCOEF/DIFFUS EXTERNAL FCUSP,JCUSP,SOLOUT c ------ FILE DE DONNEES ---------- OPEN(8,FILE='res_rad5') write(8,2045) 2045 format(1x,'results for radau on cusp') REWIND 8 C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP sumh=0.0D+0 summ=0.0D+0 c -------- INITIAL CONSTANTES -------- NNERV=32 DIFFUS=1.D0*NNERV*NNERV/144.D0 N=3*NNERV C --- COMPUTE THE JACOBIAN SWITCH IJAC=0 C --- JACOBIAN IS A FULL MATRIX MLJAC=3 MUJAC=3 C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 C ---------- VAL INIT ------------ X=0.D0 XEND=1.1D0 ANERV=NNERV DEL=2.D0*3.14159265358979324D0/ANERV DO 25 INERV=1,NNERV Y(3*INERV-2)=0.D0 Y(3*INERV-1)=-2.D0*COS(INERV*DEL) 25 Y(3*INERV)=2.D0*SIN(INERV*DEL) C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO 10 I=1,20 10 WORK(I)=0.D0 DO 12 I=1,20 12 IWORK(I)=0 C --- ENDPOINT OF INTEGRATION XEND=1.1D0 it1=mclock() C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) CALL RADAU5(N,FCUSP,X,Y,XEND,H, & RTOL,ATOL,ITOL, & FCUSP,IJAC,MLJAC,MUJAC, & FCUSP,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION true(1)=-1.335038235173363825 true(2)=-0.141920661299976116 true(3)=2.189999851122752954 true(4)=-1.290165517136865728 true(5)=0.292210513241939329 true(6)=2.524498007953815718 true(7)=-1.206268463248866837 true(8)=0.702876002804259450 true(9)=2.603037671957832137 true(10)=-1.081173370722796200 true(11)=1.054547339698463687 true(12)=2.403900155664309457 true(13)=-0.922551477213655165 true(14)=1.326991956338080918 true(15)=2.009305096775369514 true(16)=-0.743049818521982534 true(17)=1.516881284521938182 true(18)=1.537256339189765457 true(19)=-0.555201077072863929 true(20)=1.632603197056899916 true(21)=1.077437487481636876 true(22)=-0.369158363066035655 true(23)=1.687674223256960258 true(24)=0.673204001909134905 true(25)=-0.192671593795137051 true(26)=1.695724385342398648 true(27)=0.333758479535656796 true(28)=-0.030615931836238017 true(29)=1.667262708083148298 true(30)=0.050978698243957534 true(31)=0.117513584875619235 true(32)=1.607508563419502269 true(33)=-0.190604748914752998 true(34)=0.259898961244445617 true(35)=1.514823442340568479 true(36)=-0.411323787318084589 true(37)=0.411809029672400558 true(38)=1.379804789392094260 true(39)=-0.638121744946811883 true(40)=0.590441346230457812 true(41)=1.185589061664514451 true(42)=-0.905945997151089418 true(43)=0.803741778414404449 true(44)=0.910756427168161885 true(45)=-1.251345457775128863 true(46)=1.037877442048335801 true(47)=0.545036626743780133 true(48)=-1.683821753687386274 true(49)=1.239043542405442416 true(50)=0.169981336507011738 true(51)=-2.112958754094543822 true(52)=1.406385681620871097 true(53)=-0.235380986562835855 true(54)=-2.450796096861493262 true(55)=1.524334200774267799 true(56)=-0.633461856049010260 true(57)=-2.576413161518578492 true(58)=1.588649099727842025 true(59)=-0.986582203794960052 true(60)=-2.442161394270367795 true(61)=1.606022353430074771 true(62)=-1.269240297385074114 true(63)=-2.104018859237057304 true(64)=1.588788794126354791 true(65)=-1.473056296837721718 true(66)=-1.670122571729852451 true(67)=1.549115780473624505 true(68)=-1.603417743271582846 true(69)=-1.233609811984700428 true(70)=1.495889929838369103 true(71)=-1.672805947347039028 true(72)=-0.844976238622025912 true(73)=1.434154221021214460 true(74)=-1.695067644864453216 true(75)=-0.518751841694241591 true(76)=1.365334914988092461 true(77)=-1.681659890115195530 true(78)=-0.249120054639391870 true(79)=1.286403800980685275 true(80)=-1.639285126097438457 true(81)=-0.019980596158691364 true(82)=1.184974025791693888 true(83)=-1.567910985925968939 true(84)=0.194039539947058944 true(85)=1.011140518164431143 true(86)=-1.455860565434940201 true(87)=0.436843623543739620 true(88)=-1.349821324547813729 true(89)=-1.223845158570813537 true(90)=0.809099908070360854 true(91)=-1.355008974443540401 true(92)=-0.926131110369012061 true(93)=1.232945832067604962 true(94)=-1.352261107347051711 true(95)=-0.559070645046367311 true(96)=1.716745798614099647 sum=0.0D+0 do 270 k=1,96 error(k)=dabs(y(k)-true(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/96.0D+0) sumh=sumh+sum summ=max(sum,summ) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 20 CONTINUE C --- PRINT STATISTICS WRITE(6,*)TARRAY(1),summ,sumh summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' write (6,*) y(1),y(n/2),y(n) WRITE (8,*)(IWORK(J),J=14,20) WRITE(8,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (8,91) (IWORK(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) write (8,*) y(1),y(n/2),y(n) C -------- NEW TOLERANCE --- IF (TARRAY(1).GT.500.) STOP TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N) COMMON /INTERN/XOUT COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL write (6,*) 'stat',NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F8.4,' Y =',2F18.7,' NSTEP =',I4) RETURN END SUBROUTINE FCUSP(N,T,Y,F,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),F(N) COMMON/NERVES/NNERV COMMON/DIFFCOEF/DIFFUS c----------- DO 25 INERV=1,NNERV X=Y(3*INERV-2) A=Y(3*INERV-1) B=Y(3*INERV) IF(INERV.EQ.1)THEN XRIGHT=Y(3*NNERV-2) ARIGHT=Y(3*NNERV-1) BRIGHT=Y(3*NNERV) ELSE XRIGHT=Y(3*INERV-5) ARIGHT=Y(3*INERV-4) BRIGHT=Y(3*INERV-3) END IF IF(INERV.EQ.NNERV)THEN XLEFT=Y(1) ALEFT=Y(2) BLEFT=Y(3) ELSE XLEFT=Y(3*INERV+1) ALEFT=Y(3*INERV+2) BLEFT=Y(3*INERV+3) END IF XDOT=-10000.D0*(B+X*(A+X*X)) U=(X-0.7D0)*(X-1.3D0) V=U/(U+0.1D0) ADOT=B+0.07D0*V BDOT=(1.D0*(1.D0-A*A)*B-A)-0.4D0*X+0.035D0*V F(3*INERV-2)=XDOT+DIFFUS*(XLEFT-2.D0*X+XRIGHT) F(3*INERV-1)=ADOT+DIFFUS*(ALEFT-2.D0*A+ARIGHT) F(3*INERV) =BDOT+DIFFUS*(BLEFT-2.D0*B+BRIGHT) 25 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT KS PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (BANDED JACOBIAN) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH,ND=2*NH-2) PARAMETER (IJAC=1,MLJAC=0,MUJAC=0,IMAS=0) PARAMETER (LJAC=MLJAC+MUJAC+1,LE=2*MLJAC+MUJAC+1) PARAMETER (LWORK=ND*(LJAC+3*LE+12)+20,LIWORK=3*ND+20) C --- DECLARATIONS DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),U(N) dimension true(nd),error(nd) COMMON/TRANS/QQ,UZERO REAL*4 TARRAY(2) EXTERNAL FKS,JKS,SOLOUT C --- DATA FOR THE PROBLEM QQ=0.025D0 c ------ FILE DE DONNEES ---------- c OPEN(8,FILE='res_rad5') c REWIND 8 write(6,3050) 3050 format(1x,'results for ks by radau') C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=7 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- OUTPUT DURING INTEGRATION --- IOUT=0 C --- INITIAL VALUES T=0.0D0 C --- INITIAL POSITIONS IN Y-SPACE ---- AN=N DO I=1,N DELX=1.D0/AN X=DELX*(I-1) U1=MIN(X-0.0D0,0.1D0-X) U2=20.D0*(X-0.2D0)*(0.3D0-X) U3=MIN(X-0.6D0,0.7D0-X) U4=MIN(X-0.9D0,1.0D0-X) U(I)=16.D0*MAX(0.D0,U1,U2,U3,U4) END DO C --- FOURIER TRANSFORM --- CALL REALFT(U,NH,+1) DO I=1,N U(I)=U(I)/AN END DO C --- INITIAL POSITION IN FOURIER MODES --- UZERO=U(1) DO I=1,ND Y(I)=U(I+2) END DO C --- WRITE DATA TO FILE PLUS CONTROL --- U(1)=UZERO U(2)=0.D0 DO I=3,N U(I)=Y(I-2) END DO CALL REALFT(U,NH,-1) DO I=1,N U(I)=2.D0*U(I) END DO C --- ENDPOINT OF INTEGRATION TEND=100.0D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 IWORK(I)=0 WORK(I)=0.D0 END DO C --- CALL OF THE SUBROUTINE RADAU5 c CALL DTIME(TARRAY) it1=mclock() CALL RADAU5(ND,FKS,T,Y,TEND,H, & RTOL,ATOL,ITOL, & JKS,IJAC,MLJAC,MUJAC, & JKS,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c CALL DTIME(TARRAY) true(1)=0.0387294818017374 true(2)=0.0011594920936697 true(3)=0.0161135569739015 true(4)=-0.0497378990194494 true(5)=0.0161388094706610 true(6)=-0.0790294186417376 true(7)=-0.0562945063709753 true(8)=-0.0194501629216551 true(9)=-0.0317390967305054 true(10)=-0.0256301653146683 true(11)=-0.0896908152487471 true(16)=0.0452570608746347 true(21)=0.0110762324013310 true(26)=0.0095676261122614 true(31)=0.1089375008451597 true(36)=-0.1201229482736044 true(41)=-0.1354295218545231 true(46)=-0.0146538030311935 true(51)=-0.0092223615970260 true(56)=-0.2760611800802050 true(61)=-0.0174686339809433 true(66)=0.1982033914963349 true(71)=-0.1768843940011891 true(76)=0.3090308065592430 true(81)=-0.0952958329765616 true(86)=-0.1063846133920493 true(91)=0.0390151203721951 true(96)=0.0451950368314852 true(101)=0.0477521089259518 true(151)=0.0112038077853544 true(201)=0.0003656188545739 true(251)=-0.0001443415867781 true(301)=-0.0000090034942150 true(351)=0.0000016423881092 true(401)=0.0000001490634819 true(451)=0.0000000334389501 true(501)=-0.0000000020970060 true(551)=-0.0000000017295829 true(601)=0.0000000000353463 true(651)=0.0000000000036143 true(701)=0.0000000000004182 true(751)=-0.0000000000000002 true(801)=-0.0000000000000070 true(851)=-0.0000000000000012 true(901)=0.0000000000000001 true(951)=0.0000000000000000 true(1001)=0.0000000000000000 summ=0.0D+0 do 270 k=1,10 error(k)=dabs(y(k)-true(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue do 271 k=11,100,5 error(k)=dabs(y(k)-true(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 271 continue do 272 k=101,1024,50 error(k)=dabs(y(k)-true(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 272 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh ID=0 WRITE (6,*)(IWORK(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (IWORK(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC IF (TARRAY(1).GT.1000.) STOP 30 CONTINUE STOP END SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT IF(MOD(NR,10).EQ.1) WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F7.4,' Y =',2E18.10,' NSTEP =',I4) RETURN END C *********************************************************** include 'fft.f' cc This module calculates the r.h.s for the KS equation cc cc d/dt u = - (dx^2+ dx^4)u - u u' cc cc in its form with fourier modes. cc I.e., the system cc cc dt u_n = (n^2 q^2 - n^4 q^4) u_n cc - (i q n/2) sum_{n_1+n_2=n} u_{n_1} u_{n_2} cc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine FKS(ND,T,Y,F,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH) DIMENSION y(ND),f(ND) DIMENSION U(N) COMMON/TRANS/QQ,UZERO c --- copy y to u U(1)=UZERO U(2)=0.D0 DO I=3,N U(I)=Y(I-2) END DO CALL REALFT(U,NH,-1) DO I=1,N U(I)=2.D0*U(I) END DO C --- SQUARE ------ DO I=1,N U(I)=U(I)**2/2.D0 END DO C --- TRANSFORM BACK --- CALL REALFT(U,NH,+1) AN=N DO I=1,N U(I)=U(I)/AN END DO do J=1,NH-1 DIAG=(QQ*J)**2*(1.D0-(QQ*J)**2) F(2*J-1)=DIAG*Y(2*J-1) +QQ*J*U(2*J+2) F(2*J )=DIAG*Y(2*J ) -QQ*J*U(2*J+1) enddo end c-------------------------------------------------------------- subroutine JKS(ND,x,y,dfy,ldfy,RPAR,IPAR) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION y(ND),DFY(LDFY,ND) COMMON/TRANS/QQ,UZERO do J=1,NH-1 DIAG=(QQ*J)**2*(1.D0-(QQ*J)**2) DFY(1,2*J-1)=DIAG DFY(1,2*J )=DIAG enddo end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 C * * * * * * * * * * * * * * * * * * * * * * * * * compile equation compile //venus/user/hairer/hailib/time cfeh driver_radau5 ~book2/appendix/radau5 ~book2/appendix/decsol equation time IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=15,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) dimension yend(nd),error(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) common/const/pi EXTERNAL FHIRES,JHIRES,SOLOUT c ------ FILE DE DONNEES ---------- write(6,3050) 3050 format(1x,'results on Hires problem') pi=4.0D+0*datan(1.0D+0) C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=1 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP sumh=0.0D+0 summ=0.0D+0 C --- DIMENSION OF THE SYSTEM N=15 C --- COMPUTE THE JACOBIAN IJAC=1 C --- JACOBIAN IS A FULL MATRIX MLJAC=N C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 C --- INITIAL VALUES X=0.0D0 do 171 ijk=1,15 y(ijk)=0.0D+0 171 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL*1.D-4 atol=rtol ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=1.0d-3 c CALL DTIME(TARRAY) it1=mclock() DO 20 I=1,1 C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) CALL RADAU5(N,FHIRES,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JHIRES,IJAC,MLJAC,MUJAC, & FHIRES,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c DO 15 K=1,8 c 15 WRITE (8,*)Y(K) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO XEND=XEND+100.D0 c CALL DTIME(TARRAY) yend(1)=-0.170799033d-1 yend(2)=-0.6660979d-2 yend(3)=+0.27531919d0 yend(4)=-0.39115732d0 yend(5)=-0.38851731d0 yend(6)=+0.27795920d0 yend(7)=+0.1114600281d0 yend(8)=+0.2979129627d-6 yend(9)=-0.31427403d-7 yend(10)=+0.70165883d-3 yend(11)=+0.85207538d-3 yend(12)=-0.77741454d-3 yend(13)=-0.77631967d-3 yend(14)=+0.78439426d-4 yend(15)=+0.252322784d-4 c sum=0.0D+0 do 270 k=1,15 error(k)=dabs(yend(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 continue it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) WRITE(6,*)TARRAY(1),summ,sumh/11.0d+0 summ=summ*atol sumh=sumh*atol/(11.0D+0) write(6,*) summ,sumh WRITE(6,*)NFE,NJE,NSTEP,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- WRITE (6,*)(ISTAT(J),J=14,20) WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) RETURN END C SUBROUTINE FHIRES (N, tn, Y, DY, RPAR, IPAR) INTEGER N DOUBLE PRECISION tn, Y(*), DY(*) c c implicit double precision (a-h,o-z) double precision y0, yN double precision kcns,pi,mi,dif,dx,amult,prod,rat common /const/pi double precision + c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc, + gamma,delta,sumer, + uin1,uin2,ud1,ud2,ud3,ud4,qud1,qud2,qud3,qud4 parameter (c=1.6d-8,cs=1d-9,cp=1d-8,r=25d3,rp=50d0, + lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4, + rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2, + gamma=40.67286402d-9,delta=17.7493332d0) uin1=0.5d0*sin(2.0d3*pi*tn) uin2=2.0d0*sin(2.0d4*pi*tn) ud1=+y(3)-y(5)-y(7)-uin2 ud2=-y(4)+y(6)-y(7)-uin2 ud3=+y(4)+y(5)+y(7)+uin2 ud4=-y(3)-y(6)+y(7)+uin2 qud1=gamma*(exp(delta*ud1)-1.0d0) qud2=gamma*(exp(delta*ud2)-1.0d0) qud3=gamma*(exp(delta*ud3)-1.0d0) qud4=gamma*(exp(delta*ud4)-1.0d0) dy(1)=(y(8)-0.5d0*y(10)+0.5d0*y(11)+y(14)-y(1)/r)/c dy(2)=(y(9)-0.5d0*y(12)+0.5d0*y(13)+y(15)-y(2)/r)/c dy(3)=(y(10)-qud1+qud4)/cs dy(4)=(-y(11)+qud2-qud3)/cs dy(5)=(y(12)+qud1-qud3)/cs dy(6)=(-y(13)-qud2+qud4)/cs dy(7)=(-y(7)/rp+qud1+qud2-qud3-qud4)/cp dy(8)=-y(1)/lh dy(9)=-y(2)/lh dy(10)=(0.5d0*y(1)-y(3)-rg2*y(10))/ls2 dy(11)=(-0.5d0*y(1)+y(4)-rg3*y(11))/ls3 dy(12)=(0.5d0*y(2)-y(5)-rg2*y(12))/ls2 dy(13)=(-0.5d0*y(2)+y(6)-rg3*y(13))/ls3 dy(14)=(-y(1)+uin1-(ri+rg1)*y(14))/ls1 dy(15)=(-y(2)-(rc+rg1)*y(15))/ls1 return end SUBROUTINE JHIRES(N,tn,Y,pd,LDFY,RPAR,IPAR) C -------- JACOBIAN FOR HIRES PROBLEM -------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),pd(LDFY,N) C------ METTRE A ZERO ------- c c implicit double precision (a-h,o-z) double precision kcns,pi,mi,dif,dx,amult,prod,rat,a1,a2,a3,u common /const/pi double precision + c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc, + gamma,delta,sumer, + uin2,ud1,ud2,ud3,ud4,qpud1,qpud2,qpud3,qpud4 parameter (c=1.6d-8,cs=1d-9,cp=1d-8,r=25d3,rp=50d0, + lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4, + rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2, + gamma=40.67286402d-9,delta=17.7493332d0) do 177 i=1,n do 177 j=1,n pd(i,j)=0.0D+0 177 continue uin2=2.0d0*sin(2.0d4*pi*tn) ud1=+y(3)-y(5)-y(7)-uin2 ud2=-y(4)+y(6)-y(7)-uin2 ud3=+y(4)+y(5)+y(7)+uin2 ud4=-y(3)-y(6)+y(7)+uin2 qpud1=gamma*delta*exp(delta*ud1) qpud2=gamma*delta*exp(delta*ud2) qpud3=gamma*delta*exp(delta*ud3) qpud4=gamma*delta*exp(delta*ud4) pd(1,1)=-1.0d0/(c*r) pd(1,8)=-1.0d0/c pd(1,10)=-0.5d0/c pd(1,11)=-pd(1,10) pd(1,14)=pd(1,8) pd(2,2)=pd(1,1) pd(2,9)=pd(1,8) pd(2,12)=pd(1,10) pd(2,13)=pd(1,11) pd(2,15)=pd(1,14) pd(3,3)=(-qpud1-qpud4)/cs pd(3,5)=qpud1/cs pd(3,6)=-qpud4/cs pd(3,7)=(qpud1+qpud4)/cs pd(3,10)=1.0d0/cs pd(4,4)=(-qpud2-qpud3)/cs pd(4,5)=-qpud3/cs pd(4,6)=qpud2/cs pd(4,7)=(-qpud2-qpud3)/cs pd(4,11)=-1.0d0/cs pd(5,3)=qpud1/cs pd(5,4)=-qpud3/cs pd(5,5)=(-qpud1-qpud3)/cs pd(5,7)=(-qpud1-qpud3)/cs pd(5,12)=1.0d0/cs pd(6,3)=-qpud4/cs pd(6,4)=qpud2/cs pd(6,6)=(-qpud2-qpud4)/cs pd(6,7)=(qpud2+qpud4)/cs pd(6,13)=-1.0d0/cs pd(7,3)=(qpud1+qpud4)/cp pd(7,4)=(-qpud2-qpud3)/cp pd(7,5)=(-qpud1-qpud3)/cp pd(7,6)=(qpud2+qpud4)/cp pd(7,7)=(-qpud1-qpud2-qpud3-qpud4-1.0d0/rp)/cp pd(8,1)=-1.0d0/lh pd(9,2)=pd(8,1) pd(10,1)=0.5d0/ls2 pd(10,3)=-1.0d0/ls2 pd(10,10)=-rg2/ls2 pd(11,1)=-0.5d0/ls3 pd(11,4)=1.0d0/ls3 pd(11,11)=-rg3/ls3 pd(12,2)=pd(10,1) pd(12,5)=pd(10,3) pd(12,12)=pd(10,10) pd(13,2)=pd(11,1) pd(13,6)=pd(11,4) pd(13,13)=pd(11,11) pd(14,1)=-1.0d0/ls1 pd(14,14)=-(ri+rg1)/ls1 pd(15,2)=pd(14,1) pd(15,15)=-(rc+rg1)/ls1 return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT VDPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=6,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) + ,true(nd),error(nd) common/add/sumh,rtol,atol C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) EXTERNAL FVANDER,JVANDER,SOLOUT c ------ FILE DE DONNEES ---------- write(6,3500) 3500 format(1x,'radau5 on B5') RPAR=1.D-6 C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+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 EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 iout=1 C --- INITIAL VALUES X=0.0D0 do 291 i=1,6 y(i)=1.0D+0 291 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=20.0D0 it1=mclock() C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) 220 continue CALL RADAU5(N,FVANDER,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JVANDER,IJAC,MLJAC,MUJAC, & FVANDER,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c WRITE (8,*) Y(1) c WRITE (8,*) Y(2) write(6,*) y(1),y(2),x,h C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO sum=0.0D+0 param=1000.0D+0 true(1)=dexp(-10.0*x)*(cos(param*x)+sin(param*x)) true(2)=dexp(-10.0*x)*(cos(param*x)-sin(param*x)) true(3)=dexp(-4.0D+0*x) true(4)=dexp(-x) true(5)=dexp(-0.5D+0*x) true(6)=dexp(-0.1D+0*x) do 270 k=1,6 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue summ=dsqrt(sum/dble(nd)) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) WRITE(6,*)TARRAY(1),summ,sumh summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh WRITE (6,*)(ISTAT(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT common/add/sumh,rtol,atol dimension error(6),true(6) c WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) param=1000.0D+0 true(1)=dexp(-10.0*x)*(cos(param*x)+sin(param*x)) true(2)=dexp(-10.0*x)*(cos(param*x)-sin(param*x)) true(3)=dexp(-4.0D+0*x) true(4)=dexp(-x) true(5)=dexp(-0.5D+0*x) true(6)=dexp(-0.1D+0*x) sum=0.0D+0 do 270 k=1,6 error(k)=true(k)-y(k) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/6.0D+0) if(sum.gt.sumh) sumh=sum c write(6,*) sum,sumh return END SUBROUTINE FVANDER(N,X,Y,F,EPS,IPAR) C --- RIGHT-HAND SIDE OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),F(N) f(1)=-10.0D+0*y(1)+1000.0D+0*y(2) f(2)=-1000.0D+0*y(1)-10.0D+0*y(2) f(3)=-4.0D+0*y(3) f(4)=-y(4) f(5)=-0.5D+0*y(5) f(6)=-0.1D+0*y(6) RETURN END C SUBROUTINE JVANDER(N,X,Y,DFY,LDFY,EPS,IPAR) C --- JACOBIAN OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) DFY(1,1)=-10.0D0 DFY(1,2)=1000.D0 DFY(2,1)=-1000.0D+0 DFY(2,2)=-10.0D+0 dfy(3,3)=-4.0D+0 dfy(4,4)=-1.0D+0 dfy(5,5)=-0.5D+0 dfy(6,6)=-0.1D+0 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT KS PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (BANDED JACOBIAN) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH,ND=2*NH-2) PARAMETER (IJAC=1,MLJAC=0,MUJAC=0,IMAS=0) PARAMETER (LJAC=MLJAC+MUJAC+1,LE=2*MLJAC+MUJAC+1) PARAMETER (LWORK=ND*(LJAC+3*LE+12)+20,LIWORK=3*ND+20) C --- DECLARATIONS DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),U(N) dimension true(nd),error(nd) COMMON/TRANS/QQ,UZERO REAL*4 TARRAY(2) EXTERNAL FKS,JKS,SOLOUT C --- DATA FOR THE PROBLEM QQ=0.025D0 c ------ FILE DE DONNEES ---------- OPEN(6,FILE='radks.res') REWIND 6 write(6,3050) 3050 format(1x,'results for ks by radau') C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=2 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- OUTPUT DURING INTEGRATION --- IOUT=0 C --- INITIAL VALUES T=0.0D0 C --- INITIAL POSITIONS IN Y-SPACE ---- AN=N DO I=1,N DELX=1.D0/AN X=DELX*(I-1) U1=MIN(X-0.0D0,0.1D0-X) U2=20.D0*(X-0.2D0)*(0.3D0-X) U3=MIN(X-0.6D0,0.7D0-X) U4=MIN(X-0.9D0,1.0D0-X) U(I)=16.D0*MAX(0.D0,U1,U2,U3,U4) END DO C --- FOURIER TRANSFORM --- CALL REALFT(U,NH,+1) DO I=1,N U(I)=U(I)/AN END DO C --- INITIAL POSITION IN FOURIER MODES --- UZERO=U(1) DO I=1,ND Y(I)=U(I+2) END DO C --- WRITE DATA TO FILE PLUS CONTROL --- U(1)=UZERO U(2)=0.D0 DO I=3,N U(I)=Y(I-2) END DO CALL REALFT(U,NH,-1) DO I=1,N U(I)=2.D0*U(I) END DO C --- ENDPOINT OF INTEGRATION TEND=100.0D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 IWORK(I)=0 WORK(I)=0.D0 END DO C --- CALL OF THE SUBROUTINE RADAU5 c CALL DTIME(TARRAY) it1=mclock() CALL RADAU5(ND,FKS,T,Y,TEND,H, & RTOL,ATOL,ITOL, & JKS,IJAC,MLJAC,MUJAC, & JKS,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c CALL DTIME(TARRAY) true(1)=0.0387294818017374 true(2)=0.0011594920936697 true(3)=0.0161135569739015 true(4)=-0.0497378990194494 true(5)=0.0161388094706610 true(6)=-0.0790294186417376 true(7)=-0.0562945063709753 true(8)=-0.0194501629216551 true(9)=-0.0317390967305054 true(10)=-0.0256301653146683 true(11)=-0.0896908152487471 true(16)=0.0452570608746347 true(21)=0.0110762324013310 true(26)=0.0095676261122614 true(31)=0.1089375008451597 true(36)=-0.1201229482736044 true(41)=-0.1354295218545231 true(46)=-0.0146538030311935 true(51)=-0.0092223615970260 true(56)=-0.2760611800802050 true(61)=-0.0174686339809433 true(66)=0.1982033914963349 true(71)=-0.1768843940011891 true(76)=0.3090308065592430 true(81)=-0.0952958329765616 true(86)=-0.1063846133920493 true(91)=0.0390151203721951 true(96)=0.0451950368314852 true(101)=0.0477521089259518 true(151)=0.0112038077853544 true(201)=0.0003656188545739 true(251)=-0.0001443415867781 true(301)=-0.0000090034942150 true(351)=0.0000016423881092 true(401)=0.0000001490634819 true(451)=0.0000000334389501 true(501)=-0.0000000020970060 true(551)=-0.0000000017295829 true(601)=0.0000000000353463 true(651)=0.0000000000036143 true(701)=0.0000000000004182 true(751)=-0.0000000000000002 true(801)=-0.0000000000000070 true(851)=-0.0000000000000012 true(901)=0.0000000000000001 true(951)=0.0000000000000000 true(1001)=0.0000000000000000 summ=0.0D+0 do 270 k=1,10 error(k)=dabs(y(k)-true(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue do 271 k=11,100,5 error(k)=dabs(y(k)-true(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 271 continue do 272 k=101,1024,50 error(k)=dabs(y(k)-true(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 272 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh us=log10(summ) ut=log10(sumh) write(6,*) us,ut ID=0 WRITE (6,*)(IWORK(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (IWORK(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC IF (TARRAY(1).GT.1000.) STOP 30 CONTINUE STOP END SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT IF(MOD(NR,10).EQ.1) WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F7.4,' Y =',2E18.10,' NSTEP =',I4) RETURN END C *********************************************************** include 'fft.f' cc This module calculates the r.h.s for the KS equation cc cc d/dt u = - (dx^2+ dx^4)u - u u' cc cc in its form with fourier modes. cc I.e., the system cc cc dt u_n = (n^2 q^2 - n^4 q^4) u_n cc - (i q n/2) sum_{n_1+n_2=n} u_{n_1} u_{n_2} cc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine FKS(ND,T,Y,F,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH) DIMENSION y(ND),f(ND) DIMENSION U(N) COMMON/TRANS/QQ,UZERO c --- copy y to u U(1)=UZERO U(2)=0.D0 DO I=3,N U(I)=Y(I-2) END DO CALL REALFT(U,NH,-1) DO I=1,N U(I)=2.D0*U(I) END DO C --- SQUARE ------ DO I=1,N U(I)=U(I)**2/2.D0 END DO C --- TRANSFORM BACK --- CALL REALFT(U,NH,+1) AN=N DO I=1,N U(I)=U(I)/AN END DO do J=1,NH-1 DIAG=(QQ*J)**2*(1.D0-(QQ*J)**2) F(2*J-1)=DIAG*Y(2*J-1) +QQ*J*U(2*J+2) F(2*J )=DIAG*Y(2*J ) -QQ*J*U(2*J+1) enddo end c-------------------------------------------------------------- subroutine JKS(ND,x,y,dfy,ldfy,RPAR,IPAR) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION y(ND),DFY(LDFY,ND) COMMON/TRANS/QQ,UZERO do J=1,NH-1 DIAG=(QQ*J)**2*(1.D0-(QQ*J)**2) DFY(1,2*J-1)=DIAG DFY(1,2*J )=DIAG enddo end C SUBROUTINE FTIGE(NN,T,TH,F,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION F(NN),TH(150),U(150),V(150),W(150) DIMENSION ALPHA(150),BETA(150),STH(150),CTH(150) COMMON/NNNN/N,NNCOM,NSQ,NQUATR,DELTAS C ----- CALCUL DES TH(I) ET DES SIN ET COS ------------- DO 22 I=2,N THDIFF=TH(I)-TH(I-1) STH(I)=DSIN(THDIFF) 22 CTH(I)=DCOS(THDIFF) C -------- CALCUL DU COTE DROIT DU SYSTEME LINEAIRE ----- IF(T.GT.3.14159265358979324D0)THEN C --------- CASE T GREATER PI ------------ C ---------- I=1 ------------ TERM1=(-3.D0*TH(1)+TH(2))*NQUATR V(1)=TERM1 C -------- I=2,..,N-1 ----------- DO 32 I=2,N-1 TERM1=(TH(I-1)-2.D0*TH(I)+TH(I+1))*NQUATR 32 V(I)=TERM1 C ----------- I=N ------------- TERM1=(TH(N-1)-TH(N))*NQUATR V(N)=TERM1 ELSE C --------- CASE T LESS EQUAL PI ------------ FABS=1.5D0*DSIN(T)*DSIN(T) FX=-FABS FY= FABS C ---------- I=1 ------------ TERM1=(-3.D0*TH(1)+TH(2))*NQUATR TERM2=NSQ*(FY*DCOS(TH(1))-FX*DSIN(TH(1))) V(1)=TERM1+TERM2 C -------- I=2,..,N-1 ----------- DO 34 I=2,N-1 TERM1=(TH(I-1)-2.D0*TH(I)+TH(I+1))*NQUATR TERM2=NSQ*(FY*DCOS(TH(I))-FX*DSIN(TH(I))) 34 V(I)=TERM1+TERM2 C ----------- I=N ------------- TERM1=(TH(N-1)-TH(N))*NQUATR TERM2=NSQ*(FY*DCOS(TH(N))-FX*DSIN(TH(N))) V(N)=TERM1+TERM2 END IF C -------- COMPUTE PRODUCT DV=W ------------ W(1)=STH(2)*V(2) DO 43 I=2,N-1 43 W(I)=-STH(I)*V(I-1)+STH(I+1)*V(I+1) W(N)=-STH(N)*V(N-1) C -------- TERME 3 ----------------- DO 435 I=1,N 435 W(I)=W(I)+TH(N+I)*TH(N+I) C ------- SOLVE SYSTEM CW=W --------- ALPHA(1)=1.D0 DO 44 I=2,N ALPHA(I)=2.D0 44 BETA(I-1)=-CTH(I) ALPHA(N)=3.D0 DO 45 I=N-1,1,-1 Q=BETA(I)/ALPHA(I+1) W(I)=W(I)-W(I+1)*Q 45 ALPHA(I)=ALPHA(I)-BETA(I)*Q W(1)=W(1)/ALPHA(1) DO 46 I=2,N 46 W(I)=(W(I)-BETA(I-1)*W(I-1))/ALPHA(I) C -------- COMPUTE U=CV+DW --------- U(1)=V(1)-CTH(2)*V(2)+STH(2)*W(2) DO 47 I=2,N-1 47 U(I)=2.D0*V(I)-CTH(I)*V(I-1)-CTH(I+1)*V(I+1) & -STH(I)*W(I-1)+STH(I+1)*W(I+1) U(N)=3.D0*V(N)-CTH(N)*V(N-1)-STH(N)*W(N-1) C -------- PUT DERIVATIVES IN RIGHT PLACE ------------- DO 54 I=1,N F(I)=TH(N+I) 54 F(N+I)=U(I) RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT TIGE PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=80,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) COMMON/NNNN/NCOM,NNCOM,NSQ,NQUATR,DELTAS DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK) dimension true(nd),error(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) EXTERNAL FTIGE,SOLOUT c ------ FILE DE DONNEES ---------- open(6,file='radaubeam.res') write(6,3050) 3050 format(1x,'radau on the beam problem') C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP sumh=0.0D+0 summ=0.0D+0 C ---------- CONSTANTS -------------- N=40 NN=2*N NCOM=N NSQ=N*N NQUATR=NSQ*NSQ NNCOM=NN AN=N DELTAS=1.D0/AN C --- SET DEFAULT VALUES DO 12 I=1,20 WORK(I)=0.D0 12 IWORK(I)=0 c WORK(3)=0.1 c WORK(4)=0.3 c WORK(5)=0.99 c WORK(6)=2.0 C --------- INITIAL VALUES ------------- T=0.D0 DO 1 I=1,NN 1 Y(I)=0.D0 TEND=5.D0 H=1.D-3 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --------- SWITCH FOR OUTPUT -------- IOUT=0 C --- SECOND ORDER OPTION IJAC=0 c IWORK(9)=N c MLJAC=N MLJAC=NN C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM IMAS=0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) CALL RADAU5(NN,FTIGE,T,Y,TEND,H, & RTOL,ATOL,ITOL, & FTIGE,IJAC,MLJAC,MUJAC, & FTIGE,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) c CALL DTIME(TARRAY) true(1)=-0.005792366591294675 true(2)=-0.016952985507199259 true(3)=-0.027691033129713322 true(4)=-0.038008156558781729 true(5)=-0.047906168597422688 true(6)=-0.057387104352737008 true(7)=-0.066453273134522699 true(8)=-0.075107305819780661 true(9)=-0.083352197654124544 true(10)=-0.091191346546446469 true(11)=-0.098628587001297248 true(12)=-0.105668220037774708 true(13)=-0.112315039540924422 true(14)=-0.1 18574355272698475 true(15)=-0.124452012875526880 true(16)=-0.129954411326390999 true(17)=-0.135088518061004200 true(18)=-0.139861881919410397 true(19)=-0.144282644101482929 true(20)=-0.148359547246256976 true(21)=-0.152101942900106414 true(22)=-0.155519797806080921 true(23)=-0.158623699341992299 true(24)=-0.161424860370167541 true(25)=-0.163935123819275499 true(26)=-0.166166967344037066 true(27)=-0.168133508177817718 true(28)=-0.169848508060189926 true(29)=-0.171326378244038509 true(30)=-0.172582184746215274 true(31)=-0.173631653797526901 true(32)=-0.174491177383960691 true(33)=-0.175177818786287100 true(34)=-0.175709317871242317 true(35)=-0.176104096022807288 true(36)=-0.176381260717507812 true(37)=-0.176560609756417469 true(38)=-0.176662635226010517 true(39)=-0.176708527080694206 true(40)=-0.176720176107510191 true(41)= 0.037473626808570053 true(42)=0.109911788012810762 true(43)=0.179836047447039129 true(44)=0.247242730557127186 true(45)=0.312129382035491301 true(46)=0.374494737701689822 true(47)=0.434338607372647125 true(48)=0.491662035432760524 true(49)=0.546467785483476383 true(50)=0.598760970245279030 true(51)=0.648549361126755851 true(52)=0.695843516905088648 true(53)=0.740657266848912124 true(54)=0.783008174791347177 true(55)=0.822917665884869456 true(56)=0.860411030561688098 true(57)=0.895517550233742218 true(58)=0.928270826293034365 true(59)=0.958708933474210358 true(60)=0.986874782150222219 true(61)=1.012816579967983789 true(62)=1.036587736684594479 true(63)=1.058246826485315033 true(64)=1.077857811432700289 true(65)=1.095490221995530989 true(66)=1.111219164319120026 true(67)=1.125125269269998022 true(68)=1.137294526582397119 true(69)=1.147818025203744592 true(70)=1.156792131966898566 true(71)=1.164318845152484938 true(72)=1.170505992580311363 true(73)=1.175467424328008220 true(74)=1.179323003206967714 true(75)=1.182198586301326345 true(76)=1.184226111211404704 true(77)=1.185543909813440450 true(78)=1.186297084230907673 true(79)=1.186637618874913665 true(80)=1.186724615129383839 sum=0.0D+0 do 270 k=1,80 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/80.0D+0) sumh=sumh+sum summ=max(summ,sum) 20 continue it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh ID=0 c WRITE(6,9921)(Y(I),I=1,NN) 9921 FORMAT(1X,F22.16) WRITE (6,*) (IWORK(I),I=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (IWORK(I),I=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- IF(TARRAY(1).GT.1000.)STOP TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT VDPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=49999,LIWORK=3*ND+20) parameter(ijac=0,mljac=1,mujac=1,imas=0) parameter(ljac=mljac+mujac+1,le=2*mljac+mujac+1) parameter(lwork=nd*(ljac+3*le+12)+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) + ,true(nd),error(nd) common/add/sumh,rtol,atol C -------- END PARAMETER COMMON /const/pi,dx,ntop REAL*4 TARRAY(2) EXTERNAL FVANDER,JVANDER,SOLOUT open(6,file='radauhard.bigres250') c ------ FILE DE DONNEES ---------- write(6,3500) 3500 format(1x,'radau5 on hard parabolic') RPAR=1.D-6 C --- LOOP FOR DIFFERENT TOLERANCES ntop=250 n=nd write(6,*) ntop,nd pi=4.0*DATAN(1.0D+0) c write(6,*) ntop c c dx=1.0d0/(n+1.0d0) NTOLMN=2 NTOLMX=10 NTOLDF=2 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM C --- COMPUTE THE JACOBIAN ANALYTICALLY C --- JACOBIAN IS A FULL MATRIX C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 iout=1 C --- INITIAL VALUES X=0.0D0 do 201 j=1,n sumer=0.0D+0 do 21 jj=1,ntop sumer=sumer+dsin(pi*jj*j*dx) 21 continue y(j)=sumer 201 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=2.0D0 it1=mclock() C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) 220 continue CALL RADAU5(N,FVANDER,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JVANDER,IJAC,MLJAC,MUJAC, & FVANDER,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c WRITE (8,*) Y(1) c WRITE (8,*) Y(2) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO sum=0.0D+0 c tn=xend do 31 j=1,n true(j)=0.0D+0 do 33 i=1,ntop true(j)=true(j)+dexp(-tn*pi*pi*i*i)*dsin(dx*pi*i*j) 33 continue 31 continue c c do 270 k=1,49999 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue summ=dsqrt(sum/dble(nd)) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh us=log10(summ) ut=log10(sumh) write(6,*) us,ut WRITE (6,*)(ISTAT(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT common/add/sumh,rtol,atol dimension error(49999),true(49999) COMMON /const/pi,dx,ntop c WRITE (6,99) X,Y(1),Y(2),NR-1 tn=x do 31 j=1,n sumer=0.0D+0 do 33 i=1,ntop sumer=sumer+dexp(-tn*pi*pi*i*i)*dsin(dx*pi*i*j) 33 continue true(j)=sumer 31 continue 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) sum=0.0D+0 do 270 k=1,49999 error(k)=true(k)-y(k) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/49999.0D+0) if(sum.gt.sumh) sumh=sum c write(6,*) sum,sumh return END SUBROUTINE FVANDER(N,X,Y,F,EPS,IPAR) C --- RIGHT-HAND SIDE OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),F(N) COMMON /const/pi,dx,ntop tn=x amult=2.0d0*dx prod=dx*dx c SUMER=0.0D+0 do 103 i=1,ntop sumer=sumer+dexp(-pi*pi*tn*i*i)*dsin(pi*dx*i) 103 continue F(1)=(y(2)-2.0d0*y(1))/(dx*dx)+y(1)**2-0.5d+0*y(1)*sumer- + 0.5d+0*sumer**2 do 101 i=2,n-1 sumer=0.0d+0 do 104 ii=1,ntop sumer=sumer+dexp(-pi*pi*tn*ii*ii)*dsin(pi*i*ii*dx) 104 continue f(i)=(y(i-1)-2.0d0*y(i)+y(i+1))/prod + +y(I)**2-0.5D+0*y(i)*sumer-0.5D+0*sumer**2 101 continue sumer=0.0D+0 do 105 II=1,ntop sumer=sumer+dexp(-pi*pi*tn*ii*ii)*dsin(pi*n*dx*ii) 105 continue F(n)=(y(n-1)-2.0d0*y(n))/(dx*dx) + +y(n)**2-0.5D+0*y(n)*sumer-0.5D+0*sumer**2 return end C SUBROUTINE JVANDER(N,X,Y,DFY,LDFY,EPS,IPAR) C --- JACOBIAN OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) COMMON /const/pi,dx,ntop tn=x prod=dx*dx amult=2.0d0*dx a1=1.0d0/prod a2=-2.0d0/prod do 141 j=2,n dfy(1,j)=a1 141 continue u=a2 do 142 j=1,n sumer=0.0D+0 do 143 k=1,ntop sumer=sumer+dexp(-pi*pi*tn*k*k)*dsin(pi*j*dx*k) 143 continue dfy(2,j)=u+2.0D+0*y(j)-0.5D+0*sumer 142 continue do 144 j=1,n-1 dfy(3,j)=a1 144 continue return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT VDPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=49999,LIWORK=3*ND+20) parameter(ijac=0,mljac=1,mujac=1,imas=0) parameter(ljac=mljac+mujac+1,le=2*mljac+mujac+1) parameter(lwork=nd*(ljac+3*le+12)+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) + ,true(nd),error(nd) common/add/sumh,rtol,atol C -------- END PARAMETER COMMON /const/pi,dx,dif,ntop REAL*4 TARRAY(2) EXTERNAL FVANDER,JVANDER,SOLOUT open(6,file='radauflood.hardres') c ------ FILE DE DONNEES ---------- write(6,3500) 3500 format(1x,'radau5 on flood problem') RPAR=1.D-6 C --- LOOP FOR DIFFERENT TOLERANCES ntop=20 dif=2.0D-5 dif=4.0D+0/500000.0D+0 n=nd write(6,*) dif,nd pi=4.0*DATAN(1.0D+0) dx=4.0d0/(n+1.0d0) NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM C --- COMPUTE THE JACOBIAN ANALYTICALLY C --- JACOBIAN IS A FULL MATRIX C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 iout=1 C --- INITIAL VALUES X=0.0D0 do 201 j=1,n y(j)=dexp(-(-2.0D+0+j*dx)**2) 201 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=4.0D0 it1=mclock() C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) 220 continue CALL RADAU5(N,FVANDER,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JVANDER,IJAC,MLJAC,MUJAC, & FVANDER,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c WRITE (8,*) Y(1) c WRITE (8,*) Y(2) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO sum=0.0D+0 c tn=xend do 31 j=1,n true(j)= + (1.0d0/dsqrt(1.0d0+4.0d0*dif*tn))* & dexp(-((-2.0d0+j*dx-tn)**2)/(1.0d0+4.0d0*dif*tn)) 31 continue c c do 270 k=1,49999 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue summ=dsqrt(sum/dble(nd)) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh us=log10(summ) ut=log10(sumh) write(6,*) us,ut WRITE (6,*)(ISTAT(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT common/add/sumh,rtol,atol dimension error(49999),true(49999) COMMON /const/pi,dx,dif,ntop c WRITE (6,99) X,Y(1),Y(2),NR-1 tn=x do 31 j=1,n true(j)=(1.0d0/dsqrt(1.0d0+4.0d0*dif*tn))* & dexp(-((-2.0d0+j*dx-tn)**2)/(1.0d0+4.0d0*dif*tn)) 31 continue 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) sum=0.0D+0 do 270 k=1,49999 error(k)=true(k)-y(k) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/49999.0D+0) if(sum.gt.sumh) sumh=sum c write(6,*) sum,sumh return END SUBROUTINE FVANDER(N,X,Y,F,EPS,IPAR) C --- RIGHT-HAND SIDE OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),F(N) COMMON /const/pi,dx,dif,ntop tn=x amult=2.0d0*dx prod=dx*dx c y0=(1.0d0/dsqrt(1.d0+4.0d0*dif*tn))*dexp(-((-2.0d0-tn)**2)/ & (1.0d0+4.0d0*dif*tn)) yN=(1.0d0/dsqrt(1.d0+4.d0*dif*tn))*dexp(-((2.d0-tn)**2)/ & (1.0d0+4.0d0*dif*tn)) f(1)=dif*(y0-2.d0*y(1)+y(2))/(dx*dx)-(y(2)-y0)/(2.d0*dx) rat=dif/prod do 30 i=2,n-1 f(i)=rat*(y(i-1)-2.0d0*y(i)+y(i+1))- & (y(i+1)-y(i-1))/(amult) 30 continue f(n)=dif*(y(n-1)-2.d0*y(n)+yN)/(dx*dx)-(yN-y(n-1))/(2.d0*dx) return end C SUBROUTINE JVANDER(N,X,Y,DFY,LDFY,EPS,IPAR) C --- JACOBIAN OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) COMMON /const/pi,dx,dif,ntop tn=x prod=dx*dx amult=2.0d0*dx a1=1.0d0/prod a2=-2.0d0/prod do 141 j=2,n dfy(1,j)=a1 141 continue u=a2 do 142 j=1,n sumer=0.0D+0 do 143 k=1,ntop sumer=sumer+dexp(-pi*pi*tn*k*k)*dsin(pi*j*dx*k) 143 continue dfy(2,j)=u+2.0D+0*y(j)-0.5D+0*sumer 142 continue do 144 j=1,n-1 dfy(3,j)=a1 144 continue return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 C * * * * * * * * * * * * * * * * * * * * * * * * * compile equation compile //venus/user/hairer/hailib/time cfeh driver_radau5 ~book2/appendix/radau5 ~book2/appendix/decsol equation time IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=15,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) dimension yend(nd),error(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) common/const/pi EXTERNAL FHIRES,JHIRES,SOLOUT c ------ FILE DE DONNEES ---------- open(6,file='ring on radau') write(6,3050) 3050 format(1x,'results on Ring problem') pi=4.0D+0*datan(1.0D+0) C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=1 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP sumh=0.0D+0 summ=0.0D+0 C --- DIMENSION OF THE SYSTEM N=15 C --- COMPUTE THE JACOBIAN IJAC=1 C --- JACOBIAN IS A FULL MATRIX MLJAC=N C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 C --- INITIAL VALUES X=0.0D0 do 171 ijk=1,15 y(ijk)=0.0D+0 171 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL*1.D-4 atol=rtol ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=1.0d-3 c CALL DTIME(TARRAY) it1=mclock() DO 20 I=1,1 C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) CALL RADAU5(N,FHIRES,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JHIRES,IJAC,MLJAC,MUJAC, & FHIRES,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c DO 15 K=1,8 c 15 WRITE (8,*)Y(K) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO XEND=XEND+100.D0 c CALL DTIME(TARRAY) yend(1)=-0.170799033d-1 yend(2)=-0.6660979d-2 yend(3)=+0.27531919d0 yend(4)=-0.39115732d0 yend(5)=-0.38851731d0 yend(6)=+0.27795920d0 yend(7)=+0.1114600281d0 yend(8)=+0.2979129627d-6 yend(9)=-0.31427403d-7 yend(10)=+0.70165883d-3 yend(11)=+0.85207538d-3 yend(12)=-0.77741454d-3 yend(13)=-0.77631967d-3 yend(14)=+0.78439426d-4 yend(15)=+0.252322784d-4 c sum=0.0D+0 do 270 k=1,15 error(k)=dabs(yend(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 continue it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) WRITE(6,*)TARRAY(1),summ,sumh/11.0d+0 summ=summ*atol sumh=sumh*atol/(11.0D+0) write(6,*) summ,sumh WRITE(6,*)NFE,NJE,NSTEP,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- WRITE (6,*)(ISTAT(J),J=14,20) WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) RETURN END C SUBROUTINE FHIRES (N, tn, Y, DY, RPAR, IPAR) INTEGER N DOUBLE PRECISION tn, Y(*), DY(*) c c implicit double precision (a-h,o-z) double precision y0, yN double precision kcns,pi,mi,dif,dx,amult,prod,rat common /const/pi double precision + c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc, + gamma,delta,sumer, + uin1,uin2,ud1,ud2,ud3,ud4,qud1,qud2,qud3,qud4 parameter (c=1.6d-8,cs=1d-9,cp=1d-8,r=25d3,rp=50d0, + lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4, + rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2, + gamma=40.67286402d-9,delta=17.7493332d0) uin1=0.5d0*sin(2.0d3*pi*tn) uin2=2.0d0*sin(2.0d4*pi*tn) ud1=+y(3)-y(5)-y(7)-uin2 ud2=-y(4)+y(6)-y(7)-uin2 ud3=+y(4)+y(5)+y(7)+uin2 ud4=-y(3)-y(6)+y(7)+uin2 qud1=gamma*(exp(delta*ud1)-1.0d0) qud2=gamma*(exp(delta*ud2)-1.0d0) qud3=gamma*(exp(delta*ud3)-1.0d0) qud4=gamma*(exp(delta*ud4)-1.0d0) dy(1)=(y(8)-0.5d0*y(10)+0.5d0*y(11)+y(14)-y(1)/r)/c dy(2)=(y(9)-0.5d0*y(12)+0.5d0*y(13)+y(15)-y(2)/r)/c dy(3)=(y(10)-qud1+qud4)/cs dy(4)=(-y(11)+qud2-qud3)/cs dy(5)=(y(12)+qud1-qud3)/cs dy(6)=(-y(13)-qud2+qud4)/cs dy(7)=(-y(7)/rp+qud1+qud2-qud3-qud4)/cp dy(8)=-y(1)/lh dy(9)=-y(2)/lh dy(10)=(0.5d0*y(1)-y(3)-rg2*y(10))/ls2 dy(11)=(-0.5d0*y(1)+y(4)-rg3*y(11))/ls3 dy(12)=(0.5d0*y(2)-y(5)-rg2*y(12))/ls2 dy(13)=(-0.5d0*y(2)+y(6)-rg3*y(13))/ls3 dy(14)=(-y(1)+uin1-(ri+rg1)*y(14))/ls1 dy(15)=(-y(2)-(rc+rg1)*y(15))/ls1 return end SUBROUTINE JHIRES(N,tn,Y,pd,LDFY,RPAR,IPAR) C -------- JACOBIAN FOR HIRES PROBLEM -------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),pd(LDFY,N) C------ METTRE A ZERO ------- c c implicit double precision (a-h,o-z) double precision kcns,pi,mi,dif,dx,amult,prod,rat,a1,a2,a3,u common /const/pi double precision + c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc, + gamma,delta,sumer, + uin2,ud1,ud2,ud3,ud4,qpud1,qpud2,qpud3,qpud4 parameter (c=1.6d-8,cs=1d-9,cp=1d-8,r=25d3,rp=50d0, + lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4, + rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2, + gamma=40.67286402d-9,delta=17.7493332d0) do 177 i=1,n do 177 j=1,n pd(i,j)=0.0D+0 177 continue uin2=2.0d0*sin(2.0d4*pi*tn) ud1=+y(3)-y(5)-y(7)-uin2 ud2=-y(4)+y(6)-y(7)-uin2 ud3=+y(4)+y(5)+y(7)+uin2 ud4=-y(3)-y(6)+y(7)+uin2 qpud1=gamma*delta*exp(delta*ud1) qpud2=gamma*delta*exp(delta*ud2) qpud3=gamma*delta*exp(delta*ud3) qpud4=gamma*delta*exp(delta*ud4) pd(1,1)=-1.0d0/(c*r) pd(1,8)=-1.0d0/c pd(1,10)=-0.5d0/c pd(1,11)=-pd(1,10) pd(1,14)=pd(1,8) pd(2,2)=pd(1,1) pd(2,9)=pd(1,8) pd(2,12)=pd(1,10) pd(2,13)=pd(1,11) pd(2,15)=pd(1,14) pd(3,3)=(-qpud1-qpud4)/cs pd(3,5)=qpud1/cs pd(3,6)=-qpud4/cs pd(3,7)=(qpud1+qpud4)/cs pd(3,10)=1.0d0/cs pd(4,4)=(-qpud2-qpud3)/cs pd(4,5)=-qpud3/cs pd(4,6)=qpud2/cs pd(4,7)=(-qpud2-qpud3)/cs pd(4,11)=-1.0d0/cs pd(5,3)=qpud1/cs pd(5,4)=-qpud3/cs pd(5,5)=(-qpud1-qpud3)/cs pd(5,7)=(-qpud1-qpud3)/cs pd(5,12)=1.0d0/cs pd(6,3)=-qpud4/cs pd(6,4)=qpud2/cs pd(6,6)=(-qpud2-qpud4)/cs pd(6,7)=(qpud2+qpud4)/cs pd(6,13)=-1.0d0/cs pd(7,3)=(qpud1+qpud4)/cp pd(7,4)=(-qpud2-qpud3)/cp pd(7,5)=(-qpud1-qpud3)/cp pd(7,6)=(qpud2+qpud4)/cp pd(7,7)=(-qpud1-qpud2-qpud3-qpud4-1.0d0/rp)/cp pd(8,1)=-1.0d0/lh pd(9,2)=pd(8,1) pd(10,1)=0.5d0/ls2 pd(10,3)=-1.0d0/ls2 pd(10,10)=-rg2/ls2 pd(11,1)=-0.5d0/ls3 pd(11,4)=1.0d0/ls3 pd(11,11)=-rg3/ls3 pd(12,2)=pd(10,1) pd(12,5)=pd(10,3) pd(12,12)=pd(10,10) pd(13,2)=pd(11,1) pd(13,6)=pd(11,4) pd(13,13)=pd(11,11) pd(14,1)=-1.0d0/ls1 pd(14,14)=-(ri+rg1)/ls1 pd(15,2)=pd(14,1) pd(15,15)=-(rc+rg1)/ls1 return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 C * * * * * * * * * * * * * * * * * * * * * * * * * compile equation compile //venus/user/hairer/hailib/time cfeh driver_radau5 ~book2/appendix/radau5 ~book2/appendix/decsol equation time IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=15,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) dimension yend(nd),error(nd) C -------- END PARAMETER LIST -------- REAL*4 TARRAY(2) common/const/pi EXTERNAL FHIRES,JHIRES,SOLOUT c ------ FILE DE DONNEES ---------- open(6,file='radhardring.res') write(6,3050) 3050 format(1x,'results on Ring problem') pi=4.0D+0*datan(1.0D+0) C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=3 NTOLMX=10 NTOLDF=1 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP sumh=0.0D+0 summ=0.0D+0 C --- DIMENSION OF THE SYSTEM N=15 C --- COMPUTE THE JACOBIAN IJAC=1 C --- JACOBIAN IS A FULL MATRIX MLJAC=N C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 C --- INITIAL VALUES X=0.0D0 do 171 ijk=1,15 y(ijk)=0.0D+0 171 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL*1.D-4 atol=rtol ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION iwork(2)=10000000 XEND=1.0d-3 c CALL DTIME(TARRAY) it1=mclock() DO 20 I=1,1 C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) CALL RADAU5(N,FHIRES,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JHIRES,IJAC,MLJAC,MUJAC, & FHIRES,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c DO 15 K=1,8 c 15 WRITE (8,*)Y(K) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO XEND=XEND+100.D0 c CALL DTIME(TARRAY) c yend(1)=-0.233905735842070174E-01 yend(2)= -0.736748548598215210E-02 yend(3)=0.258295670891733276 yend(4)=-0.406446572164108511 yend(5)=-0.403945566551075719 yend(6)=0.260796676505321678 yend(7)=0.110676186126208484 yend(8)=0.293990434238965346E-06 yend(9)=-0.284002993248608379E-07 yend(10)=0.726719826722729682E-03 yend(11)=0.792948719688146885E-03 yend(12)= -0.725528349561929123E-03 yend(13)=-0.794140196849026266E-03 yend(14)=0.708849541688199994E-04 yend(15)=0.239005907525775679E-04 sum=0.0D+0 do 270 k=1,15 error(k)=dabs(yend(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) write(6,*) y(k) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 continue it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) WRITE(6,*)TARRAY(1),summ,sumh/11.0d+0 summ=summ*atol sumh=sumh*atol/(11.0D+0) write(6,*) summ,sumh WRITE(6,*)NFE,NJE,NSTEP,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- WRITE (6,*)(ISTAT(J),J=14,20) WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) RETURN END C SUBROUTINE FHIRES (N, tn, Y, DY, RPAR, IPAR) INTEGER N DOUBLE PRECISION tn, Y(*), DY(*) c c implicit double precision (a-h,o-z) double precision y0, yN double precision kcns,pi,mi,dif,dx,amult,prod,rat common /const/pi double precision + c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc, + gamma,delta,sumer, + uin1,uin2,ud1,ud2,ud3,ud4,qud1,qud2,qud3,qud4 parameter (c=1.6d-8,cs=2.0d-12,cp=1d-8,r=25d3,rp=50d0, + lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4, + rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2, + gamma=40.67286402d-9,delta=17.7493332d0) uin1=0.5d0*sin(2.0d3*pi*tn) uin2=2.0d0*sin(2.0d4*pi*tn) ud1=+y(3)-y(5)-y(7)-uin2 ud2=-y(4)+y(6)-y(7)-uin2 ud3=+y(4)+y(5)+y(7)+uin2 ud4=-y(3)-y(6)+y(7)+uin2 qud1=gamma*(exp(delta*ud1)-1.0d0) qud2=gamma*(exp(delta*ud2)-1.0d0) qud3=gamma*(exp(delta*ud3)-1.0d0) qud4=gamma*(exp(delta*ud4)-1.0d0) dy(1)=(y(8)-0.5d0*y(10)+0.5d0*y(11)+y(14)-y(1)/r)/c dy(2)=(y(9)-0.5d0*y(12)+0.5d0*y(13)+y(15)-y(2)/r)/c dy(3)=(y(10)-qud1+qud4)/cs dy(4)=(-y(11)+qud2-qud3)/cs dy(5)=(y(12)+qud1-qud3)/cs dy(6)=(-y(13)-qud2+qud4)/cs dy(7)=(-y(7)/rp+qud1+qud2-qud3-qud4)/cp dy(8)=-y(1)/lh dy(9)=-y(2)/lh dy(10)=(0.5d0*y(1)-y(3)-rg2*y(10))/ls2 dy(11)=(-0.5d0*y(1)+y(4)-rg3*y(11))/ls3 dy(12)=(0.5d0*y(2)-y(5)-rg2*y(12))/ls2 dy(13)=(-0.5d0*y(2)+y(6)-rg3*y(13))/ls3 dy(14)=(-y(1)+uin1-(ri+rg1)*y(14))/ls1 dy(15)=(-y(2)-(rc+rg1)*y(15))/ls1 return end SUBROUTINE JHIRES(N,tn,Y,pd,LDFY,RPAR,IPAR) C -------- JACOBIAN FOR HIRES PROBLEM -------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),pd(LDFY,N) C------ METTRE A ZERO ------- c c implicit double precision (a-h,o-z) double precision kcns,pi,mi,dif,dx,amult,prod,rat,a1,a2,a3,u common /const/pi double precision + c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc, + gamma,delta,sumer, + uin2,ud1,ud2,ud3,ud4,qpud1,qpud2,qpud3,qpud4 parameter (c=1.6d-8,cs=2.0d-12,cp=1d-8,r=25d3,rp=50d0, + lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4, + rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2, + gamma=40.67286402d-9,delta=17.7493332d0) do 177 i=1,n do 177 j=1,n pd(i,j)=0.0D+0 177 continue uin2=2.0d0*sin(2.0d4*pi*tn) ud1=+y(3)-y(5)-y(7)-uin2 ud2=-y(4)+y(6)-y(7)-uin2 ud3=+y(4)+y(5)+y(7)+uin2 ud4=-y(3)-y(6)+y(7)+uin2 qpud1=gamma*delta*exp(delta*ud1) qpud2=gamma*delta*exp(delta*ud2) qpud3=gamma*delta*exp(delta*ud3) qpud4=gamma*delta*exp(delta*ud4) pd(1,1)=-1.0d0/(c*r) pd(1,8)=-1.0d0/c pd(1,10)=-0.5d0/c pd(1,11)=-pd(1,10) pd(1,14)=pd(1,8) pd(2,2)=pd(1,1) pd(2,9)=pd(1,8) pd(2,12)=pd(1,10) pd(2,13)=pd(1,11) pd(2,15)=pd(1,14) pd(3,3)=(-qpud1-qpud4)/cs pd(3,5)=qpud1/cs pd(3,6)=-qpud4/cs pd(3,7)=(qpud1+qpud4)/cs pd(3,10)=1.0d0/cs pd(4,4)=(-qpud2-qpud3)/cs pd(4,5)=-qpud3/cs pd(4,6)=qpud2/cs pd(4,7)=(-qpud2-qpud3)/cs pd(4,11)=-1.0d0/cs pd(5,3)=qpud1/cs pd(5,4)=-qpud3/cs pd(5,5)=(-qpud1-qpud3)/cs pd(5,7)=(-qpud1-qpud3)/cs pd(5,12)=1.0d0/cs pd(6,3)=-qpud4/cs pd(6,4)=qpud2/cs pd(6,6)=(-qpud2-qpud4)/cs pd(6,7)=(qpud2+qpud4)/cs pd(6,13)=-1.0d0/cs pd(7,3)=(qpud1+qpud4)/cp pd(7,4)=(-qpud2-qpud3)/cp pd(7,5)=(-qpud1-qpud3)/cp pd(7,6)=(qpud2+qpud4)/cp pd(7,7)=(-qpud1-qpud2-qpud3-qpud4-1.0d0/rp)/cp pd(8,1)=-1.0d0/lh pd(9,2)=pd(8,1) pd(10,1)=0.5d0/ls2 pd(10,3)=-1.0d0/ls2 pd(10,10)=-rg2/ls2 pd(11,1)=-0.5d0/ls3 pd(11,4)=1.0d0/ls3 pd(11,11)=-rg3/ls3 pd(12,2)=pd(10,1) pd(12,5)=pd(10,3) pd(12,12)=pd(10,10) pd(13,2)=pd(11,1) pd(13,6)=pd(11,4) pd(13,13)=pd(11,11) pd(14,1)=-1.0d0/ls1 pd(14,14)=-(ri+rg1)/ls1 pd(15,2)=pd(14,1) pd(15,15)=-(rc+rg1)/ls1 return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT VDPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=1999,LIWORK=3*ND+20) parameter(ijac=0,mljac=1,mujac=1,imas=0) parameter(ljac=mljac+mujac+1,le=2*mljac+mujac+1) parameter(lwork=nd*(ljac+3*le+12)+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) + ,true(nd),error(nd) common/add/sumh,rtol,atol C -------- END PARAMETER COMMON /const/pi,dx,ntop REAL*4 TARRAY(2) EXTERNAL FVANDER,JVANDER,SOLOUT open(6,file='radauhard.nobigres20') c ------ FILE DE DONNEES ---------- write(6,3500) 3500 format(1x,'radau5 on easy parabolic') RPAR=1.D-6 C --- LOOP FOR DIFFERENT TOLERANCES ntop=20 n=nd write(6,*) ntop,nd pi=4.0*DATAN(1.0D+0) c write(6,*) ntop c c dx=1.0d0/(n+1.0d0) NTOLMN=2 NTOLMX=10 NTOLDF=2 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM C --- COMPUTE THE JACOBIAN ANALYTICALLY C --- JACOBIAN IS A FULL MATRIX C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 c iout=1 C --- INITIAL VALUES X=0.0D0 do 201 j=1,n sumer=0.0D+0 do 21 jj=1,ntop sumer=sumer+dsin(pi*jj*j*dx) 21 continue y(j)=sumer 201 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=2.0D0 it1=mclock() C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) 220 continue CALL RADAU5(N,FVANDER,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JVANDER,IJAC,MLJAC,MUJAC, & FVANDER,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c WRITE (8,*) Y(1) c WRITE (8,*) Y(2) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO sum=0.0D+0 c tn=xend do 31 j=1,n true(j)=0.0D+0 do 33 i=1,ntop true(j)=true(j)+dexp(-tn*pi*pi*i*i)*dsin(dx*pi*i*j) 33 continue 31 continue c c do 270 k=1,1999 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue summ=dsqrt(sum/dble(nd)) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh us=log10(summ) ut=log10(sumh) write(6,*) us,ut WRITE (6,*)(ISTAT(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT common/add/sumh,rtol,atol dimension error(1999),true(1999) COMMON /const/pi,dx,ntop c WRITE (6,99) X,Y(1),Y(2),NR-1 tn=x do 31 j=1,n sumer=0.0D+0 do 33 i=1,ntop sumer=sumer+dexp(-tn*pi*pi*i*i)*dsin(dx*pi*i*j) 33 continue true(j)=sumer 31 continue 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) sum=0.0D+0 do 270 k=1,1999 error(k)=true(k)-y(k) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/1999.0D+0) if(sum.gt.sumh) sumh=sum c write(6,*) sum,sumh return END SUBROUTINE FVANDER(N,X,Y,F,EPS,IPAR) C --- RIGHT-HAND SIDE OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),F(N) COMMON /const/pi,dx,ntop tn=x amult=2.0d0*dx prod=dx*dx c SUMER=0.0D+0 do 103 i=1,ntop sumer=sumer+dexp(-pi*pi*tn*i*i)*dsin(pi*dx*i) 103 continue F(1)=(y(2)-2.0d0*y(1))/(dx*dx)+y(1)**2-0.5d+0*y(1)*sumer- + 0.5d+0*sumer**2 do 101 i=2,n-1 sumer=0.0d+0 do 104 ii=1,ntop sumer=sumer+dexp(-pi*pi*tn*ii*ii)*dsin(pi*i*ii*dx) 104 continue f(i)=(y(i-1)-2.0d0*y(i)+y(i+1))/prod + +y(I)**2-0.5D+0*y(i)*sumer-0.5D+0*sumer**2 101 continue sumer=0.0D+0 do 105 II=1,ntop sumer=sumer+dexp(-pi*pi*tn*ii*ii)*dsin(pi*n*dx*ii) 105 continue F(n)=(y(n-1)-2.0d0*y(n))/(dx*dx) + +y(n)**2-0.5D+0*y(n)*sumer-0.5D+0*sumer**2 return end C SUBROUTINE JVANDER(N,X,Y,DFY,LDFY,EPS,IPAR) C --- JACOBIAN OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) COMMON /const/pi,dx,ntop tn=x prod=dx*dx amult=2.0d0*dx a1=1.0d0/prod a2=-2.0d0/prod do 141 j=2,n dfy(1,j)=a1 141 continue u=a2 do 142 j=1,n sumer=0.0D+0 do 143 k=1,ntop sumer=sumer+dexp(-pi*pi*tn*k*k)*dsin(pi*j*dx*k) 143 continue dfy(2,j)=u+2.0D+0*y(j)-0.5D+0*sumer 142 continue do 144 j=1,n-1 dfy(3,j)=a1 144 continue return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT VDPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=1999,LIWORK=3*ND+20) parameter(ijac=0,mljac=1,mujac=1,imas=0) parameter(ljac=mljac+mujac+1,le=2*mljac+mujac+1) parameter(lwork=nd*(ljac+3*le+12)+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) + ,true(nd),error(nd) common/add/sumh,rtol,atol C -------- END PARAMETER COMMON /const/pi,dx,dif,ntop REAL*4 TARRAY(2) EXTERNAL FVANDER,JVANDER,SOLOUT open(6,file='radauflood.easynores') c ------ FILE DE DONNEES ---------- write(6,3500) 3500 format(1x,'radau5 on flood problem') RPAR=1.D-6 C --- LOOP FOR DIFFERENT TOLERANCES ntop=20 dif=2.0D-5 c dif=4.0D+0/500000.0D+0 n=nd write(6,*) dif,nd pi=4.0*DATAN(1.0D+0) dx=4.0d0/(n+1.0d0) NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM C --- COMPUTE THE JACOBIAN ANALYTICALLY C --- JACOBIAN IS A FULL MATRIX C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 c iout=1 C --- INITIAL VALUES X=0.0D0 do 201 j=1,n y(j)=dexp(-(-2.0D+0+j*dx)**2) 201 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=4.0D0 it1=mclock() C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) 220 continue CALL RADAU5(N,FVANDER,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JVANDER,IJAC,MLJAC,MUJAC, & FVANDER,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c WRITE (8,*) Y(1) c WRITE (8,*) Y(2) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO sum=0.0D+0 c tn=xend do 31 j=1,n true(j)= + (1.0d0/dsqrt(1.0d0+4.0d0*dif*tn))* & dexp(-((-2.0d0+j*dx-tn)**2)/(1.0d0+4.0d0*dif*tn)) 31 continue c c do 270 k=1,1999 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue summ=dsqrt(sum/dble(nd)) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh us=log10(summ) ut=log10(sumh) write(6,*) us,ut WRITE (6,*)(ISTAT(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT common/add/sumh,rtol,atol dimension error(1999),true(1999) COMMON /const/pi,dx,dif,ntop c WRITE (6,99) X,Y(1),Y(2),NR-1 tn=x do 31 j=1,n true(j)=(1.0d0/dsqrt(1.0d0+4.0d0*dif*tn))* & dexp(-((-2.0d0+j*dx-tn)**2)/(1.0d0+4.0d0*dif*tn)) 31 continue 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) sum=0.0D+0 do 270 k=1,1999 error(k)=true(k)-y(k) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/1999.0D+0) if(sum.gt.sumh) sumh=sum c write(6,*) sum,sumh return END SUBROUTINE FVANDER(N,X,Y,F,EPS,IPAR) C --- RIGHT-HAND SIDE OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),F(N) COMMON /const/pi,dx,dif,ntop tn=x amult=2.0d0*dx prod=dx*dx c y0=(1.0d0/dsqrt(1.d0+4.0d0*dif*tn))*dexp(-((-2.0d0-tn)**2)/ & (1.0d0+4.0d0*dif*tn)) yN=(1.0d0/dsqrt(1.d0+4.d0*dif*tn))*dexp(-((2.d0-tn)**2)/ & (1.0d0+4.0d0*dif*tn)) f(1)=dif*(y0-2.d0*y(1)+y(2))/(dx*dx)-(y(2)-y0)/(2.d0*dx) rat=dif/prod do 30 i=2,n-1 f(i)=rat*(y(i-1)-2.0d0*y(i)+y(i+1))- & (y(i+1)-y(i-1))/(amult) 30 continue f(n)=dif*(y(n-1)-2.d0*y(n)+yN)/(dx*dx)-(yN-y(n-1))/(2.d0*dx) return end C SUBROUTINE JVANDER(N,X,Y,DFY,LDFY,EPS,IPAR) C --- JACOBIAN OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) COMMON /const/pi,dx,dif,ntop tn=x prod=dx*dx amult=2.0d0*dx a1=1.0d0/prod a2=-2.0d0/prod do 141 j=2,n dfy(1,j)=a1 141 continue u=a2 do 142 j=1,n sumer=0.0D+0 do 143 k=1,ntop sumer=sumer+dexp(-pi*pi*tn*k*k)*dsin(pi*j*dx*k) 143 continue dfy(2,j)=u+2.0D+0*y(j)-0.5D+0*sumer 142 continue do 144 j=1,n-1 dfy(3,j)=a1 144 continue return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT VDPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=49999,LIWORK=3*ND+20) parameter(ijac=0,mljac=1,mujac=1,imas=0) parameter(ljac=mljac+mujac+1,le=2*mljac+mujac+1) parameter(lwork=nd*(ljac+3*le+12)+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),ISTAT(20) + ,true(nd),error(nd) common/add/sumh,rtol,atol C -------- END PARAMETER COMMON /const/pi,dx,ntop REAL*4 TARRAY(2) EXTERNAL FVANDER,JVANDER,SOLOUT open(6,file='radauhard.bigres5') c ------ FILE DE DONNEES ---------- write(6,3500) 3500 format(1x,'radau5 on hard parabolic') RPAR=1.D-6 C --- LOOP FOR DIFFERENT TOLERANCES ntop=5 n=nd write(6,*) ntop,nd pi=4.0*DATAN(1.0D+0) c write(6,*) ntop c c dx=1.0d0/(n+1.0d0) NTOLMN=2 NTOLMX=10 NTOLDF=1 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM C --- COMPUTE THE JACOBIAN ANALYTICALLY C --- JACOBIAN IS A FULL MATRIX C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 iout=1 C --- INITIAL VALUES X=0.0D0 do 201 j=1,n sumer=0.0D+0 do 21 jj=1,ntop sumer=sumer+dsin(pi*jj*j*dx) 21 continue y(j)=sumer 201 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-6 C --- SET DEFAULT VALUES DO I=1,20 WORK(I)=0.D0 IWORK(I)=0 ISTAT(I)=0 END DO C --- ENDPOINT OF INTEGRATION XEND=2.0D0 it1=mclock() C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) 220 continue CALL RADAU5(N,FVANDER,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JVANDER,IJAC,MLJAC,MUJAC, & FVANDER,IMAS,MLMAS,MUMAS, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- PRINT SOLUTION c WRITE (8,*) Y(1) c WRITE (8,*) Y(2) C --- PRINT STATISTICS DO J=14,20 ISTAT(J)=ISTAT(J)+IWORK(J) END DO sum=0.0D+0 c tn=xend do 31 j=1,n true(j)=0.0D+0 do 33 i=1,ntop true(j)=true(j)+dexp(-tn*pi*pi*i*i)*dsin(dx*pi*i*j) 33 continue 31 continue c c do 270 k=1,49999 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue summ=dsqrt(sum/dble(nd)) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh us=log10(summ) ut=log10(sumh) write(6,*) us,ut WRITE (6,*)(ISTAT(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (ISTAT(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),CONT(LRC) COMMON /INTERN/XOUT common/add/sumh,rtol,atol dimension error(49999),true(49999) COMMON /const/pi,dx,ntop c WRITE (6,99) X,Y(1),Y(2),NR-1 tn=x do 31 j=1,n sumer=0.0D+0 do 33 i=1,ntop sumer=sumer+dexp(-tn*pi*pi*i*i)*dsin(dx*pi*i*j) 33 continue true(j)=sumer 31 continue 99 FORMAT(1X,'X =',F6.2,' Y =',2F18.7,' NSTEP =',I4) sum=0.0D+0 do 270 k=1,49999 error(k)=true(k)-y(k) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/49999.0D+0) if(sum.gt.sumh) sumh=sum c write(6,*) sum,sumh return END SUBROUTINE FVANDER(N,X,Y,F,EPS,IPAR) C --- RIGHT-HAND SIDE OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),F(N) COMMON /const/pi,dx,ntop tn=x amult=2.0d0*dx prod=dx*dx c SUMER=0.0D+0 do 103 i=1,ntop sumer=sumer+dexp(-pi*pi*tn*i*i)*dsin(pi*dx*i) 103 continue F(1)=(y(2)-2.0d0*y(1))/(dx*dx)+y(1)**2-0.5d+0*y(1)*sumer- + 0.5d+0*sumer**2 do 101 i=2,n-1 sumer=0.0d+0 do 104 ii=1,ntop sumer=sumer+dexp(-pi*pi*tn*ii*ii)*dsin(pi*i*ii*dx) 104 continue f(i)=(y(i-1)-2.0d0*y(i)+y(i+1))/prod + +y(I)**2-0.5D+0*y(i)*sumer-0.5D+0*sumer**2 101 continue sumer=0.0D+0 do 105 II=1,ntop sumer=sumer+dexp(-pi*pi*tn*ii*ii)*dsin(pi*n*dx*ii) 105 continue F(n)=(y(n-1)-2.0d0*y(n))/(dx*dx) + +y(n)**2-0.5D+0*y(n)*sumer-0.5D+0*sumer**2 return end C SUBROUTINE JVANDER(N,X,Y,DFY,LDFY,EPS,IPAR) C --- JACOBIAN OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(LDFY,N) COMMON /const/pi,dx,ntop tn=x prod=dx*dx amult=2.0d0*dx a1=1.0d0/prod a2=-2.0d0/prod do 141 j=2,n dfy(1,j)=a1 141 continue u=a2 do 142 j=1,n sumer=0.0D+0 do 143 k=1,ntop sumer=sumer+dexp(-pi*pi*tn*k*k)*dsin(pi*j*dx*k) 143 continue dfy(2,j)=u+2.0D+0*y(j)-0.5D+0*sumer 142 continue do 144 j=1,n-1 dfy(3,j)=a1 144 continue return end