C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=2,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (NWKJAC=400,NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJF,SOLSF,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FVANDER,JVANDER,SPRMOM c ------ FILE DE DONNEES ---------- C --- LOOP FOR DIFFERENT TOLERANCES write(6,3500) 3500 format(1x,'vdp by sprint') 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 --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 C --- INITIAL VALUES X=0.0D0 Y(1)=2.0D0 Y(2)=0.0D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION XEND=1.0D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 2, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF DO 20 I=1,11 C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FVANDER,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJF, WKJAC, NWKJAC, JACPVT, JVANDER, SOLSF,SPRMOM,WKMON, 3 NWKMON) C --- PRINT SOLUTION c WRITE (6,*) Y(1) c WRITE (6,*) Y(2) ID=0 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(ndim)) 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,*)ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV C WRITE(6,99)N,T,HLAST,H,Y(1),Y(2),Y(3) C 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',3F15.8) RETURN END c C SUBROUTINE FVANDER( N, T, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE EPS=1.0D-6 R(1)=Y(2) - YDOT(1) PROD=1.D0-Y(1)*Y(1) R(2)=(PROD*Y(2)-Y(1))/EPS - YDOT(2) END IF RETURN END C SUBROUTINE JVANDER(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) FAC=H*D EPS=1.0D-6 DFY(1,1)= -1.D0 DFY(1,2)=FAC DFY(2,1)=FAC*(-2.0D0*Y(1)*Y(2)-1.0D0)/EPS DFY(2,2)=FAC*(1.0D0-Y(1)**2)/EPS -1.D0 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=3,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (NWKJAC=400,NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJF,SOLSF,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FROBER,JROBER,SPRMOM c ------ FILE DE DONNEES ---------- write(6,3050) 3050 format(1x,'results on rob by sprint') 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 --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 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.d-6*RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION XEND=1.0D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 2, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF DO 20 I=1,12 C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FROBER,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJF, WKJAC, NWKJAC, JACPVT, JROBER, SOLSF,SPRMOM,WKMON, 3 NWKMON) C --- PRINT SOLUTION c WRITE (6,*) Y(1) c WRITE (6,*) Y(2) c WRITE (6,*) Y(3) ID=0 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(ndim)) 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,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE(6,*)TARRAY(1),summ,sumh/12.0D+0 summ=summ*atol sumh=sumh*atol/12.0D+0 write(6,*) SUMM,SUMH WRITE(6,*)ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV C WRITE(6,99)N,T,HLAST,H,Y(1),Y(2),Y(3) C 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',3F15.8) RETURN END c C SUBROUTINE FROBER( N, T, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE R(1)=-0.04D0*Y(1)+1.D4*Y(2)*Y(3) - YDOT(1) R(3)=3.D7*Y(2)*Y(2) - YDOT(3) R(2)=0.04D0*Y(1)-1.D4*Y(2)*Y(3)-3.D7*Y(2)*Y(2) - YDOT(2) END IF RETURN END C SUBROUTINE JROBER(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) FAC=H*D PROD1=1.0D4*Y(2) PROD2=1.0D4*Y(3) PROD3=6.0D7*Y(2) DFY(1,1)=FAC*(-0.04D0 )-1.D0 DFY(1,2)=FAC*(PROD2 ) DFY(1,3)=FAC*(PROD1 ) DFY(2,1)=FAC*(0.04D0 ) DFY(2,2)=FAC*(-PROD2-PROD3 )-1.D0 DFY(2,3)=FAC*(-PROD1 ) DFY(3,1)=0.D0 DFY(3,2)=FAC*(PROD3 ) DFY(3,3)=0.D0 -1.D0 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=8,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (NWKJAC=NDIM*NDIM+2,NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJF,SOLSF,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FHIRES,JHIRES c ------ FILE DE DONNEES ---------- C --- LOOP FOR DIFFERENT TOLERANCES write(6,3050) 3050 format(1x,'results on hires by sprint') 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 -------- N=NDIM C --- SET DEFAULT VALUES DO 12 I = 1,21 12 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 14 I = 5,14+N INFORM(I) = 0 14 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ISET=1 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=1.0D-4*RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION XEND=321.8122D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 2, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF DO 20 I=1,2 C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FHIRES,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJF, WKJAC, NWKJAC, JACPVT, JHIRES, SOLSF,SPRMON,WKMON, 3 NWKMON) C --- PRINT SOLUTION c DO 15 K=1,8 c 15 WRITE (6,*)Y(K) XEND=XEND+100.D0 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(ndim)) 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,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE(6,*)TARRAY(1),summ,sumh/2.0D+0 summ=summ*atol sumh=sumh*atol/2.0D+0 write(6,*) summ,sumh WRITE(6,*)NFE,NJE,NSTEP,ID,ID,NDEC,NBSOL C -------- NEW TOLERANCE --- ID=0 WRITE(6,*)ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c C SUBROUTINE FHIRES( N, T, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) COMMON/NERVES/NNERV COMMON/DIFFCOEF/DIFFUS IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE c----------- R(1) = -1.71D0*Y(1) + 0.43D0*Y(2) + 8.32D0*Y(3)+0.0007D0-YDOT(1) R(2) = 1.71D0*Y(1) - 8.75D0*Y(2) -YDOT(2) R(3) = -10.03D0*Y(3) + 0.43D0*Y(4) + 0.035D0*Y(5) -YDOT(3) R(4) = 8.32D0*Y(2) + 1.71D0*Y(3) - 1.12D0*Y(4) -YDOT(4) R(5) = -1.745D0*Y(5) + 0.43D0*Y(6) + 0.43D0*Y(7) -YDOT(5) R(6) = -280.0D0*Y(6)*Y(8) + 0.69D0*Y(4) + 1.71D0*Y(5) - & 0.43D0*Y(6) + 0.69D0*Y(7) -YDOT(6) R(7) = 280.0D0*Y(6)*Y(8) - 1.81D0*Y(7) -YDOT(7) R(8) = - 280.0D0*Y(6)*Y(8) + 1.81D0*Y(7) -YDOT(8) END IF RETURN END SUBROUTINE JHIRES(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) FAC=H*D DO 1 I=1,N DO 1 J=1,N 1 DFY(I,J)=0.D0 C DFY(1,1)=FAC*( -1.71D0 ) -1.D0 DFY(1,2)=FAC*( 0.43D0 ) DFY(1,3)=FAC*( + 8.32D0 ) C DFY(2,1)=FAC*( 1.71D0 ) DFY(2,2)=FAC*( - 8.75D0 ) -1.D0 C DFY(3,3)=FAC*( -10.03D0 ) -1.D0 DFY(3,4)=FAC*( 0.43D0 ) DFY(3,5)=FAC*( + 0.035D0 ) C DFY(4,2)=FAC*( 8.32D0 ) DFY(4,3)=FAC*( + 1.71D0 ) DFY(4,4)=FAC*( - 1.12D0 ) -1.D0 C DFY(5,5)=FAC*( -1.745D0 ) -1.D0 DFY(5,6)=FAC*( + 0.43D0 ) DFY(5,7)=FAC*( + 0.43D0 ) C DFY(6,4)=FAC*( + 0.69D0 ) DFY(6,5)=FAC*( + 1.71D0 ) DFY(6,6)=FAC*( - .43D0 -280.0D0*Y(8)) -1.D0 DFY(6,7)=FAC*( + 0.69D0 ) DFY(6,8)=FAC*( -280.0D0*Y(6) ) C DFY(7,6)=FAC*( 280.0D0*Y(8) ) DFY(7,7)=FAC*( -1.81D0 ) -1.D0 DFY(7,8)=FAC*( 280.0D0*Y(6) ) C DFY(8,6)=FAC*( -280.0D0*Y(8) ) DFY(8,7)=FAC*( 1.81D0 ) DFY(8,8)=FAC*( -280.0D0*Y(6) ) -1.D0 C RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=3,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (NWKJAC=400,NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJF,SOLSF,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FOREGON,JOREGON,SPRMOM c ------ FILE DE DONNEES ---------- write(6,3050) 3050 format(1x,'sprint results on oreg') 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 --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 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.d-6*RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION XEND=30.0D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 2, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF DO 20 I=1,12 C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FOREGON,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJF, WKJAC, NWKJAC, JACPVT, JOREGON, SOLSF,SPRMOM,WKMON, 3 NWKMON) C --- PRINT SOLUTION c WRITE (6,*) Y(1) c WRITE (6,*) Y(2) c WRITE (6,*) Y(3) XEND=XEND+30.D0 c :r bit.f 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(ndim)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE 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 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' summ=summ*atol sumh=sumh*atol/12.0D+0 write(6,*) summ,sumh ID=0 WRITE(6,*)ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV C WRITE(6,99)N,T,HLAST,H,Y(1),Y(2) C 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',2F16.5) RETURN END c C SUBROUTINE FOREGON( N, T, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE R(1)=77.27D0*(Y(2)+Y(1)*(1.D0-8.375D-6*Y(1)-Y(2)))- YDOT(1) R(2)=(Y(3)-(1.D0+Y(1))*Y(2))/77.27D0 - YDOT(2) R(3)=0.161D0*(Y(1)-Y(3)) - YDOT(3) END IF RETURN END C SUBROUTINE JOREGON(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) EPS=1.0D-6 FAC=H*D DFY(1,1)=FAC*(77.27D0*(1.D0-2.D0*8.375D-6*Y(1)-Y(2)))-1.D0 DFY(1,2)=FAC*(77.27D0*(1.D0-Y(1)) ) DFY(1,3)=0.D0 DFY(2,1)=FAC*(-Y(2)/77.27D0 ) DFY(2,2)=FAC*(-(1.D0+Y(1))/77.27D0 )-1.D0 DFY(2,3)=FAC*(1.D0/77.27D0 ) DFY(3,1)=FAC*(.161D0 ) DFY(3,2)=.0D0 DFY(3,3)=FAC*(-.161D0 )-1.D0 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=4,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (NWKJAC=400,NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJF,SOLSF,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FE5,JE5,SPRMOM c ------ FILE DE DONNEES ---------- write(6,3050) 3050 format(1x,'e5 by sprint') 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 C --- DIMENSION OF THE SYSTEM N=4 C --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 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=1 C --- ENDPOINT OF INTEGRATION XEND=10.0D0 c CALL DTIME(TARRAY) it1= mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 2, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF DO 20 I=1,7 C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FE5,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJF, WKJAC, NWKJAC, JACPVT, JE5, SOLSF,SPRMOM,WKMON, 3 NWKMON) C --- PRINT SOLUTION c WRITE (6,*) Y(1) c WRITE (6,*) Y(2) c WRITE (6,*) Y(3) c WRITE (6,*) Y(4) ID=0 XEND=XEND*100.D0 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(ndim)) 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 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' ID=0 c CALL DTIME(TARRAY) WRITE(6,*)ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV C WRITE(6,99)N,T,HLAST,H,Y(1),Y(2),Y(3) C 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',3F15.8) RETURN END c C SUBROUTINE FE5( N, T, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE 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) F2=PROD1-PROD3 F4=PROD2-PROD4 R(1)=-PROD1-PROD2 - YDOT(1) R(2)=F2 - YDOT(2) R(4)=F4 - YDOT(4) R(3)=F2-F4 - YDOT(3) END IF RETURN END C SUBROUTINE JE5(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) FAC=H*D A=7.89D-10 B=1.1D7 CM=1.13D9 C=1.13D3 DFY(1,1)=FAC*(-A-B*Y(3) )-1.D0 DFY(1,2)=0.0D0 DFY(1,3)=FAC*(-B*Y(1) ) DFY(1,4)=0.0D0 DFY(2,1)=FAC*(A ) DFY(2,2)=FAC*(-CM*Y(3) )-1.D0 DFY(2,3)=FAC*(-CM*Y(2) ) DFY(2,4)=0.0D0 DFY(3,1)=FAC*(A-B*Y(3)) DFY(3,2)=FAC*(-CM*Y(3) ) DFY(3,3)=FAC*(-B*Y(1)-CM*Y(2)) -1.D0 DFY(3,4)=FAC*C DFY(4,1)=FAC*(B*Y(3)) DFY(4,2)=0.0D0 DFY(4,3)=FAC*(B*Y(1)) DFY(4,4)=FAC*(-C) -1.D0 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=80,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (NWKJAC=NDIM*NDIM+2,NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJF,SOLSF,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- COMMON/NNNN/NCOM,NNCOM,NSQ,NQUATR,DELTAS EXTERNAL FTIGE,JTIGE c ------ FILE DE DONNEES ---------- write(6,3050) 3050 format(1x,' sprint on beam') 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 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,21 12 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 14 I = 5,14+N INFORM(I) = 0 14 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ISET=0 C --------- INITIAL VALUES ------------- T=0.D0 DO 1 I=1,NN 1 Y(I)=0.D0 TEND=5.D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 2, NN, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF C --- CALL OF SPRINT CALL SPRINT( NN, T, TEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FTIGE,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJF, WKJAC, NWKJAC, JACPVT, JTIGE, SOLSF,SPRMON,WKMON, 3 NWKMON) 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 --- PRINT SOLUTION WRITE(6,9921)(Y(I),I=1,NN) 9921 FORMAT(1X,F22.16) ID=0 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE (6,*) ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' IF(TARRAY(1).GT.300.)GOTO 66 C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE 66 CONTINUE STOP END C SUBROUTINE FTIGE( NN, T, TH, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION YDOT(NN), R(NN), RWORK(NR) DIMENSION TH(150),U(150),V(150),W(150) DIMENSION ALPHA(150),BETA(150),STH(150),CTH(150) COMMON/NNNN/N,NNCOM,NSQ,NQUATR,DELTAS IF(IRES .EQ. -1)THEN DO 10 I = 1,NN 10 R(I) = - YDOT(I) ELSE c----------- 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 R(I)=TH(N+I) -YDOT(I) 54 R(N+I)=U(I) -YDOT(N+I) END IF RETURN END SUBROUTINE JTIGE(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) WRITE(6,*)' OHO, FROM WHERE DO I COME ???' STOP END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * compile equation_sprint compile //venus/user/hairer/hailib/time cfeh driver_sprint ~fortran/sprint/sprint_collection equation_sprint time IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(MX=8,MY=5,MACHS1=2,MACHS2=4) PARAMETER (NDIM=2*MX*MY,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (NWKJAC=NDIM*NDIM+2,NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJF,SOLSF,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FPLATE,JPLATE 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,3050) 3050 format('plate by sprint') 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 -------- N=NDIM C --- SET DEFAULT VALUES DO 12 I = 1,21 12 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 14 I = 5,14+N INFORM(I) = 0 14 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ISET=1 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=1 C --- ENDPOINT OF INTEGRATION XEND=7.D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 2, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FPLATE,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJF, WKJAC, NWKJAC, JACPVT, JPLATE, SOLSF,SPRMON,WKMON, 3 NWKMON) 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 ID=0 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE(6,*)ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE 77 STOP END C SUBROUTINE FPLATE( N, X, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE c----------- C -------- LA BOUCLE ------- DO 1 I=1,NX DO 1 J=1,NY K=I+NX*(J-1) C -------- DERIVEE DEUXIEME ---- R(K)=Y(K+NDEMI) -YDOT(K) 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 R(K+NDEMI)=-OMEGA*Y(K+NDEMI)-FAC*UC+FORCE*WEIGHT-YDOT(K+NDEMI) 1 CONTINUE END IF RETURN END SUBROUTINE JPLATE(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT HFD=H*D 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)=HFD*(1.D0) C ------ POINT CENTRAL --- DFY(K+NDEMI,K)=HFD*(-FAC*16.D0) IF(I.GT.1)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-HFD*(FAC) DFY(K+NDEMI,K-1)=HFD*(FAC*8.D0) END IF IF(I.LT.NX)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-HFD*(FAC) DFY(K+NDEMI,K+1)=HFD*(FAC*8.D0) END IF IF(J.GT.1)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-HFD*(FAC) DFY(K+NDEMI,K-NX)=HFD*(FAC*8.D0) END IF IF(J.LT.NY)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-HFD*(FAC) DFY(K+NDEMI,K+NX)=HFD*(FAC*8.D0) END IF IF(I.GT.1 .AND.J.GT.1 )DFY(K+NDEMI,K-NX-1)=HFD*(-FAC*2.D0) IF(I.LT.NX.AND.J.GT.1 )DFY(K+NDEMI,K-NX+1)=HFD*(-FAC*2.D0) IF(I.GT.1 .AND.J.LT.NY)DFY(K+NDEMI,K+NX-1)=HFD*(-FAC*2.D0) IF(I.LT.NX.AND.J.LT.NY)DFY(K+NDEMI,K+NX+1)=HFD*(-FAC*2.D0) IF(I.GT.2)DFY(K+NDEMI,K-2)=HFD*(-FAC) IF(I.LT.NXM1)DFY(K+NDEMI,K+2)=HFD*(-FAC) IF(J.GT.2)DFY(K+NDEMI,K-2*NX)=HFD*(-FAC) IF(J.LT.NYM1)DFY(K+NDEMI,K+2*NX)=HFD*(-FAC) DFY(K+NDEMI,K+NDEMI)= HFD*(-OMEGA) 2 CONTINUE C ------- SOUSTRAIRE LA DIAGONALE -------- DO 3 I=1,N DFY(I,I)=DFY(I,I)-1.D0 3 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=96,NRW=21+3*NDIM,NWKRES=4,NY2DIM=20) PARAMETER (NWKJAC=10*NDIM+2,NWKMON=NDIM) PARAMETER (NIWJAC=NDIM,NYSAVE=NDIM) dimension error(ndim),true(ndim) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJB,SOLSB,MONITR,SBLEND,SPGEAR,STHETA,SPRMOM C -------- END PARAMETER LIST ---------- COMMON/NERVES/NNERV COMMON/DIFFCOEF/DIFFUS EXTERNAL FCUSP,JCUSP c ------ FILE DE DONNEES ---------- write(6,3029) 3029 format(1x,'sprint results for cusp 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 -------- INITIAL CONSTANTES -------- NNERV=32 DIFFUS=1.D0*NNERV*NNERV/144.D0 N=3*NNERV C --- SET DEFAULT VALUES DO 12 I = 1,21 12 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 14 I = 5,14+N INFORM(I) = 0 14 CONTINUE SNORM = 'AVERL2' ML=3 MU=3 C --- COMPUTE THE JACOBIAN ISET=0 C ---------- VAL INIT ------------ X=0.D0 XEND=1.1D0 ANERV=NNERV DEL=2.D0*3.14159265358979324D0/ANERV DO 18 INERV=1,NNERV Y(3*INERV-2)=0.D0 Y(3*INERV-1)=-2.D0*COS(INERV*DEL) 18 Y(3*INERV)=2.D0*SIN(INERV*DEL) C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION XEND=1.1D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 1, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FCUSP,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJB, WKJAC, NWKJAC, JACPVT, JCUSP, SOLSB,SPRMOM,WKMON, 3 NWKMON) 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) ID=0 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE (6,*) ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- IF(TARRAY(1).GT.100.)STOP TOLST=TOLST*TOLFC 30 CONTINUE STOP END C c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV c WRITE(6,99)N,T,HLAST,H,Y(1),Y(2) c 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',2F13.8) RETURN END C C SUBROUTINE FCUSP( N, T, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) COMMON/NERVES/NNERV COMMON/DIFFCOEF/DIFFUS IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE 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 R(3*INERV-2)=XDOT+DIFFUS*(XLEFT-2.D0*X+XRIGHT)- YDOT(3*INERV-2) R(3*INERV-1)=ADOT+DIFFUS*(ALEFT-2.D0*A+ARIGHT)- YDOT(3*INERV-1) R(3*INERV) =BDOT+DIFFUS*(BLEFT-2.D0*B+BRIGHT)- YDOT(3*INERV) 25 CONTINUE END IF RETURN END SUBROUTINE JCUSP(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) WRITE(6,*)' OHO, FROM WHERE DO I COME ???' STOP END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=1000,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (MU=2,ML=2,NWKJAC=(2*ML+MU+1)*NDIM+2) PARAMETER (NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJB,SOLSB,MONITR,SBLEND,SPGEAR,STHETA,SPRMOM C -------- END PARAMETER LIST ---------- EXTERNAL FBRUS,JBRUS COMMON/PARAM/N,N2,GAMMA,GAMMA2 c ------ FILE DE DONNEES ---------- write(6,3050) 3050 format(1x,'bruss results by sprint') 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 -------- PI=3.14159265358979324D0 N=500 C --- DIMENSION OF THE SYSTEM N2=2*N USDELQ=(DBLE(N+1))**2 GAMMA=0.02D0*USDELQ GAMMA2=2.D0*GAMMA C --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 ISET=0 C --- INITIAL VALUES X=0.D0 DO 1 I=1,N ANP1=N+1 XI=I/ANP1 Y(2*I)=3.D0 1 Y(2*I-1)=1.D0+0.5D0*DSIN(2.D0*PI*XI) C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION XEND=10.D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 1, N2, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF C --- CALL OF SPRINT CALL SPRINT( N2, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FBRUS,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJB, WKJAC, NWKJAC, JACPVT, JBRUS, SOLSB,SPRMOM,WKMON, 3 NWKMON) c CALL DTIME(TARRAY) true(1)=0.9949197002317599 true(8)=3.0213845767604077 true(15)=0.9594350193986054 true(22)=3.0585989778165419 true(29)=0.9243010095428502 true(36)=3.0952478919989637 true(43)=0.8897959106772672 true(50)=3.1310118289054731 true(57)=0.8561653620284367 true(64)=3.1656101198770159 true(71)=0.8236197147449046 true(78)=3.1988043370624344 true(85)=0.7923328094811884 true(92)=3.2303999530641514 true(99)=0.7624421042573115 true(106)=3.2602463873623941 true(113)=0.7340499750795348 true(120)=3.2882356529108807 true(127)=0.7072259700779899 true(134)=3.3142998590079271 true(141)=0.6820097782458483 true(148)=3.3384078449410937 true(155)=0.6584146743834650 true(162)=3.3605612157873943 true(169)=0.6364312187752559 true(176)=3.3807900316323134 true(183)=0.6160310186921587 true(190)=3.3991483695914764 true(197)=0.5971703941198909 true(204)=3.4157099395342736 true(211)=0.5797938277687891 true(218)=3.4305638938070224 true(225)=0.5638371159206763 true(232)=3.4438109320334580 true(239)=0.5492301695479158 true(246)=3.4555597666485198 true(253)=0.5358994429426996 true(260)=3.4659239846027008 true(267)=0.5237699892215797 true(274)=3.4750193162238476 true(281)=0.5127671585747183 true(288)=3.4829613034792271 true(295)=0.5028179665048467 true(302)=3.4898633463634923 true(309)=0.4938521662914935 true(316)=3.4958350971335204 true(323)=0.4858030633656755 true(330)=3.5009811668111510 true(337)=0.4786081100251151 TRUE(344)=3.5054001059792705 true(351)=0.4722093177200750 true(358)=3.5091836216744015 true(365)=0.4665535216425440 true(372)=3.5124159935026285 true(379)=0.4615925290790646 true(386)=3.5151736544621075 true(393)=0.4572831793403656 true(400)=3.5175249049438184 true(407)=0.4535873393501199 true(414)=3.5195297317024448 true(421)=0.4504718553589467 true(428)=3.5212397070273984 true(435)=0.4479084778719241 true(442)=3.5226979467564341 true(449)=0.4458737738041973 true(456)=3.5239391090719634 true(463)=0.4443490371324889 true(470)=3.5249894191569453 true(477)=0.4433202068820853 true(484)=3.5258667077466495 true(491)=0.4427777991494095 true(498)=3.5265804544017270 true(505)=0.4427168579654424 true(512)=3.5271318289682063 true(519)=0.4431369281018266 true(526)=3.5275137272135266 true(533)=0.4440420513508381 true(540)=3.5277107990730161 true(547)=0.4454407863109616 true(554)=3.5276994703501980 true(561)=0.4473462502188303 true(568)=3.5274479611304068 true(575)=0.4497761798232572 true(582)=3.5269163066324394 true(589)=0.4527530066369863 true(596)=3.5260563887768472 true(603)=0.4563039400688689 true(610)=3.5248119894251024 true(617)=0.4604610498812091 true(624)=3.5231188790654930 true(631)=0.4652613370907894 true(638)=3.5209049576992761 true(645)=0.4707467798082714 true(652)=3.5180904678044698 true(659)=0.4769643375804777 true(666)=3.5145883024867057 true(673)=0.4839658945842979 true(680)=3.5103044351908528 true(687)=0.4918081185812277 true(694)=3.5051385005173827 true(701)=0.5005522089940899 true(708)=3.4989845585737802 true(715)=0.5102635039989190 true(722)=3.4917320776245013 true(729)=0.5210109134090777 true(736)=3.4832671712209993 true(743)=0.5328661417420937 true(750)=3.4734741260299615 true(757)=0.5459026646938675 true(764)=3.4622372546582585 true(771)=0.5601944229089820 true(778)=3.4494431032230182 true(785)=0.5758142001453760 true(792)=3.4349830354873361 true(799)=0.5928316594749734 true(806)=3.4187562033108012 true(813)=0.6113110218368440 true(820)=3.4006728962523969 true(827)=0.6313083867734524 true(834)=3.3806582409098729 true(841)=0.6528687160104193 true(848)=3.3586561928427350 true(855)=0.6760225267555723 true(862)=3.3346337311179157 true(869)=0.7007823726569721 true(876)=3.3085851288057930 true(883)=0.7271392249346637 true(890)=3.2805361342349380 true(897)=0.7550589020044152 true(904)=3.2505478606008622 true(911)=0.7844787296769868 true(918)=3.2187201496972175 true(925)=0.8153046416214843 true(932)=3.1851941538893653 true(939)=0.8474089465959840 true(946)=3.1501538739882800 true(953)=0.8806289904192589 true(960)=3.1138264039027113 true(967)=0.9147669230929857 true(974)=3.0764806689389470 true(981)=0.9495907429372025 true(988)=3.0384245041548366 true(995)=0.9848367306701233 sum=0.0D+0 do 270 k=1,995,7 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/143.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 ID=0 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE (6,*) ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV c WRITE(6,99)N,T,HLAST,H,Y(1),Y(2) c 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',2F13.8) RETURN END C C SUBROUTINE FBRUS( N, T, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) COMMON/PARAM/NHALF,N2,GAMMA,GAMMA2 IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE I=1 IU=2*I-1 IV=2*I UI=Y(IU) VI=Y(IV) UIM=1.D0 VIM=3.D0 UIP=Y(IU+2) VIP=Y(IV+2) PROD=UI*UI*VI R(IU)=1.D0+PROD-4.D0*UI+GAMMA*(UIM-2.D0*UI+UIP) -YDOT(IU) R(IV)=3.D0*UI-PROD+GAMMA*(VIM-2.D0*VI+VIP) -YDOT(IV) DO 5 I=2,NHALF-1 IU=2*I-1 IV=2*I UI=Y(IU) VI=Y(IV) UIM=Y(IU-2) VIM=Y(IV-2) UIP=Y(IU+2) VIP=Y(IV+2) PROD=UI*UI*VI R(IU)=1.D0+PROD-4.D0*UI+GAMMA*(UIM-2.D0*UI+UIP) -YDOT(IU) R(IV)=3.D0*UI-PROD+GAMMA*(VIM-2.D0*VI+VIP) -YDOT(IV) 5 CONTINUE I=NHALF IU=2*I-1 IV=2*I UI=Y(IU) VI=Y(IV) UIM=Y(IU-2) VIM=Y(IV-2) UIP=1.D0 VIP=3.D0 PROD=UI*UI*VI R(IU)=1.D0+PROD-4.D0*UI+GAMMA*(UIM-2.D0*UI+UIP) -YDOT(IU) R(IV)=3.D0*UI-PROD+GAMMA*(VIM-2.D0*VI+VIP) -YDOT(IV) END IF RETURN END C SUBROUTINE JBRUS(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) COMMON/PARAM/NHALF,N2,GAMMA,GAMMA2 FAC=H*D DO 1 I=1,NHALF IU=2*I-1 IV=2*I UI=Y(IU) VI=Y(IV) UIVI=UI*VI UI2=UI*UI DFY(3,IU)=FAC*(2.D0*UIVI-4.D0-GAMMA2) -1.D0 DFY(2,IV)=FAC*UI2 DFY(4,IU)=FAC*(3.D0-2.D0*UIVI) DFY(3,IV)=FAC*(-UI2-GAMMA2) -1.D0 DFY(2,IU)=0.D0 DFY(4,IV)=0.D0 1 CONTINUE FAC1=FAC*GAMMA DO 2 I=1,N2-2 DFY(1,I+2)=FAC1 DFY(5,I) =FAC1 2 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=15,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (NWKJAC=400,NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension yend(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) common/const/pi EXTERNAL PREPJF,SOLSF,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FVANDER,JVANDER,SPRMOM c ------ FILE DE DONNEES ---------- C --- LOOP FOR DIFFERENT TOLERANCES write(6,3500) 3500 format(1x,'vdp by sprint') pi=4.0D+0*datan(1.0D+0) 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=15 C --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 C --- INITIAL VALUES X=0.0D0 do 151 ijk=1,15 y(ijk)=0.0D+0 151 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION XEND=1.0D-3 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 2, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF DO 20 I=1,1 C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FVANDER,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJF, WKJAC, NWKJAC, JACPVT, JVANDER, SOLSF,SPRMOM,WKMON, 3 NWKMON) C --- PRINT SOLUTION c WRITE (6,*) Y(1) c WRITE (6,*) Y(2) ID=0 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 nd=15 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,*)ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV C WRITE(6,99)N,T,HLAST,H,Y(1),Y(2),Y(3) C 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',3F15.8) RETURN END c C SUBROUTINE FVANDER( N, tn, Y, YDOT, dy, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), dy(N), RWORK(NR) 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) if(ires.eq.-1) then do 10 i=1,n dy(i)=-ydot(i) 10 continue else 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 do 167 i=1,15 dy(i)=dy(i)-ydot(i) 167 continue endif return end C SUBROUTINE JVANDER(N,tn,Y,YDOT,H,D,ML,MU,pd,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),pd(LDFY,N) 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) fac=h*d 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 do 192 i=1,15 do 192 j=1,15 pd(i,j)=pd(i,j)*fac 192 continue do 191 i=1,15 pd(i,i)=pd(i,i)-1.0D+0 191 continue return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT 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 (NRW=21+3*ND,NWKRES=4,NY2DIM=14) PARAMETER (MU=0,ML=0,NWKJAC=(2*ML+MU+1)*ND+2) PARAMETER (NWKMON=3,NIWJAC=ND,NYSAVE=ND) C --- DECLARATIONS DIMENSION Y(ND),YDOTI(ND),RWORK(NRW),INFORM(14+ND) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON),U(N) dimension error(nd),true(nd) CHARACTER*6 SNORM COMMON/TRANS/QQ,UZERO REAL*4 TARRAY(2) EXTERNAL PREPJB,SOLSB,MONITR,SBLEND,SPGEAR,STHETA,SPRMOM EXTERNAL FKS,JKS C --- DATA FOR THE PROBLEM QQ=0.025D0 c ------ FILE DE DONNEES ---------- write(6,3050) 3050 format(1x,'ks by sprint') 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 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 --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 C ---------- VAL INIT ------------ X=0.D0 XEND=100.D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=1 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 1, ND, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF C --- CALL OF SPRINT CALL SPRINT( ND, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FKS,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJB, WKJAC, NWKJAC, JACPVT, JKS, SOLSB,SPRMOM,WKMON, 3 NWKMON) 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 ID=0 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE (6,*) ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC IF (TARRAY(1).GT.1000.) STOP 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV c WRITE(6,99)N,T,HLAST,H,Y(1),Y(2) c 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',2F13.8) 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, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH) DIMENSION Y(ND), YDOT(ND), R(ND), RWORK(NR) DIMENSION U(N) COMMON/TRANS/QQ,UZERO IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE 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) R(2*J-1)=DIAG*Y(2*J-1) +QQ*J*U(2*J+2) - YDOT(2*J-1) R(2*J )=DIAG*Y(2*J ) -QQ*J*U(2*J+1) - YDOT(2*J ) enddo END IF end c-------------------------------------------------------------- SUBROUTINE JKS(ND,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION y(ND),DFY(LDFY,ND),YDOT(ND) COMMON/TRANS/QQ,UZERO FAC=H*D do J=1,NH-1 DIAG=FAC*(QQ*J)**2*(1.D0-(QQ*J)**2) DFY(1,2*J-1)=DIAG - 1.0D0 DFY(1,2*J )=DIAG - 1.0D0 enddo end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=6,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (NWKJAC=400,NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJF,SOLSF,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FVANDER,JVANDER,SPRMOM c ------ FILE DE DONNEES ---------- C --- LOOP FOR DIFFERENT TOLERANCES write(6,3500) 3500 format(1x,'B5 by sprint') 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 nd=6 C --- DIMENSION OF THE SYSTEM N=6 C --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 2 c c this specifies the mode of running. =2 is one step. = 1 is usual INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 C --- INITIAL VALUES X=0.0D0 do 791 i=1,6 y(i)=1.0D+0 791 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION XEND=20.0D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 2, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF 220 continue C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FVANDER,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJF, WKJAC, NWKJAC, JACPVT, JVANDER, SOLSF,SPRMOM,WKMON, 3 NWKMON) C --- PRINT SOLUTION c WRITE (6,*) Y(1) c WRITE (6,*) Y(2) c write(6,*) x ID=0 true(1)=dexp(-10.0D+0*x)*(cos(1000.0D+0*x)+sin(1000.0D+0*x)) true(2)=dexp(-10.0D+0*x)*(cos(1000.0D+0*x)-sin(1000.0D+0*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 345 ik=1,6 error(ik)=true(ik)-y(ik) sum=sum+((error(ik)/(rtol*dabs(y(ik))+atol))**2) 345 continue sum=dsqrt(sum/6.0D+0) if(sum.gt.summ) summ=sum c write(6,*) sum,summ if(x.lt.xend) GOTO 220 955 continue xout=x aa=cos(1000.0D+0*xout) bb=sin(1000.0D+0*xout) cc=dexp(-10.0D+0*xout) true(1)=cc*(aa+bb) true(2)=cc*(aa-bb) true(3)=dexp(-4.0D+0*xout) true(4)=dexp(-xout) true(5)=dexp(-0.5D+0*xout) true(6)=dexp(-0.1D+0*xout) sum=0.0D+0 do 270 k=1,6 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 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,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV C WRITE(6,99)N,T,HLAST,H,Y(1),Y(2),Y(3) C 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',3F15.8) RETURN END c C SUBROUTINE FVANDER( N, T, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE r(1)=-10.0D+0*y(1)+1000.0D+0*y(2) -ydot(1) r(2)=-1000.0D+0*y(1)-10.0D+0*y(2) -ydot(2) r(3)=-4.0D+0*y(3)-ydot(3) r(4)=-y(4)-ydot(4) r(5)=-0.5D+0*y(5)-ydot(5) r(6)=-0.1D+0*y(6)-ydot(6) END IF RETURN END C SUBROUTINE JVANDER(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) FAC=H*D DFY(1,1)= -fac*10.0D+0 -1.D0 DFY(1,2)=FAC*1000.0D+0 DFY(2,1)=FAC*(-1000.0D+0) DFY(2,2)=-10.0D+0*FAC -1.D0 dfy(3,3)=-4.0D+0*fac-1.0D+0 dfy(4,4)=-1.0D+0*fac-1.0D+0 dfy(5,5)=-0.5D+0*fac-1.0D+0 dfy(6,6)=-0.1D+0*fac-1.0D+0 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT 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 (NRW=21+3*ND,NWKRES=4,NY2DIM=14) PARAMETER (MU=0,ML=0,NWKJAC=(2*ML+MU+1)*ND+2) PARAMETER (NWKMON=3,NIWJAC=ND,NYSAVE=ND) C --- DECLARATIONS DIMENSION Y(ND),YDOTI(ND),RWORK(NRW),INFORM(14+ND) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON),U(N) dimension error(nd),true(nd) CHARACTER*6 SNORM COMMON/TRANS/QQ,UZERO REAL*4 TARRAY(2) EXTERNAL PREPJB,SOLSB,MONITR,SBLEND,SPGEAR,STHETA,SPRMOM EXTERNAL FKS,JKS open(6,file='sprintks.res') C --- DATA FOR THE PROBLEM QQ=0.025D0 c ------ FILE DE DONNEES ---------- write(6,3050) 3050 format(1x,'ks by sprint') 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 --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 C ---------- VAL INIT ------------ X=0.D0 XEND=100.D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=1 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 1, ND, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF C --- CALL OF SPRINT CALL SPRINT( ND, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FKS,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJB, WKJAC, NWKJAC, JACPVT, JKS, SOLSB,SPRMOM,WKMON, 3 NWKMON) 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 ID=0 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE (6,*) ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- TOLST=TOLST*TOLFC IF (TARRAY(1).GT.1000.) STOP 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV c WRITE(6,99)N,T,HLAST,H,Y(1),Y(2) c 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',2F13.8) 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, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH) DIMENSION Y(ND), YDOT(ND), R(ND), RWORK(NR) DIMENSION U(N) COMMON/TRANS/QQ,UZERO IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE 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) R(2*J-1)=DIAG*Y(2*J-1) +QQ*J*U(2*J+2) - YDOT(2*J-1) R(2*J )=DIAG*Y(2*J ) -QQ*J*U(2*J+1) - YDOT(2*J ) enddo END IF end c-------------------------------------------------------------- SUBROUTINE JKS(ND,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION y(ND),DFY(LDFY,ND),YDOT(ND) COMMON/TRANS/QQ,UZERO FAC=H*D do J=1,NH-1 DIAG=FAC*(QQ*J)**2*(1.D0-(QQ*J)**2) DFY(1,2*J-1)=DIAG - 1.0D0 DFY(1,2*J )=DIAG - 1.0D0 enddo end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=80,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (NWKJAC=NDIM*NDIM+2,NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJF,SOLSF,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- COMMON/NNNN/NCOM,NNCOM,NSQ,NQUATR,DELTAS EXTERNAL FTIGE,JTIGE c ------ FILE DE DONNEES ---------- open(6,file='sprintbeam.res') write(6,3050) 3050 format(1x,' sprint on beam') 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,21 12 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 100000 DO 14 I = 5,14+N INFORM(I) = 0 14 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ISET=0 C --------- INITIAL VALUES ------------- T=0.D0 DO 1 I=1,NN 1 Y(I)=0.D0 TEND=5.D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 2, NN, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF C --- CALL OF SPRINT CALL SPRINT( NN, T, TEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FTIGE,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJF, WKJAC, NWKJAC, JACPVT, JTIGE, SOLSF,SPRMON,WKMON, 3 NWKMON) 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 us=log10(summ) ut=log10(sumh) write(6,*) us,ut ID=0 C --- PRINT SOLUTION c WRITE(6,9921)(Y(I),I=1,NN) 9921 FORMAT(1X,F22.16) ID=0 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE (6,*) ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' IF(TARRAY(1).GT.300.)GOTO 66 C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE 66 CONTINUE STOP END C SUBROUTINE FTIGE( NN, T, TH, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION YDOT(NN), R(NN), RWORK(NR) DIMENSION TH(150),U(150),V(150),W(150) DIMENSION ALPHA(150),BETA(150),STH(150),CTH(150) COMMON/NNNN/N,NNCOM,NSQ,NQUATR,DELTAS IF(IRES .EQ. -1)THEN DO 10 I = 1,NN 10 R(I) = - YDOT(I) ELSE c----------- 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 R(I)=TH(N+I) -YDOT(I) 54 R(N+I)=U(I) -YDOT(N+I) END IF RETURN END SUBROUTINE JTIGE(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) WRITE(6,*)' OHO, FROM WHERE DO I COME ???' STOP END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=15,NRW=21+3*NDIM,NWKRES=4,NY2DIM=14) PARAMETER (NWKJAC=400,NWKMON=3,NIWJAC=NDIM,NYSAVE=NDIM) C --- DECLARATIONS DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension yend(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) common/const/pi EXTERNAL PREPJF,SOLSF,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FVANDER,JVANDER,SPRMOM c ------ FILE DE DONNEES ---------- C --- LOOP FOR DIFFERENT TOLERANCES open(6,file='sprbigring.res') write(6,3500) 3500 format(1x,'vdp by sprint') pi=4.0D+0*datan(1.0D+0) NTOLMN=2 NTOLMX=12 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=15 C --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 INFORM(3) = 1 INFORM(4) = 5000000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 C --- INITIAL VALUES X=0.0D0 do 151 ijk=1,15 y(ijk)=0.0D+0 151 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION XEND=1.0D-3 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 2, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF DO 20 I=1,1 C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FVANDER,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJF, WKJAC, NWKJAC, JACPVT, JVANDER, SOLSF,SPRMOM,WKMON, 3 NWKMON) C --- PRINT SOLUTION c WRITE (6,*) Y(1) c WRITE (6,*) Y(2) ID=0 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 nd=15 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 summ=summ*atol sumh=sumh*atol us=log10(summ) write(6,*) summ,sumh write(6,*) us,us WRITE(6,*)ID,ID,ID,ID,ID,ID,ID WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV C WRITE(6,99)N,T,HLAST,H,Y(1),Y(2),Y(3) C 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',3F15.8) RETURN END c C SUBROUTINE FVANDER( N, tn, Y, YDOT, dy, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), dy(N), RWORK(NR) 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) if(ires.eq.-1) then do 10 i=1,n dy(i)=-ydot(i) 10 continue else 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 do 167 i=1,15 dy(i)=dy(i)-ydot(i) 167 continue endif return end C SUBROUTINE JVANDER(N,tn,Y,YDOT,H,D,ML,MU,pd,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),pd(LDFY,N) 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) fac=h*d 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 do 192 i=1,15 do 192 j=1,15 pd(i,j)=pd(i,j)*fac 192 continue do 191 i=1,15 pd(i,i)=pd(i,i)-1.0D+0 191 continue return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=49999,NRW=21+3*NDIM,NWKRES=4,NY2DIM=20) PARAMETER (NWKMON=ndim,NIWJAC=NDIM,NYSAVE=NDIM) parameter(mu=1,ml=1,nwkjac=(2*ml+mu+1)*ndim+2) DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJB,SOLSB,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FVANDER,JVANDER,SPRMOM c ------ FILE DE DONNEES ---------- COMMON /const/pi,dx,ntop open (6,file='sprhard.hardnores50') ntop=50 pi=4.0D+0*DATAN(1.0D+0) dx=1.0d0/(ndim+1.0d0) write(6,*) ntop,ndim C --- LOOP FOR DIFFERENT TOLERANCES write(6,3500) 3500 format(1x,'hard parabolic by sprint') 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 nd=49999 C --- DIMENSION OF THE SYSTEM N=49999 C --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 c c this specifies the mode of running. =2 is one step. = 1 is usual INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 iset=0 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=1 C --- ENDPOINT OF INTEGRATION XEND=2.0D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 1, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF 220 continue C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FVANDER,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJB, WKJAC, NWKJAC, JACPVT, JVANDER, SOLSB,SPRMOM,WKMON, 3 NWKMON) C --- PRINT SOLUTION c WRITE (6,*) Y(1) c WRITE (6,*) Y(2) c write(6,*) x ID=0 t=x do 12000 j=1,n sumer=0.0D+0 do 12001 jj=1,ntop sumer=sumer+dexp(-pi*pi*t*jj*jj)*dsin(pi*jj*j*dx) 12001 continue true(j)=sumer 12000 continue sum=0.0D+0 do 270 k=1,49999 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 it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) c 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 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV C WRITE(6,99)N,T,HLAST,H,Y(1),Y(2),Y(3) C 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',3F15.8) RETURN END c C SUBROUTINE FVANDER( N, T, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) COMMON /const/pi,dx,ntop tn=t IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE 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 r(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 r(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 r(n)=(y(n-1)-2.0d0*y(n))/(dx*dx) + +y(n)**2-0.5D+0*y(n)*sumer-0.5D+0*sumer**2 do 446 i=1,n r(i)=r(i)-ydot(i) 446 continue endif return end C SUBROUTINE JVANDER(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) common /const/pi,dx,ntop fac=h*d 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*fac 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 dfy(2,j)=dfy(2,j)*fac-1.0D+0 142 continue do 144 j=1,n-1 dfy(3,j)=a1*fac 144 continue return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=1999,NRW=21+3*NDIM,NWKRES=4,NY2DIM=20) PARAMETER (NWKMON=ndim,NIWJAC=NDIM,NYSAVE=NDIM) parameter(mu=1,ml=1,nwkjac=(2*ml+mu+1)*ndim+2) DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJB,SOLSB,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FVANDER,JVANDER,SPRMOM c ------ FILE DE DONNEES ---------- COMMON /const/pi,dx,ntop,dif open (6,file='sprflood.easynores') ntop=250 pi=4.0D+0*DATAN(1.0D+0) dif=2.0D-5 c dif=4.0D+0/500000.0D+0 dx=4.0d0/(ndim+1.0d0) write(6,*) dif,ndim C --- LOOP FOR DIFFERENT TOLERANCES write(6,3500) 3500 format(1x,'hard parabolic by sprint') 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 nd=1999 C --- DIMENSION OF THE SYSTEM N=1999 C --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 1 c c this specifies the mode of running. =2 is one step. = 1 is usual INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 iset=0 C --- INITIAL VALUES X=0.0D0 do 40 j=1,n y(j)=dexp(-(-2.0d0+j*dx)**2) 40 continue C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=1 C --- ENDPOINT OF INTEGRATION XEND=4.0D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 1, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF 220 continue C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FVANDER,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJB, WKJAC, NWKJAC, JACPVT, JVANDER, SOLSB,SPRMOM,WKMON, 3 NWKMON) C --- PRINT SOLUTION c WRITE (6,*) Y(1) c WRITE (6,*) Y(2) c write(6,*) x ID=0 t=x tn=x err=0.0d0 do 9930 j=1,n err=err+((y(j)-(1.0d0/dsqrt(1.0d0+4.0d0*dif*tn))* & dexp(-((-2.0d0+j*dx-tn)**2)/(1.0d0+4.0d0*dif*tn)))**2) + /((rtol*dabs(y(j))+atol)**2) 9930 continue sum=dsqrt(err/n) if(sum.gt.summ) summ=sum c write(6,*) sum,summ if(x.lt.xend) GOTO 220 955 continue t=x tn=x c err=0.0d0 c do 3300 j=1,n err=err+((y(j)-(1.0d0/dsqrt(1.0d0+4.0d0*dif*tn))* & dexp(-((-2.0d0+j*dx-tn)**2)/(1.0d0+4.0d0*dif*tn)))**2) + /((rtol*dabs(y(j))+atol)**2) 3300 continue sum=dsqrt(err/dble(nd)) sumh=sumh+sum it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) c 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 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV C WRITE(6,99)N,T,HLAST,H,Y(1),Y(2),Y(3) C 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',3F15.8) RETURN END c C SUBROUTINE FVANDER( N, T, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) COMMON /const/pi,dx,ntop,dif tn=t IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE amult=2.0d0*dx prod=dx*dx 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)) r(1)=dif*(y0-2.d0*y(1)+y(2))/(dx*dx)-(y(2)-y0)/(2.d0*dx) rat=dif/prod do 3044 i=2,n-1 r(i)=rat*(y(i-1)-2.0d0*y(i)+y(i+1))- & (y(i+1)-y(i-1))/(amult) 3044 continue r(n)=dif*(y(n-1)-2.d0*y(n)+yN)/(dx*dx)-(yN-y(n-1))/(2.d0*dx) do 446 i=1,n r(i)=r(i)-ydot(i) 446 continue endif return end C SUBROUTINE JVANDER(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) common /const/pi,dx,ntop,dif fac=h*d 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*fac 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 dfy(2,j)=dfy(2,j)*fac-1.0D+0 142 continue do 144 j=1,n-1 dfy(3,j)=a1*fac 144 continue return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR SPRINT C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=49999,NRW=21+3*NDIM,NWKRES=4,NY2DIM=20) PARAMETER (NWKMON=ndim,NIWJAC=NDIM,NYSAVE=NDIM) parameter(mu=1,ml=1,nwkjac=(2*ml+mu+1)*ndim+2) DIMENSION Y(NDIM),YDOTI(NDIM),RWORK(NRW),INFORM(14+NDIM) DIMENSION WKRES(NWKRES),YSAVE(NYSAVE,NY2DIM),WKJAC(NWKJAC) DIMENSION JACPVT(NIWJAC),WKMON(NWKMON) dimension true(ndim),error(ndim) CHARACTER*6 SNORM REAL*4 TARRAY(2) EXTERNAL PREPJB,SOLSB,MONITR,SBLEND,SPGEAR,STHETA,SPRMON C -------- END PARAMETER LIST ---------- EXTERNAL FVANDER,JVANDER,SPRMOM c ------ FILE DE DONNEES ---------- COMMON /const/pi,dx,ntop open (6,file='sprhard.hardres50') ntop=50 pi=4.0D+0*DATAN(1.0D+0) dx=1.0d0/(ndim+1.0d0) write(6,*) ntop,ndim C --- LOOP FOR DIFFERENT TOLERANCES write(6,3500) 3500 format(1x,'hard parabolic by sprint') 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 nd=49999 C --- DIMENSION OF THE SYSTEM N=49999 C --- SET DEFAULT VALUES DO 18 I = 1,21 18 RWORK(I)=0.D0 INFORM(1) = 0 INFORM(2) = 2 c c this specifies the mode of running. =2 is one step. = 1 is usual INFORM(3) = 1 INFORM(4) = 100000 DO 19 I = 5,14+N INFORM(I) = 0 19 CONTINUE SNORM = 'AVERL2' C --- COMPUTE THE JACOBIAN ANALYTICALLY ISET=1 iset=0 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=1 C --- ENDPOINT OF INTEGRATION XEND=2.0D0 c CALL DTIME(TARRAY) it1=mclock() C --- CALL THE SETUP ROUTINES AND SPRINT IF THE SETUPS ARE O.K. CALL MATSET( 1, N, ML, MU, NWKJAC, NIWJAC, ISET) IF(ISET .LT. 0)THEN WRITE(6,92) 92 FORMAT(' SETUP FAILED IN MATSET') STOP END IF MAXORD = 0 CALL BLDSET( NY2DIM, IRET, MAXORD) IF(IRET .LT. 0)THEN WRITE(6,91) 91 FORMAT(' SETUP FAILED IN INTEGRATOR INITIALISATION') STOP END IF 220 continue C --- CALL OF SPRINT CALL SPRINT( N, X, XEND, Y, YDOTI, RWORK, NRW, RTOL, ATOL 1 ,ITOL, INFORM, SNORM, FVANDER,WKRES,NWKRES,SBLEND,YSAVE,NYSAVE 2 , PREPJB, WKJAC, NWKJAC, JACPVT, JVANDER, SOLSB,SPRMOM,WKMON, 3 NWKMON) C --- PRINT SOLUTION c WRITE (6,*) Y(1) c WRITE (6,*) Y(2) c write(6,*) x ID=0 t=x do 12000 j=1,n sumer=0.0D+0 do 12001 jj=1,ntop sumer=sumer+dexp(-pi*pi*t*jj*jj)*dsin(pi*jj*j*dx) 12001 continue true(j)=sumer 12000 continue sum=0.0D+0 do 270 k=1,49999 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 it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) c 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 write(6,*) '*** steps=',inform(6),' fns =',inform(7),'jacs=', + inform(8),'***' WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE SPRMOM( N, T, HLAST, H, Y, YDOT, YSAVE, NYH, R, ACOR, 1 RESWK, NRESWK, WKMON, NWKMON, IMON, INLN, HMIN, HMXI) C********************************************************************** C DUMMY EXAMPLE MONITOR FOR SPRINT . C********************************************************************** INTEGER N, NYH, NRESWK, NWKMON, IMON, INLN, ITRACE, IDEV DOUBLE PRECISION T, HLAST, H, Y(1), YDOT(1), YSAVE(NYH,1), R(1), 1 ACOR(1), RESWK(NRESWK), WKMON(NWKMON), HMIN, HMXI COMMON /SDEV2/ ITRACE, IDEV C WRITE(6,99)N,T,HLAST,H,Y(1),Y(2),Y(3) C 99 FORMAT(I6,' X=',F10.5,' HLAST,H=',2F13.9,' Y=',3F15.8) RETURN END c C SUBROUTINE FVANDER( N, T, Y, YDOT, R, IRES, RWORK, NR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), YDOT(N), R(N), RWORK(NR) COMMON /const/pi,dx,ntop tn=t IF(IRES .EQ. -1)THEN DO 10 I = 1,N 10 R(I) = - YDOT(I) ELSE 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 r(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 r(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 r(n)=(y(n-1)-2.0d0*y(n))/(dx*dx) + +y(n)**2-0.5D+0*y(n)*sumer-0.5D+0*sumer**2 do 446 i=1,n r(i)=r(i)-ydot(i) 446 continue endif return end C SUBROUTINE JVANDER(N,X,Y,YDOT,H,D,ML,MU,DFY,LDFY) C --- JACOBIAN OF VAN DER POL'S EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),YDOT(N),DFY(LDFY,N) common /const/pi,dx,ntop fac=h*d 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*fac 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 dfy(2,j)=dfy(2,j)*fac-1.0D+0 142 continue do 144 j=1,n-1 dfy(3,j)=a1*fac 144 continue return end