C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON THE VANDERPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * CCC include '/home/hairer/fortran/programme/dmebdf/dmebdf.f' IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR LSODE (FULL JACOBIAN) PARAMETER (ND=2,LWORK=(31+2*ND)*ND+1,LIWORK=ND) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd) + ,error(nd) REAL*4 TARRAY(2) COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,MAXSTEP c ------ FILE DE DONNEES ---------- write(6,2050) 2050 format(1x,'vdp by mebdf') 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) maxstep=100000 DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM N=2 C --- SET DEFAULT VALUES MF=21 INDEX=1 LOUT=6 C --- INITIAL VALUES X=0.0D0 Y(1)=2.0D0 Y(2)=0.0D0 XEND=11.0D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ccc ITOL=2 ITOL=2 MAXDER=7 C --- MAXIMAL NUMBER OF STEPS c IWORK(6)=10000 C --- ENDPOINT OF INTEGRATION C XEND=1.0D0 H=1.0D-6 CCC CALL DTIME(TARRAY) it1=mclock() XOUT=1.0D0 DO 20 I=1,11 220 CONTINUE C --- CALL OF THE SUBROUTINE CALL DRIVER(N,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,MAXDER,ITOL,RTOL,ATOL) IF (INDEX.EQ.1) THEN INDEX=0 GOTO 220 ELSE XOUT=XOUT+1.D0 END IF if(i.eq.1) then true(1)=-0.1863646254808130E+01 true(2)= 0.7535430865435460E+00 else if(i.eq.2) then true(1)= 0.1706167732170456E+01 true(2)= -0.8928097010248257E+00 else if(i.eq.3) then true(1)= -0.1510606936744095E+01 true(2)= 0.1178380000730945E+01 else if(i.eq.4) then true(1)= 0.1194414677695905E+01 true(2)= -0.2799585996540082E+01 else if(i.eq.5) then true(1)= 0.1890428596416747E+01 true(2)= -0.7345118680166940E+00 else if(i.eq.6) then true(1)=-0.1737716306805883E+01 true(2)= 0.8604008653025923E+00 else if(i.eq.7) then true(1)= 0.1551614645548223E+01 true(2)= -0.1102382892533321E+01 else if(i.eq.8) then true(1)=-0.1278631984330405E+01 true(2)= 0.2013890883009155E+01 else if(i.eq.9) then true(1)= -0.1916552949489830E+01 true(2)= 0.7169573003463228E+00 else if(i.eq.10) then true(1)= 0.1768163792391936E+01 true(2)=-0.8315276407898496E+00 else if(i.eq.11) then true(1)=-0.1590150544829062E+01 true(2)= 0.1040279389212485E+01 endif sum=0.0D+0 do 270 k=1,2 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 continue it2=mclock() tarray(1)=(it2-it1)/100.0D+0 write(6,*) x,y(1),y(2) WRITE(6,*)TARRAY(1),summ,sumh/11.0d+0 summ=summ*atol sumh=sumh*atol/(11.0D+0) write(6,*) summ,sumh us=log10(summ) uv=log10(sumh) write(6,*) us,uv WRITE(6,*)NFE,NJE,NSTEP,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE F(N,X,Y,DF) C --- RIGHT-HAND SIDE OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DF(N) EPS=1.D-6 DF(1)=Y(2) PROD=1.D0-Y(1)*Y(1) DF(2)=(PROD*Y(2)-Y(1))/EPS RETURN END C SUBROUTINE PDERV(X,Y,DFY,N,MEBAND) C --- JACOBIAN OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(MEBAND,N) EPS=1.D-6 DFY(1,1)=0.D0 DFY(1,2)=1.D0 DFY(2,1)=(-2.0D0*Y(1)*Y(2)-1.0D0)/EPS DFY(2,2)=(1.0D0-Y(1)**2)/EPS RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON THE VANDERPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * CCC include '/home/hairer/fortran/programme/dmebdf/dmebdf.f' IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR LSODE (FULL JACOBIAN) PARAMETER (ND=3,LWORK=(31+2*ND)*ND+1,LIWORK=ND) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd), + error(nd) REAL*4 TARRAY(2) COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,MAXSTEP c ------ FILE DE DONNEES ---------- write(6,2050) 2050 format(1x,'mebdf results for robertson') 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) maxstep=100000 DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM N=3 C --- SET DEFAULT VALUES MF=21 INDEX=1 LOUT=6 C --- INITIAL VALUES X=0.0D0 Y(1)=1.0D0 Y(2)=0.0D0 y(3)=0.0D+0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL*1.0D-6 ccc ITOL=2 ITOL=2 MAXDER=7 C --- MAXIMAL NUMBER OF STEPS c IWORK(6)=10000 C --- ENDPOINT OF INTEGRATION XEND=1.0d+11 H=1.0D-6 CCC CALL DTIME(TARRAY) it1=mclock() XOUT=1.0D+0 DO 20 I=1,12 220 CONTINUE C --- CALL OF THE SUBROUTINE CALL DRIVER(N,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,MAXDER,ITOL,RTOL,ATOL) IF (INDEX.EQ.1) THEN INDEX=0 GOTO 220 ELSE c if (mod(ntol,4).eq.1) write (6,*) xout,y(1),y(2) do 151 k=1,3 c write(6,*) y(k) 151 continue c write(6,*) h XOUT=XOUT*10.0D+0 END IF if(i.eq.1) then true(1)= 0.9664597373330035E+00 true(2)=0.3074626578578675E-04 true(3)=0.3350951640121071E-01 else if(i.eq.2) then true(1)=0.8413699238414729E+00 true(2)=0.1623390937990473E-04 true(3)=0.1586138422491472E+00 else if(i.eq.3) then true(1)=0.6172348823960878E+00 true(2)=0.6153591274639123E-05 true(3)=0.3827589640126376E+00 else if(i.eq.4) then true(1)=0.3368745306607069E+00 true(2)=0.2013702318261393E-05 true(3)=0.6631234556369748E+00 else if(i.eq.5) then true(1)=0.1073004285378040E+00 true(2)=0.4800166972571660E-06 true(3)=0.8926990914454987E+00 else if(i.eq.6) then true(1)=0.1786592114209946E-01 true(2)=0.7274751468436319E-07 true(3)=0.9821340061103859E+00 else if(i.eq.7) then true(1)= 0.2031483924973415E-02 true(2)=0.8142277783356159E-08 true(3)=0.9979685079327488E+00 else if(i.eq.8) then true(1)=0.2076093439016395E-03 true(2)=0.8306077485067610E-09 true(3)=0.9997923898254906E+00 else if(i.eq.9) then true(1)=0.2082417512179460E-04 true(2)=0.8329841429908955E-10 true(3)=0.9999791757415798E+00 else if(i.eq.10) then true(1)=0.2083229471647004E-05 true(2)=0.8332935037760723E-11 true(3)=0.9999979167621954E+00 else if(i.eq.11) then true(1)=0.2083328471883087E-06 true(2)=0.8333315602809495E-12 true(3)=0.9999997916663195E+00 else if(i.eq.12) then true(1)=0.2083340149701284E-07 true(2)=0.8333360770334744E-13 true(3)=0.9999999791665152E+00 end if sum=0.0d+0 do 270 k=1,3 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh/12.0D+0 summ=summ*atol sumh=sumh*atol/12.0D+0 write(6,*) SUMM,SUMH us=log10(summ) ut=log10(sumh) write(6,*) us,ut ID=0 WRITE(6,*)NFE,NJE,NSTEP,ID,ID,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END SUBROUTINE F(N,X,Y,dy) C --- RIGHT-HAND SIDE OF ROBERTSON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),dy(N) dy(1)=-0.04D0*Y(1)+1.D4*Y(2)*Y(3) dy(3)=3.D7*Y(2)*Y(2) dy(2)=-dy(1)-dy(3) RETURN END C SUBROUTINE PDERV(X,Y,DFY,N,MEBND) C --- JACOBIAN OF ROBERTSON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(MEBND,N) PROD1=1.0D4*Y(2) PROD2=1.0D4*Y(3) PROD3=6.0D7*Y(2) DFY(1,1)=-0.04D0 DFY(1,2)=PROD2 DFY(1,3)=PROD1 DFY(2,1)=0.04D0 DFY(2,2)=-PROD2-PROD3 DFY(2,3)=-PROD1 DFY(3,1)=0.D0 DFY(3,2)=PROD3 DFY(3,3)=0.D0 RETURN END c C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON THE VANDERPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * CCC include '/home/hairer/fortran/programme/dmebdf/dmebdf.f' IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR LSODE (FULL JACOBIAN) PARAMETER (ND=8,LWORK=(31+2*ND)*ND+1,LIWORK=ND) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd), + error(nd) REAL*4 TARRAY(2) COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,MAXSTEP c ------ FILE DE DONNEES ---------- write(6,2050) 2050 format(1x,'mebdf results for hires') 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) maxstep=100000 DO 30 NTOL=1,NRLOOP sumh=0.0D+0 summ=0.0D+0 C --- DIMENSION OF THE SYSTEM N=8 C --- SET DEFAULT VALUES MF=21 INDEX=1 LOUT=6 C --- INITIAL VALUES X=0.0D0 Y(1)=1.0D0 Y(2)=0.0D0 y(3)=0.0D+0 y(4)=0.0D+0 y(5)=0.0D+0 y(6)=0.0D+0 y(7)=0.0D+0 y(8)=0.0057d+0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL*1.0D-4 atol=rtol ccc ITOL=2 ITOL=2 MAXDER=7 C --- MAXIMAL NUMBER OF STEPS c IWORK(6)=10000 C --- ENDPOINT OF INTEGRATION XEND=421.8122d+0 H=1.0D-6 CCC CALL DTIME(TARRAY) it1=mclock() XOUT=xend-100.0D+0 DO 20 I=1,2 220 CONTINUE C --- CALL OF THE SUBROUTINE CALL DRIVER(N,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,MAXDER,ITOL,RTOL,ATOL) IF (INDEX.EQ.1) THEN INDEX=0 GOTO 220 ELSE XOUT=XOUT+100.D0 END IF if(i.eq.1) then true(1)=0.000737131257332567 true(2)=0.000144248572631618 true(3)=0.000058887297409676 true(4)=0.001175651343283149 true(5)=0.002386356198831330 true(6)=0.006238968252742796 true(7)=0.002849998395185769 true(8)=0.002850001604814231 else if(i.eq.2) then true(1)=0.000670305503581864 true(2)=0.000130996846986347 true(3)=0.000046862231597733 true(4)=0.001044668020551705 true(5)=0.000594883830951485 true(6)=0.001399628833942774 true(7)=0.001014492757718480 true(8)=0.004685507242281520 end if sum=0.0D+0 do 270 k=1,8 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh/2.0D+0 summ=summ*atol sumh=sumh*atol/2.0D+0 us=log10(summ) ut=log10(sumh) write(6,*) us,ut write(6,*) summ,sumh WRITE(6,*)NFE,NJE,NSTEP,ID,ID,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE F (N, X, Y, DY) INTEGER N DOUBLE PRECISION X, Y, DY DIMENSION Y (8), DY (8) DY (1) = -1.71D0*Y(1) + 0.43D0*Y(2) + 8.32D0*Y(3) + 0.0007D0 DY (2) = 1.71D0*Y(1) - 8.75D0*Y(2) DY (3) = -10.03D0*Y(3) + 0.43D0*Y(4) + 0.035D0*Y(5) DY (4) = 8.32D0*Y(2) + 1.71D0*Y(3) - 1.12D0*Y(4) DY (5) = -1.745D0*Y(5) + 0.43D0*Y(6) + 0.43D0*Y(7) DY (6) = -280.0D0*Y(6)*Y(8) + 0.69D0*Y(4) + 1.71D0*Y(5) - & 0.43D0*Y(6) + 0.69D0*Y(7) DY (7) = 280.0D0*Y(6)*Y(8) - 1.81D0*Y(7) DY (8) = -DY (7) RETURN END SUBROUTINE pderv(X,Y,DFY,N,mebnd) C -------- JACOBIAN FOR HIRES PROBLEM -------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(MEBND,N) C------ METTRE A ZERO ------- DO 1 I=1,N DO 1 J=1,N 1 DFY(I,J)=0.D0 C DFY(1,1)= -1.71D0 DFY(1,2)= 0.43D0 DFY(1,3)= + 8.32D0 C DFY(2,1)= 1.71D0 DFY(2,2)= - 8.75D0 C DFY(3,3)= -10.03D0 DFY(3,4)= 0.43D0 DFY(3,5)= + 0.035D0 C DFY(4,2)= 8.32D0 DFY(4,3)= + 1.71D0 DFY(4,4)= - 1.12D0 C DFY(5,5)= -1.745D0 DFY(5,6)= + 0.43D0 DFY(5,7)= + 0.43D0 C DFY(6,4)= + 0.69D0 DFY(6,5)= + 1.71D0 DFY(6,6)= - 0.43D0 -280.0D0*Y(8) DFY(6,7)= + 0.69D0 DFY(6,8)= -280.0D0*Y(6) C DFY(7,6)= 280.0D0*Y(8) DFY(7,7)= - 1.81D0 DFY(7,8)= 280.0D0*Y(6) C DFY(8,6)= - 280.0D0*Y(8) DFY(8,7)= 1.81D0 DFY(8,8)= - 280.0D0*Y(6) RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON THE VANDERPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * CCC include '/home/hairer/fortran/programme/dmebdf/dmebdf.f' IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR LSODE (FULL JACOBIAN) PARAMETER (ND=3,LWORK=(31+2*ND)*ND+1,LIWORK=ND) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd), + error(nd) REAL*4 TARRAY(2) COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep c ------ FILE DE DONNEES ---------- write(6,2050) 2050 format(1x,'mebdf results for orego') 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) maxstep=100000 DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM N=3 C --- SET DEFAULT VALUES MF=21 INDEX=1 LOUT=6 C --- INITIAL VALUES X=0.0D0 Y(1)=1.0D0 Y(2)=2.0D0 y(3)=3.0D+0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL*1.0D-6 ccc ITOL=2 ITOL=2 itol=1 MAXDER=7 C --- MAXIMAL NUMBER OF STEPS c IWORK(6)=10000 C --- ENDPOINT OF INTEGRATION XEND=360.0D+0 H=1.0D-6 CCC CALL DTIME(TARRAY) it1=mclock() XOUT=30.0D+0 DO 20 I=1,12 220 CONTINUE C --- CALL OF THE SUBROUTINE CALL DRIVER(N,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,MAXDER,ITOL,RTOL,ATOL) IF (INDEX.EQ.1) THEN INDEX=0 GOTO 220 ELSE do 151 k=1,3 c write(6,*) y(k) 151 continue c write(6,*) h XOUT=XOUT+30.0D+0 END IF if(i.eq.1) then true(1)=0.1000661467180497E+01 true(2)=0.1512778937348249E+04 true(3)=0.1035854312767229E+05 else if(i.eq.2) then true(1)=0.1000874625199626E+01 true(2)=0.1144336972384497E+04 true(3)=0.8372149966624639E+02 else if(i.eq.3) then true(1)=0.1001890368438751E+01 true(2)=0.5299926232295553E+03 true(3)=0.1662279579042420E+01 else if(i.eq.4) then true(1)=0.1004118022612645E+01 true(2)=0.2438326079910346E+03 true(3)=0.1008822224048647E+01 else if(i.eq.5) then true(1)=0.1008995416634061E+01 true(2)=0.1121664388662539E+03 true(3)=0.1007783229065319E+01 else if(i.eq.6) then true(1)=0.1019763472537298E+01 true(2)=0.5159761322947535E+02 true(3)=0.1016985778956374E+01 else if(i.eq.7) then true(1)= 0.1043985088527474E+01 true(2)=0.2373442027531524E+02 true(3)=0.1037691843544522E+01 else if(i.eq.8) then true(1)=0.1100849071667922E+01 true(2)=0.1091533805469020E+02 true(3)=0.1085831969810860E+01 else if(i.eq.9) then true(1)=0.1249102130020572E+01 true(2)=0.5013945178605446E+01 true(3)=0.1208326626237875E+01 else if(i.eq.10) then true(1)=0.1779724751937019E+01 true(2)=0.2281852385542403E+01 true(3)=0.1613754023671725E+01 else if(i.eq.11) then true(1)=0.1000889326903503E+01 true(2)=0.1125438585746596E+04 true(3)=0.1641049483777168E+05 else if(i.eq.12) then true(1)=0.1000814870318523E+01 true(2)=0.1228178521549889E+04 true(3)=0.1320554942846513E+03 endif sum=0.0D+0 do 270 k=1,3 error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh/12.0D+0 summ=summ*atol sumh=sumh*atol/12.0D+0 write(6,*) summ,sumh ID=0 WRITE(6,*)NFE,NJE,NSTEP,ID,ID,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE F(N,X,Y,DY) C --- RIGHT-HAND SIDE OF OREGON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DY(N) DY(1)=77.27D0*(Y(2)+Y(1)*(1.D0-8.375D-6*Y(1)-Y(2))) DY(2)=(Y(3)-(1.D0+Y(1))*Y(2))/77.27D0 DY(3)=0.161D0*(Y(1)-Y(3)) RETURN END C SUBROUTINE PDERV(X,Y,DFY,N,MEBND) C --- JACOBIAN OF OREGON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(MEBND,N) DFY(1,1)=77.27D0*(1.D0-2.D0*8.375D-6*Y(1)-Y(2)) DFY(1,2)=77.27D0*(1.D0-Y(1)) DFY(1,3)=0.D0 DFY(2,1)=-Y(2)/77.27D0 DFY(2,2)=-(1.D0+Y(1))/77.27D0 DFY(2,3)=1.D0/77.27D0 DFY(3,1)=.161D0 DFY(3,2)=.0D0 DFY(3,3)=-.161D0 RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON E5 PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * CCC include '/home/hairer/fortran/programme/dmebdf/dmebdf.f' IMPLICIT double precision (A-H,O-Z) C --- PARAMETERS FOR LSODE (FULL JACOBIAN) PARAMETER (ND=4,LWORK=(31+2*ND)*ND+1,LIWORK=ND) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd), + error(nd) double precision TARRAY(2) COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep c ------ FILE DE DONNEES ---------- write(6,2050) 2050 format(1x,'mebdf results for E5') 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) maxstep=100000 DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM N=4 C --- SET DEFAULT VALUES MF=21 INDEX=1 LOUT=6 C --- INITIAL VALUES X=0.0D0 Y(1)=1.76D-3 Y(2)=0.0D0 y(3)=0.0D+0 y(4)=0.0D+0 C --- REQUIRED TOLERANCE RTOL=TOLST atol=1.7d-24 ccc ITOL=2 ITOL=2 MAXDER=7 C --- MAXIMAL NUMBER OF STEPS c IWORK(6)=10000 C --- ENDPOINT OF INTEGRATION XEND=10000000000000.0D+0 H=1.0D-6 CCC CALL DTIME(TARRAY) it1=mclock() XOUT=10.0D+0 DO 20 I=1,7 220 CONTINUE C --- CALL OF THE SUBROUTINE CALL DRIVER(N,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,MAXDER,ITOL,RTOL,ATOL) IF (INDEX.EQ.1) THEN INDEX=0 GOTO 220 ELSE XOUT=XOUT*100.0D+0 END IF 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)) c write(6,*) true(k),y(k),error(k) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh/7.0D+0 summ=summ*rtol sumh=sumh*rtol/7.0D+0 write(6,*) summ,sumh us=log10(summ) ut=log10(sumh) write(6,*) us,ut ID=0 WRITE(6,*)NFE,NJE,NSTEP,ID,ID,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END SUBROUTINE F(N,X,Y,DF) C --- RIGHT-HAND SIDE OF ROBERTSON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DF(N) PROD1=7.89D-10*Y(1) PROD2=1.1D7*Y(1)*Y(3) PROD3=1.13D9*Y(2)*Y(3) PROD4=1.13D3*Y(4) DF(1)=-PROD1-PROD2 DF(2)=PROD1-PROD3 DF(4)=PROD2-PROD4 DF(3)=DF(2)-DF(4) RETURN END C SUBROUTINE PDERV(X,Y,DFY,N,MEBND) C --- JACOBIAN OF ROBERTSON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(MEBND,N) A=7.89D-10 B=1.1D7 CM=1.13D9 C=1.13D3 DFY(1,1)=-A-B*Y(3) DFY(1,2)=0.D0 DFY(1,3)=-B*Y(1) DFY(1,4)=0.D0 DFY(2,1)=A DFY(2,2)=-CM*Y(3) DFY(2,3)=-CM*Y(2) DFY(2,4)=0.D0 DFY(3,1)=A-B*Y(3) DFY(3,2)=-CM*Y(3) DFY(3,3)=-B*Y(1)-CM*Y(2) DFY(3,4)=C DFY(4,1)=B*Y(3) DFY(4,2)=0.D0 DFY(4,3)=B*Y(1) DFY(4,4)=-C RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON PLATE PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * CCC include '/home/hairer/fortran/programme/dmebdf/dmebdf.f' IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR LSODE (FULL JACOBIAN) PARAMETER(MX=8,MY=5,MACHS1=2,MACHS2=4) PARAMETER(ND=2*MX*MY) PARAMETER (LWORK=(31+2*ND)*ND+1,LIWORK=ND) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd), + error(nd) REAL*4 TARRAY(2) COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, + OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT maxstep=100000 NX=MX NY=MY NACHS1=MACHS1 NACHS2=MACHS2 NXM1=NX-1 NYm1=NY-1 NDEMI=NX*NY OMEGA=1000.0D+0 STIFFN=100.0D+0 WEIGHT=200.0D+0 DENOM=NX+1 DELX=2.0D+0/DENOM USH4=1.0D+0/(DELX**4) FAC=STIFFN*USH4 c ------ FILE DE DONNEES ---------- write(6,2050) 2050 format(1x,'mebdf results for PLATE') 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=ND C --- SET DEFAULT VALUES MF=21 INDEX=1 LOUT=6 C --- INITIAL VALUES X=0.0D0 do 1 i=1,N y(i)=0.0D+0 1 continue C --- REQUIRED TOLERANCE RTOL=TOLST atol=RTOL ccc ITOL=2 ITOL=2 MAXDER=7 C --- MAXIMAL NUMBER OF STEPS c IWORK(6)=10000 C --- ENDPOINT OF INTEGRATION XEND=7.0D+0 H=1.0D-6 CCC CALL DTIME(TARRAY) it1=mclock() XOUT=7.0D+0 220 CONTINUE C --- CALL OF THE SUBROUTINE CALL DRIVER(N,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,MAXDER,ITOL,RTOL,ATOL) IF (INDEX.EQ.1) THEN INDEX=0 GOTO 220 ELSE 9921 format(1x,f22.16) END IF 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 us=log10(summ) ut=log10(sumh) write(6,*) us,ut ID=0 WRITE(6,*)NFE,NJE,NSTEP,ID,ID,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE F (N, X, Y, DF) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N), DF(N) COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT c ################ print c write(6,*)' dans fplate, n,x=',n,x c############## C -------- LA BOUCLE ------- DO 1 I=1,NX DO 1 J=1,NY K=I+NX*(J-1) C -------- DERIVEE DEUXIEME ---- DF(K)=Y(K+NDEMI) C ------ POINT CENTRAL --- UC=16.D0*Y(K) IF(I.GT.1)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K-1) END IF IF(I.LT.NX)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K+1) END IF IF(J.GT.1)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K-NX) END IF IF(J.LT.NY)THEN UC=UC+Y(K) UC=UC-8.D0*Y(K+NX) END IF IF(I.GT.1 .AND.J.GT.1 )UC=UC+2.D0*Y(K-NX-1) IF(I.LT.NX.AND.J.GT.1 )UC=UC+2.D0*Y(K-NX+1) IF(I.GT.1 .AND.J.LT.NY)UC=UC+2.D0*Y(K+NX-1) IF(I.LT.NX.AND.J.LT.NY)UC=UC+2.D0*Y(K+NX+1) IF(I.GT.2)UC=UC+Y(K-2) IF(I.LT.NXM1)UC=UC+Y(K+2) IF(J.GT.2)UC=UC+Y(K-2*NX) IF(J.LT.NYM1)UC=UC+Y(K+2*NX) IF(J.EQ.NACHS1.OR.J.EQ.NACHS2)THEN XI=I*DELX FORCE=EXP(-5.D0*(X-XI-2.D0)**2)+EXP(-5.D0*(X-XI-5.D0)**2) ELSE FORCE=0.D0 END IF DF(K+NDEMI)=-OMEGA*Y(K+NDEMI)-FAC*UC+FORCE*WEIGHT 1 CONTINUE RETURN END SUBROUTINE PDERV(X,Y,DFY,N,MEBND) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(MEBND,N) COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY, & OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT c ################ print c write(6,*)' dans jac , n,x=',n,x c############## C------ METTRE A ZERO ------- DO 1 I=1,N DO 1 J=1,N 1 DFY(I,J)=0.D0 C -------- LA BOUCLE ------- DO 2 I=1,NX DO 2 J=1,NY K=I+NX*(J-1) C -------- DERIVEE DEUXIEME ---- DFY(K,K+NDEMI)=1.D0 C ------ POINT CENTRAL --- DFY(K+NDEMI,K)=-FAC*16.D0 IF(I.GT.1)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K-1)=FAC*8.D0 END IF IF(I.LT.NX)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K+1)=FAC*8.D0 END IF IF(J.GT.1)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K-NX)=FAC*8.D0 END IF IF(J.LT.NY)THEN DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC DFY(K+NDEMI,K+NX)=FAC*8.D0 END IF IF(I.GT.1 .AND.J.GT.1 )DFY(K+NDEMI,K-NX-1)=-FAC*2.D0 IF(I.LT.NX.AND.J.GT.1 )DFY(K+NDEMI,K-NX+1)=-FAC*2.D0 IF(I.GT.1 .AND.J.LT.NY)DFY(K+NDEMI,K+NX-1)=-FAC*2.D0 IF(I.LT.NX.AND.J.LT.NY)DFY(K+NDEMI,K+NX+1)=-FAC*2.D0 IF(I.GT.2)DFY(K+NDEMI,K-2)=-FAC IF(I.LT.NXM1)DFY(K+NDEMI,K+2)=-FAC IF(J.GT.2)DFY(K+NDEMI,K-2*NX)=-FAC IF(J.LT.NYM1)DFY(K+NDEMI,K+2*NX)=-FAC DFY(K+NDEMI,K+NDEMI)= -OMEGA 2 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON CUSP PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR MEBDF (FULL JACOBIAN) PARAMETER (ND=96,LWORK=50000,LIWORK=500) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd), + error(nd) COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,MAXSTEP C -------- END PARAMETER LIST -------- common/type/linear logical linear REAL*4 TARRAY(2) COMMON/NERVES/NNERV COMMON/DIFFCOEF/DIFFUS EXTERNAL F,pderv,SOLOUT linear=.false. c ------ FILE DE DONNEES ---------- write(6,2045) 2045 format(1x,'results for mebdf on cusp') C --- LOOP FOR DIFFERENT TOLERANCES maxstep=100000 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 -------- INITIAL CONSTANTES -------- NNERV=32 DIFFUS=1.D0*NNERV*NNERV/144.D0 N=3*NNERV C --- COMPUTE THE JACOBIAN SWITCH IJAC=0 C --- JACOBIAN IS A FULL MATRIX ML=3 MU=3 mbnd(1)=ml mbnd(2)=mu mbnd(3)=ml+mu+1 mbnd(4)=mbnd(3)+mu C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM IMAS=0 C --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION IOUT=0 C ---------- VAL INIT ------------ X=0.D0 XEND=1.1D0 ANERV=NNERV DEL=2.D0*3.14159265358979324D0/ANERV DO 25 INERV=1,NNERV Y(3*INERV-2)=0.D0 Y(3*INERV-1)=-2.D0*COS(INERV*DEL) 25 Y(3*INERV)=2.D0*SIN(INERV*DEL) C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ITOL=2 C --- INITIAL STEP SIZE H=1.0D-6 mf=24 index=1 lout=6 maxder=7 XEND=1.1D0 xout=xend it1=mclock() 220 continue C --- CALL OF THE SUBROUTINE RADAU5 (OR SDIRK4) call driver(n,x,h,y,xout,xend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,maxder,itol,rtol,atol) C --- PRINT SOLUTION if(index.eq.1) then index=0 goto 220 endif 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 us=log10(summ) ut=log10(sumh) write(6,*) us,ut WRITE (6,*)(IWORK(J),J=14,20) WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' WRITE (6,91) (IWORK(J),J=14,20) 91 FORMAT(' fcn=',I5,' jac=',I4,' step=',I4, & ' accpt=',I4,' rejct=',I3,' dec=',I4, & ' sol=',I5) write (8,*) y(1),y(n/2),y(n) C -------- NEW TOLERANCE --- IF (TARRAY(1).GT.500.) STOP TOLST=TOLST*TOLFC 30 CONTINUE STOP END C SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS C --- BY USING "CONTR5" (OR "CONTS4") IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N) COMMON /INTERN/XOUT COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL write (6,*) 'stat',NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL WRITE (6,99) X,Y(1),Y(2),NR-1 99 FORMAT(1X,'X =',F8.4,' Y =',2F18.7,' NSTEP =',I4) RETURN END SUBROUTINE F(N,T,Y,DF) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DF(N) COMMON/NERVES/NNERV COMMON/DIFFCOEF/DIFFUS c----------- DO 25 INERV=1,NNERV X=Y(3*INERV-2) A=Y(3*INERV-1) B=Y(3*INERV) IF(INERV.EQ.1)THEN XRIGHT=Y(3*NNERV-2) ARIGHT=Y(3*NNERV-1) BRIGHT=Y(3*NNERV) ELSE XRIGHT=Y(3*INERV-5) ARIGHT=Y(3*INERV-4) BRIGHT=Y(3*INERV-3) END IF IF(INERV.EQ.NNERV)THEN XLEFT=Y(1) ALEFT=Y(2) BLEFT=Y(3) ELSE XLEFT=Y(3*INERV+1) ALEFT=Y(3*INERV+2) BLEFT=Y(3*INERV+3) END IF XDOT=-10000.D0*(B+X*(A+X*X)) U=(X-0.7D0)*(X-1.3D0) V=U/(U+0.1D0) ADOT=B+0.07D0*V BDOT=(1.D0*(1.D0-A*A)*B-A)-0.4D0*X+0.035D0*V DF(3*INERV-2)=XDOT+DIFFUS*(XLEFT-2.D0*X+XRIGHT) DF(3*INERV-1)=ADOT+DIFFUS*(ALEFT-2.D0*A+ARIGHT) DF(3*INERV) =BDOT+DIFFUS*(BLEFT-2.D0*B+BRIGHT) 25 CONTINUE RETURN END subroutine pderv(x,y,df,n,mebnd) return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON BEAM PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * CCC include '/home/hairer/fortran/programme/dmebdf/dmebdf.f' IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR LSODE (FULL JACOBIAN) PARAMETER (ND=80,LWORK=(31+2*ND)*ND+1,LIWORK=ND) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd), + error(nd) REAL*4 TARRAY(2) COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /NNNN/NCOM,NNCOM,NSQ,NQUATR,DELTAS c ------ FILE DE DONNEES ---------- open(6,file='mebdfbeam.res3') maxstep=100000 write(6,2050) 2050 format(1x,'mebdf results for BEAM') C --- LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=4 ntoldf=1 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM N=40 C --- SET DEFAULT VALUES MF=22 INDEX=1 LOUT=6 NN=2*N NCOM=N NSQ=N*N NQUATR=NSQ*NSQ NNCOM=NN AN=N DELTAS=1.0D+0/AN C --- INITIAL VALUES X=0.0D0 do 1 i=1,NN y(i)=0.0D+0 1 continue C --- REQUIRED TOLERANCE RTOL=TOLST atol=RTOL ccc ITOL=2 ITOL=2 MAXDER=3 write(6,3070) maxder 3070 format(1x,'maximum order is', i8) C --- MAXIMAL NUMBER OF STEPS c IWORK(6)=10000 C --- ENDPOINT OF INTEGRATION XEND=5.0D+0 H=1.0D-6 CCC CALL DTIME(TARRAY) it1=mclock() XOUT=5.0D+0 220 CONTINUE C --- CALL OF THE SUBROUTINE CALL DRIVER(NN,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,MAXDER,ITOL,RTOL,ATOL) IF (INDEX.EQ.1) THEN INDEX=0 GOTO 220 ELSE 9921 format(1x,f22.16) END IF 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 WRITE(6,*)NFE,NJE,NSTEP,ID,ID,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END SUBROUTINE F(NN,T,TH,DF) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DF(NN),TH(150),U(150),V(150),W(150) DIMENSION ALPHA(150),BETA(150),STH(150),CTH(150) COMMON/NNNN/N,NNCOM,NSQ,NQUATR,DELTAS C ----- CALCUL DES TH(I) ET DES SIN ET COS ------------- DO 22 I=2,N THDIFF=TH(I)-TH(I-1) STH(I)=DSIN(THDIFF) 22 CTH(I)=DCOS(THDIFF) C -------- CALCUL DU COTE DROIT DU SYSTEME LINEAIRE ----- IF(T.GT.3.14159265358979324D0)THEN C --------- CASE T GREATER PI ------------ C ---------- I=1 ------------ TERM1=(-3.D0*TH(1)+TH(2))*NQUATR V(1)=TERM1 C -------- I=2,..,N-1 ----------- DO 32 I=2,N-1 TERM1=(TH(I-1)-2.D0*TH(I)+TH(I+1))*NQUATR 32 V(I)=TERM1 C ----------- I=N ------------- TERM1=(TH(N-1)-TH(N))*NQUATR V(N)=TERM1 ELSE C --------- CASE T LESS EQUAL PI ------------ FABS=1.5D0*DSIN(T)*DSIN(T) FX=-FABS FY= FABS C ---------- I=1 ------------ TERM1=(-3.D0*TH(1)+TH(2))*NQUATR TERM2=NSQ*(FY*DCOS(TH(1))-FX*DSIN(TH(1))) V(1)=TERM1+TERM2 C -------- I=2,..,N-1 ----------- DO 34 I=2,N-1 TERM1=(TH(I-1)-2.D0*TH(I)+TH(I+1))*NQUATR TERM2=NSQ*(FY*DCOS(TH(I))-FX*DSIN(TH(I))) 34 V(I)=TERM1+TERM2 C ----------- I=N ------------- TERM1=(TH(N-1)-TH(N))*NQUATR TERM2=NSQ*(FY*DCOS(TH(N))-FX*DSIN(TH(N))) V(N)=TERM1+TERM2 END IF C -------- COMPUTE PRODUCT DV=W ------------ W(1)=STH(2)*V(2) DO 43 I=2,N-1 43 W(I)=-STH(I)*V(I-1)+STH(I+1)*V(I+1) W(N)=-STH(N)*V(N-1) C -------- TERME 3 ----------------- DO 435 I=1,N 435 W(I)=W(I)+TH(N+I)*TH(N+I) C ------- SOLVE SYSTEM CW=W --------- ALPHA(1)=1.D0 DO 44 I=2,N ALPHA(I)=2.D0 44 BETA(I-1)=-CTH(I) ALPHA(N)=3.D0 DO 45 I=N-1,1,-1 Q=BETA(I)/ALPHA(I+1) W(I)=W(I)-W(I+1)*Q 45 ALPHA(I)=ALPHA(I)-BETA(I)*Q W(1)=W(1)/ALPHA(1) DO 46 I=2,N 46 W(I)=(W(I)-BETA(I-1)*W(I-1))/ALPHA(I) C -------- COMPUTE U=CV+DW --------- U(1)=V(1)-CTH(2)*V(2)+STH(2)*W(2) DO 47 I=2,N-1 47 U(I)=2.D0*V(I)-CTH(I)*V(I-1)-CTH(I+1)*V(I+1) & -STH(I)*W(I-1)+STH(I+1)*W(I+1) U(N)=3.D0*V(N)-CTH(N)*V(N-1)-STH(N)*W(N-1) C -------- PUT DERIVATIVES IN RIGHT PLACE ------------- DO 54 I=1,N DF(I)=TH(N+I) 54 DF(N+I)=U(I) RETURN END subroutine pderv(x,y,dfy,n,mebnd) dimension y(n),dfy(mebnd,n) return end C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON BRUSS PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * CCC include '/home/hairer/fortran/programme/dmebdf/dmebdf.f' IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR LSODE (FULL JACOBIAN) PARAMETER (ND=1000,LWORK=(31+2*ND)*ND+1,LIWORK=ND) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd), + error(nd) REAL*4 TARRAY(2) COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /param/n,n2,gamma,gamma2 c ------ FILE DE DONNEES ---------- write(6,2050) 2050 format(1x,'mebdf results for BRUSS') pi=4.0D+0*datan(1.0D+0) n=500 n2=2*n maxstep=100000 usdelq=(dble(n+1))**2 gamma=0.02D+0*usdelq gamma2=2.0D+0*gamma 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 C --- SET DEFAULT VALUES MF=23 INDEX=1 LOUT=6 C --- INITIAL VALUES X=0.0D0 do 1 i=1,N anp1=n+1 xi=i/ANP1 Y(2*I)=3.0d+0 Y(2*I-1)=1.0d+0+0.5d+0*DSIN(2.0d+0*PI*XI) 1 CONTINUE C --- REQUIRED TOLERANCE RTOL=TOLST atol=RTOL ccc ITOL=2 ITOL=2 MAXDER=7 C --- MAXIMAL NUMBER OF STEPS c IWORK(6)=10000 C --- ENDPOINT OF INTEGRATION XEND=10.0D+0 H=1.0D-6 CCC CALL DTIME(TARRAY) it1=mclock() XOUT=10.0D+0 220 CONTINUE mbnd(1)=2 mbnd(2)=2 mbnd(3)=5 mbnd(4)=7 C --- CALL OF THE SUBROUTINE CALL DRIVER(N2,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,MAXDER,ITOL,RTOL,ATOL) IF (INDEX.EQ.1) THEN INDEX=0 GOTO 220 END IF 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 errm=0.0D+0 ihold=0 do 270 k=1,995,7 error(k)=dabs(true(k)-y(k)) if(error(k).gt.errm) ihold=k if(error(k).gt.errm) errm=error(k) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue write(6,*) true(ihold),y(ihold),error(ihold),ihold 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 us=log10(summ) ut=log10(sumh) write(6,*) us,ut ID=0 WRITE(6,*)NFE,NJE,NSTEP,ID,ID,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END SUBROUTINE F(NNN,X,Y,DF) C --- RIGHT-HAND SIDE OF BRUSPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(NNN),DF(NNN) COMMON/PARAM/N,N2,GAMMA,GAMMA2 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 DF(IU)=1.D0+PROD-4.D0*UI+GAMMA*(UIM-2.D0*UI+UIP) DF(IV)=3.D0*UI-PROD+GAMMA*(VIM-2.D0*VI+VIP) DO 5 I=2,N-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 DF(IU)=1.D0+PROD-4.D0*UI+GAMMA*(UIM-2.D0*UI+UIP) DF(IV)=3.D0*UI-PROD+GAMMA*(VIM-2.D0*VI+VIP) 5 CONTINUE I=N 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 DF(IU)=1.D0+PROD-4.D0*UI+GAMMA*(UIM-2.D0*UI+UIP) DF(IV)=3.D0*UI-PROD+GAMMA*(VIM-2.D0*VI+VIP) RETURN END C SUBROUTINE pderv(X,Y,DFY,NNN,mebnd) C --- JACOBIAN OF BRUSPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(NNN),DFY(mebnd,NNN) COMMON/PARAM/N,N2,GAMMA,GAMMA2 DO 1 I=1,N IU=2*I-1 IV=2*I UI=Y(IU) VI=Y(IV) UIVI=UI*VI UI2=UI*UI DFY(3,IU)=2.D0*UIVI-4.D0-GAMMA2 DFY(2,IV)=UI2 DFY(4,IU)=3.D0-2.D0*UIVI DFY(3,IV)=-UI2-GAMMA2 DFY(2,IU)=0.D0 DFY(4,IV)=0.D0 1 CONTINUE DO 2 I=1,N2-2 DFY(1,I+2)=GAMMA DFY(5,I)=GAMMA 2 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON KS PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * CCC include '/home/hairer/fortran/programme/dmebdf/dmebdf.f' IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR LSODE (FULL JACOBIAN) parameter(mmm=9) parameter(nh=2**MMM,n=2*nh,nd=2*nh-2) PARAMETER (LWORK=(31+2)*ND+1,LIWORK=ND) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),u(n),true(nd), + error(nd) REAL*4 TARRAY(2) COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /trans/qq,uzero open(6,file='mebdf.ksres') c ------ FILE DE DONNEES ---------- write(6,2050) 2050 format(1x,'mebdf results for KS') qq=0.025d+0 maxstep=500000 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 C --- SET DEFAULT VALUES MF=23 INDEX=1 LOUT=6 mbnd(1)=0 mbnd(2)=0 mbnd(3)=1 mbnd(4)=1 C --- INITIAL VALUES X=0.0D0 t=0.0D+0 AN=N do 1 i=1,N delx=1.0D+0/an x=delx*(i-1) u1=min(x-0.0d+0,0.1d+0-x) u2=20.0d+0*(x-0.2D+0)*(0.3d+0-x) u3=min(x-0.6d+0,0.7D+0-x) u4=min(x-0.9d0,1.0d0-x) u(i)=16.0D+0*max(0.0d+0,u1,u2,u3,u4) 1 CONTINUE call realft(u,nh,1) do 465 i=1,n u(i)=u(i)/an 465 continue uzero=u(1) do 466 i=1,nd y(i)=u(i+2) 466 continue u(1)=uzero u(2)=0.0D+0 do 467 i=3,n u(i)=y(i-2) 467 continue call realft(u,nh,-1) do 468 i=1,n u(i)=2.0D+0*u(i) 468 continue C --- REQUIRED TOLERANCE RTOL=TOLST atol=RTOL ccc ITOL=2 ITOL=2 itol=1 MAXDER=7 C --- MAXIMAL NUMBER OF STEPS c IWORK(6)=10000 C --- ENDPOINT OF INTEGRATION XEND=100.0D+0 H=1.0D-6 CCC CALL DTIME(TARRAY) it1=mclock() XOUT=100.0D+0 220 CONTINUE C --- CALL OF THE SUBROUTINE CALL DRIVER(ND,t,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,MAXDER,ITOL,RTOL,ATOL) IF (INDEX.EQ.1) THEN INDEX=0 GOTO 220 ELSE 9921 format(1x,f22.16) END IF 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) write(6,*) y(k),true(k) 270 continue do 271 k=11,100,5 error(k)=dabs(y(k)-true(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 271 continue do 272 k=101,1024,50 error(k)=dabs(y(k)-true(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 272 continue sum=dsqrt(sum/dble(nd)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE CCC CALL DTIME(TARRAY) it2=mclock() tarray(1)=(it2-it1)/100.0D+0 WRITE(6,*)TARRAY(1),summ,sumh summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh us=log10(summ) ut=log10(sumh) write(6,*) us,ut ID=0 WRITE(6,*)NFE,NJE,NSTEP,ID,ID,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP 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 F(ND,T,Y,DF) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH) DIMENSION y(ND),df(ND) DIMENSION U(N) COMMON/TRANS/QQ,UZERO c --- copy y to u U(1)=UZERO U(2)=0.D0 DO I=3,N U(I)=Y(I-2) END DO CALL REALFT(U,NH,-1) DO I=1,N U(I)=2.D0*U(I) END DO C --- SQUARE ------ DO I=1,N U(I)=U(I)**2/2.D0 END DO C --- TRANSFORM BACK --- CALL REALFT(U,NH,+1) AN=N DO I=1,N U(I)=U(I)/AN END DO do J=1,NH-1 DIAG=(QQ*J)**2*(1.D0-(QQ*J)**2) DF(2*J-1)=DIAG*Y(2*J-1) +QQ*J*U(2*J+2) DF(2*J )=DIAG*Y(2*J ) -QQ*J*U(2*J+1) enddo end c-------------------------------------------------------------- subroutine PDERV(x,y,dfy,ND,mebnd) PARAMETER (MMM=9) PARAMETER (NH=2**MMM,N=2*NH) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION y(ND),DFY(mebnd,ND) COMMON/TRANS/QQ,UZERO do J=1,NH-1 DIAG=(QQ*J)**2*(1.D0-(QQ*J)**2) DFY(1,2*J-1)=DIAG DFY(1,2*J )=DIAG enddo end PROGRAM DRIMBD implicit double precision(a-h,o-z) DOUBLE PRECISION pi, h, hstart, t, tend, tol DOUBLE PRECISION y(60000),dx,hused,yend(20),ATOL,RTOL DOUBLE PRECISION mi,kcns,dif,err,errm,tout, work(3000000),SUMER double precision error(15),tarray(2) INTEGER iprob,n,mf,liwork,iwork(60000),lwork,lout,mbnd(4) INTEGER NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /ipr/iprob,ntop COMMON /const/kcns,pi,mi,dif,dx common /dim/n common/type/linear logical linear c c Workspace required by the mebdf solver. The arrays c iwork(), work() should have dimension greater than or equal c to liwork, lwork respectively. c c open(1,file='mebres.bigtop5') liwork=60000 lwork=3000000 maxstep=100000 pi=4.0D+0*datan(1.0D+0) write(*,*) write(*,*) ' tol ',' fns ',' bsubs ',' jacs ',' fact', + ' stps ',' time ',' max err',' end err' write(*,*) c ntolmn=2 ntolmx=10 ntoldf=2 nrloop=(ntolmx-ntolmn)*ntoldf+1 tolst=0.1d+0**ntolmn tolfc=0.1d+0**(1.0D+0/ntoldf) do 30 ntol=1,nrloop summ=0.0D+0 sumh=0.0D+0 n=15 mf=21 ml = 1 mu = 1 c c...mu, ml upper, lower bandwidth respectively. c mbnd(1)=ml mbnd(2)=mu mbnd(3)=mu+ml+1 mbnd(4)=2*ml+mu+1 do 51 j=1,15 y(j)=0.0d+0 51 continue t=0.0d0 c tend=1.0d-3 linear=.false. c tout=tend c c... Starting stepsize. c hstart=1.0d-5 h=hstart maxder=7 atol=tolst rtol=atol itol=2 write(6,943) itol 943 format(1x,'the value of itol is',i8) index=1 lout=6 it1 = mclock() 2 continue CALL DRIVER(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,MAXDER,itol,rtol,atol) if(index.eq.1) then index=0 goto 2 endif it2 = mclock() time = .01*float(it2-it1) err=0.0d0 t=tout c c now compute error at end point. c 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 summ=summ*atol sumh=sumh*atol write(6,*) summ,sumh us=log10(summ) ut=log10(sumh) write(6,*) us,ut WRITE(6,*)NFE,NJE,NSTEP,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c c... Right Hand side of y'=f(t,y) c c subroutine f(n,tn,y,dy) integer n double precision tn,y(*),dy(*) c c implicit double precision (a-h,o-z) double precision y0, yN double precision kcns,pi,mi,dif,dx,amult,prod,rat common /const/kcns,pi,mi,dif,dx double precision + c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc, + gamma,delta,sumer, + uin1,uin2,ud1,ud2,ud3,ud4,qud1,qud2,qud3,qud4 parameter (c=1.6d-8,cs=1d-9,cp=1d-8,r=25d3,rp=50d0, + lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4, + rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2, + gamma=40.67286402d-9,delta=17.7493332d0) uin1=0.5d0*sin(2.0d3*pi*tn) uin2=2.0d0*sin(2.0d4*pi*tn) ud1=+y(3)-y(5)-y(7)-uin2 ud2=-y(4)+y(6)-y(7)-uin2 ud3=+y(4)+y(5)+y(7)+uin2 ud4=-y(3)-y(6)+y(7)+uin2 qud1=gamma*(exp(delta*ud1)-1.0d0) qud2=gamma*(exp(delta*ud2)-1.0d0) qud3=gamma*(exp(delta*ud3)-1.0d0) qud4=gamma*(exp(delta*ud4)-1.0d0) dy(1)=(y(8)-0.5d0*y(10)+0.5d0*y(11)+y(14)-y(1)/r)/c dy(2)=(y(9)-0.5d0*y(12)+0.5d0*y(13)+y(15)-y(2)/r)/c dy(3)=(y(10)-qud1+qud4)/cs dy(4)=(-y(11)+qud2-qud3)/cs dy(5)=(y(12)+qud1-qud3)/cs dy(6)=(-y(13)-qud2+qud4)/cs dy(7)=(-y(7)/rp+qud1+qud2-qud3-qud4)/cp dy(8)=-y(1)/lh dy(9)=-y(2)/lh dy(10)=(0.5d0*y(1)-y(3)-rg2*y(10))/ls2 dy(11)=(-0.5d0*y(1)+y(4)-rg3*y(11))/ls3 dy(12)=(0.5d0*y(2)-y(5)-rg2*y(12))/ls2 dy(13)=(-0.5d0*y(2)+y(6)-rg3*y(13))/ls3 dy(14)=(-y(1)+uin1-(ri+rg1)*y(14))/ls1 dy(15)=(-y(2)-(rc+rg1)*y(15))/ls1 return end subroutine pderv(tn,y,pd,n,meband) integer n, meband double precision tn,y(*),pd(meband,n) c c implicit double precision (a-h,o-z) double precision kcns,pi,mi,dif,dx,amult,prod,rat,a1,a2,a3,u common /const/kcns,pi,mi,dif,dx double precision + c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc, + gamma,delta,sumer, + uin2,ud1,ud2,ud3,ud4,qpud1,qpud2,qpud3,qpud4 parameter (c=1.6d-8,cs=1d-9,cp=1d-8,r=25d3,rp=50d0, + lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4, + rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2, + gamma=40.67286402d-9,delta=17.7493332d0) do 177 i=1,n do 177 j=1,n pd(i,j)=0.0D+0 177 continue uin2=2.0d0*sin(2.0d4*pi*tn) ud1=+y(3)-y(5)-y(7)-uin2 ud2=-y(4)+y(6)-y(7)-uin2 ud3=+y(4)+y(5)+y(7)+uin2 ud4=-y(3)-y(6)+y(7)+uin2 qpud1=gamma*delta*exp(delta*ud1) qpud2=gamma*delta*exp(delta*ud2) qpud3=gamma*delta*exp(delta*ud3) qpud4=gamma*delta*exp(delta*ud4) pd(1,1)=-1.0d0/(c*r) pd(1,8)=-1.0d0/c pd(1,10)=-0.5d0/c pd(1,11)=-pd(1,10) pd(1,14)=pd(1,8) pd(2,2)=pd(1,1) pd(2,9)=pd(1,8) pd(2,12)=pd(1,10) pd(2,13)=pd(1,11) pd(2,15)=pd(1,14) pd(3,3)=(-qpud1-qpud4)/cs pd(3,5)=qpud1/cs pd(3,6)=-qpud4/cs pd(3,7)=(qpud1+qpud4)/cs pd(3,10)=1.0d0/cs pd(4,4)=(-qpud2-qpud3)/cs pd(4,5)=-qpud3/cs pd(4,6)=qpud2/cs pd(4,7)=(-qpud2-qpud3)/cs pd(4,11)=-1.0d0/cs pd(5,3)=qpud1/cs pd(5,4)=-qpud3/cs pd(5,5)=(-qpud1-qpud3)/cs pd(5,7)=(-qpud1-qpud3)/cs pd(5,12)=1.0d0/cs pd(6,3)=-qpud4/cs pd(6,4)=qpud2/cs pd(6,6)=(-qpud2-qpud4)/cs pd(6,7)=(qpud2+qpud4)/cs pd(6,13)=-1.0d0/cs pd(7,3)=(qpud1+qpud4)/cp pd(7,4)=(-qpud2-qpud3)/cp pd(7,5)=(-qpud1-qpud3)/cp pd(7,6)=(qpud2+qpud4)/cp pd(7,7)=(-qpud1-qpud2-qpud3-qpud4-1.0d0/rp)/cp pd(8,1)=-1.0d0/lh pd(9,2)=pd(8,1) pd(10,1)=0.5d0/ls2 pd(10,3)=-1.0d0/ls2 pd(10,10)=-rg2/ls2 pd(11,1)=-0.5d0/ls3 pd(11,4)=1.0d0/ls3 pd(11,11)=-rg3/ls3 pd(12,2)=pd(10,1) pd(12,5)=pd(10,3) pd(12,12)=pd(10,10) pd(13,2)=pd(11,1) pd(13,6)=pd(11,4) pd(13,13)=pd(11,11) pd(14,1)=-1.0d0/ls1 pd(14,14)=-(ri+rg1)/ls1 pd(15,2)=pd(14,1) pd(15,15)=-(rc+rg1)/ls1 return end PROGRAM DRIMBD DOUBLE PRECISION pi, h, hstart, t, tend, tol DOUBLE PRECISION dx,hused,yend(20),ATOL,RTOL DOUBLE PRECISION mi,kcns,dif,err,errm,tout,SUMER parameter(n=49999) parameter(lwork=(31+12)*n+1,liwork=n) INTEGER iprob,iwork(liwork),lout,mbnd(4) double precision y(n),work(lwork) INTEGER NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /ipr/iprob,ntop,rtol,atol COMMON /const/kcns,pi,mi,dif,dx common/type/linear logical linear c c Workspace required by the mebdf solver. The arrays c iwork(), work() should have dimension greater than or equal c to liwork, lwork respectively. c open(6,file='mebdfflood.hardres3') maxstep=500000 c c... Diffusion coefficient. c c dif=2.0d-05 dif=4.0d+0/500000.0d+0 c write(6,*) write(6,*) ' tol ',' fns ',' bsubs ',' jacs ',' fact', + ' stps ',' time ',' max err',' end err' write(6,*) ntolmn=2 ntolmx=8 ntoldf=2 nrloop=(ntolmx-ntolmn)*ntoldf+1 tolst=0.1d+0**ntolmn tolfc=0.1d+0**(1.0D+0/ntoldf) do 30 ntol=1,nrloop c... Number of equations. c dx=4.0d0/(n+1.0d0) write(6,*) n,dif c mf = 24 ml = 1 mu = 1 c c...mu, ml upper, lower bandwidth respectively. c mbnd(1)=ml mbnd(2)=mu mbnd(3)=mu+ml+1 mbnd(4)=2*ml+mu+1 c c now set initial conditions. c do 40 j=1,n y(j)=dexp(-(-2.0d0+j*dx)**2) 40 continue c t=0.0d0 tend=4.d0 linear=.true. c tout=tend c c... Starting stepsize. c hstart=1.0d-5 h=hstart maxder=3 atol=tolst rtol=atol itol=2 write(6,469)itol,maxder 469 format(1x,'the value of itol is',i8,'maxder is',i8) index=1 lout=6 it1 = mclock() call driver(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,maxder,itol,rtol,atol) ERRM=0.0D+0 IF(INDEX.LT.0) GOTO 107 2 IF(t.GE.tEND) GOTO 107 IF(t+H.GT.tEND) THEN INDEX=2 TOUT=tEND ELSE INDEX=3 TOUT=tEND ENDIF CALL DRIVER(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,MAXDER,itol,rtol,atol) if(index.lt.0) goto 107 err=0.0d0 c tn=t do 3330 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) + /((atol+rtol*dabs(y(j)))**2) 3330 continue c err=dsqrt(err/n)*atol errm=dmax1(errm,err) goto 2 107 continue it2 = mclock() time = .01*float(it2-it1) err=0.0d0 t=tout c c now compute error at end point. c do 130 j=1,n err=err+((y(j)-(1.d0/dsqrt(1.d0+4.d0*dif*t))* & dexp(-((-2.d0+j*dx-t)**2)/(1.0d0+4.0d0*dif*t)))**2) + /((atol+rtol*dabs(y(j)))**2) 130 continue err=dsqrt(err/n)*atol write(6,*) err,errm us=log10(err) ut=log10(errm) write(6,*) us,ut write(6,600)atol,NFE,NBSOL,NJE,NDEC,NSTEP,time,errm,err 10 continue 600 format (e9.3,I5,I6,I6,I6,I5,F8.2,E9.3,E9.3) 500 format (5(F20.16)) 601 continue tolst=tolst*tolfc 30 continue stop end c c c... Right Hand side of y'=f(t,y) c c subroutine f(n,tn,y,dy) integer n double precision tn,y(*),dy(*) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision y0, yN double precision kcns,pi,mi,dif,dx,amult,prod,rat common /const/kcns,pi,mi,dif,dx 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) c prod=dx*dx amult=2.0d0*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)) dy(1)=dif*(y0-2.d0*y(1)+y(2))/(dx*dx)-(y(2)-y0)/(2.d0*dx) rat=dif/prod do 30 i=2,n-1 dy(i)=rat*(y(i-1)-2.0d0*y(i)+y(i+1))- & (y(i+1)-y(i-1))/(amult) 30 continue dy(n)=dif*(y(n-1)-2.d0*y(n)+yN)/(dx*dx)-(yN-y(n-1))/(2.d0*dx) return end c c c... Jacobian matrix df/dy c c subroutine pderv(tn,y,pd,n,meband) integer n, meband double precision tn,y(*),pd(meband,n) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision kcns,pi,mi,dif,dx,amult,prod,rat,a1,a2,a3,u common /const/kcns,pi,mi,dif,dx 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) prod=dx*dx amult=2.0d0*dx a1=dif/prod a2=-1.0d0/amult a3=-2.0d0*a1 u=a1+a2 do 70 j=2,n pd(1,j)=u 70 continue do 80 j=1,n pd(2,j)=a3 80 continue u=a1-a2 do 90 j=1,n-1 pd(3,j)=u 90 continue return end PROGRAM DRIMBD DOUBLE PRECISION pi, h, hstart, t, tend, tol DOUBLE PRECISION dx,hused,yend(20),ATOL,RTOL DOUBLE PRECISION mi,kcns,dif,err,errm,tout, SUMER parameter(n=1999) PARAMETER (LWORK=(31+12)*N+1,LIWORK=N) INTEGER iprob,iwork(liwork),lout,mbnd(4) double precision y(n),work(lwork) INTEGER NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /ipr/iprob,ntop,atol,rtol COMMON /const/kcns,pi,mi,dif,dx common/type/linear logical linear c c Workspace required by the mebdf solver. The arrays c iwork(), work() should have dimension greater than or equal c to liwork, lwork respectively. c open(6,file='mebdfhard.easyres50') maxstep=500000 ntop=50 c pi=4.0d0*datan(1.0d0) c c... Burgers' equation parameter. c write(6,*) ntop c c... Diffusion coefficient. c dif=2.0d-05 c dif=4.0d+0/500000.0d+0 c c... Stiffness parameter. c kcns=5.0d1 c c... iprob=1 Heat flow, iprob=2 Stiff parabolic, c iprob=3 Advection-Diffusion, iprob=4 Burgers' equation. c c write(*,*) 'Problem number (1-4):' c read(*,*) iprob c write(6,*) write(6,*) ' tol ',' fns ',' bsubs ',' jacs ',' fact', + ' stps ',' time ',' max err',' end err' write(6,*) c c... Solve for tolerances 10**(-i). c do 10 i=2,10 c c... Number of equations. c c n=1999 write(6,8255) n 8255 format(1x,i8) c c... Spatial Mesh (x-mesh) stepsize. c dx=1.0d0/(n+1.0d0) c mf = 24 ml = 1 mu = 1 c c...mu, ml upper, lower bandwidth respectively. c mbnd(1)=ml mbnd(2)=mu mbnd(3)=mu+ml+1 mbnd(4)=2*ml+mu+1 c c now set initial conditions. c 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 t=0.0d0 tend=2.0d+0 linear=.false. c tout=tend c c... Starting stepsize. c hstart=1.0d-5 h=hstart tol = 10.d0**(-i) maxder=7 atol=tol rtol=tol itol=2 index=1 lout=6 it1 = mclock() call driver(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,maxder,itol,rtol,atol) ERRM=0.0D+0 IF(INDEX.LT.0) GOTO 107 2 IF(t.GE.tEND) GOTO 107 IF(t+H.GT.tEND) THEN INDEX=2 TOUT=tEND ELSE INDEX=3 TOUT=tEND ENDIF CALL DRIVER(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,MAXDER,itol,rtol,atol) if(index.lt.0) goto 107 do 120 j=1,n sumer=0.0D+0 do 121 jj=1,ntop sumer=sumer+dexp(-pi*pi*t*jj*jj)*dsin(pi*jj*j*dx) 121 continue err=err+((y(j)-sumer)/(rtol*dabs(y(j))+atol))**2 120 continue err=dsqrt(err/n)*atol if(err.gt.errm) errm=err goto 2 107 continue it2 = mclock() time = .01*float(it2-it1) err=0.0d0 t=tout c c now compute error at end point. c do 1200 j=1,n sumer=0.0D+0 do 1201 jj=1,ntop sumer=sumer+dexp(-pi*pi*t*jj*jj)*dsin(pi*jj*j*dx) 1201 continue err=err+((y(j)-sumer)/(rtol*dabs(y(j))+atol))**2 1200 continue err=dsqrt(err/n)*atol write(6,600)-i,NFE,NBSOL,NJE,NDEC,NSTEP,time,errm,err us=log10(err) ut=log10(errm) write(6,*) us,ut 10 continue 600 format (I5,I5,I6,I6,I6,I5,F8.2,E9.3,E9.3) 500 format (5(F20.16)) stop end c c c... Right Hand side of y'=f(t,y) c c subroutine f(n,tn,y,dy) integer n double precision tn,y(*),dy(*) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision y0, yN double precision kcns,pi,mi,dif,dx,amult,prod,rat common /const/kcns,pi,mi,dif,dx 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) c c... Heat Equation (non insulated rod) c 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 dy(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 dy(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 dy(n)=(y(n-1)-2.0d0*y(n))/(dx*dx) + +y(n)**2-0.5D+0*y(n)*sumer-0.5D+0*sumer**2 return end c c c... Jacobian matrix df/dy c c subroutine pderv(tn,y,pd,n,meband) integer n, meband double precision tn,y(*),pd(meband,n) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision kcns,pi,mi,dif,dx,amult,prod,rat,a1,a2,a3,u common /const/kcns,pi,mi,dif,dx 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) prod=dx*dx amult=2.0d0*dx a1=1.0d0/prod a2=-2.0d0/prod do 141 j=2,n pd(1,j)=a1 141 continue u=a2 do 142 j=1,n sumer=0.0D+0 do 143 k=1,ntop sumer=sumer+dexp(-pi*pi*tn*k*k)*dsin(pi*j*dx*k) 143 continue pd(2,j)=u+2.0D+0*y(j)-0.5D+0*sumer 142 continue do 144 j=1,n-1 pd(3,j)=a1 144 continue return end c C * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR MEBDF ON THE VANDERPOL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * CCC include '/home/hairer/fortran/programme/dmebdf/dmebdf.f' IMPLICIT REAL*8 (A-H,O-Z) C --- PARAMETERS FOR LSODE (FULL JACOBIAN) PARAMETER (ND=6,LWORK=(31+2*ND)*ND+1,LIWORK=ND) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),true(nd) + ,error(nd) REAL*4 TARRAY(2) COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep c ------ FILE DE DONNEES ---------- write(6,2050) 2050 format(1x,'B5 by mebdf') C --- LOOP FOR DIFFERENT TOLERANCES maxstep=100000 NTOLMN=2 NTOLMX=10 NTOLDF=4 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP summ=0.0D+0 sumh=0.0D+0 C --- DIMENSION OF THE SYSTEM N=6 C --- SET DEFAULT VALUES MF=21 INDEX=1 LOUT=6 C --- INITIAL VALUES X=0.0D0 do 784 i=1,6 y(i)=1.0D+0 784 continue XEND=20.0D0 C --- REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL ccc ITOL=2 ITOL=2 write(6,3700) itol 3700 format(1x,'itol is',i8) MAXDER=3 C --- MAXIMAL NUMBER OF STEPS c IWORK(6)=10000 C --- ENDPOINT OF INTEGRATION C XEND=1.0D0 H=1.0D-6 CCC CALL DTIME(TARRAY) it1=mclock() XOUT=20.0D0 220 CONTINUE C --- CALL OF THE SUBROUTINE CALL DRIVER(N,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,MAXDER,ITOL,RTOL,ATOL) if(x.gt.xout) goto 955 if(index.eq.0) goto 955 INDEX=3 c this is zero for end value and 3 for max value 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 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 us=log10(summ) ut=log10(sumh) write(6,*) us,ut WRITE(6,*)NFE,NJE,NSTEP,NDEC,NBSOL WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c SUBROUTINE F(N,X,Y,DF) C --- RIGHT-HAND SIDE OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DF(N) df(1)=-10.0D+0*y(1)+1000.0D+0*y(2) df(2)=-1000.0D+0*y(1)-10.0D+0*y(2) df(3)=-4.0D+0*y(3) df(4)=-y(4) df(5)=-0.5D+0*y(5) df(6)=-0.1D+0*y(6) RETURN END C SUBROUTINE PDERV(X,Y,DFY,N,MEBAND) C --- JACOBIAN OF VANDERPOL EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(MEBAND,N) DFY(1,1)=-10.0D0 DFY(1,2)=1000.D0 DFY(2,1)=-1000.0D+0 DFY(2,2)=-10.0D+0 dfy(3,3)=-4.0D+0 dfy(4,4)=-1.0D+0 dfy(5,5)=-0.5D+0 dfy(6,6)=-0.1D+0 RETURN END PROGRAM DRIMBD implicit double precision(a-h,o-z) DOUBLE PRECISION pi, h, hstart, t, tend, tol DOUBLE PRECISION y(60000),dx,hused,yend(20),ATOL,RTOL DOUBLE PRECISION mi,kcns,dif,err,errm,tout, work(3000000),SUMER double precision error(15),tarray(2) INTEGER iprob,n,mf,liwork,iwork(60000),lwork,lout,mbnd(4) INTEGER NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /ipr/iprob,ntop COMMON /const/kcns,pi,mi,dif,dx common /dim/n common/type/linear logical linear c c Workspace required by the mebdf solver. The arrays c iwork(), work() should have dimension greater than or equal c to liwork, lwork respectively. c open(6,file='mebdfhardring.newres4') liwork=60000 lwork=3000000 maxstep=5000000 pi=4.0D+0*datan(1.0D+0) write(6,*) write(6,*) ' tol ',' fns ',' bsubs ',' jacs ',' fact', + ' stps ',' time ',' max err',' end err' write(6,*) c ntolmn=2 ntolmx=10 ntoldf=2 nrloop=(ntolmx-ntolmn)*ntoldf+1 tolst=0.1d+0**ntolmn tolfc=0.1d+0**(1.0D+0/ntoldf) do 30 ntol=1,nrloop summ=0.0D+0 sumh=0.0D+0 n=15 mf=21 ml = 1 mu = 1 c c...mu, ml upper, lower bandwidth respectively. c mbnd(1)=ml mbnd(2)=mu mbnd(3)=mu+ml+1 mbnd(4)=2*ml+mu+1 do 51 j=1,15 y(j)=0.0d+0 51 continue t=0.0d0 c tend=1.0d-3 linear=.false. c tout=tend c c... Starting stepsize. c hstart=1.0d-5 h=hstart maxder=7 atol=tolst rtol=atol itol=2 write(6,943) itol 943 format(1x,'the value of itol is',i8) index=1 lout=6 it1 = mclock() 2 continue CALL DRIVER(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,MAXDER,itol,rtol,atol) if(index.eq.1) then index=0 goto 2 endif it2 = mclock() time = .01*float(it2-it1) err=0.0d0 t=tout c c now compute error at end point. c 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 write(6,*) summ,sumh us=log10(summ) ut=log10(sumh) write(6,*) us,ut WRITE(6,*)NFE,NJE,NSTEP,NDEC,NBSOL,npset,ncoset WRITE(6,*)' ***** TOL=',RTOL,' ELAPSED TIME=',TARRAY(1),' ****' C -------- NEW TOLERANCE --- 25 TOLST=TOLST*TOLFC 30 CONTINUE STOP END c c... Right Hand side of y'=f(t,y) c c subroutine f(n,tn,y,dy) integer n double precision tn,y(*),dy(*) c c implicit double precision (a-h,o-z) double precision y0, yN double precision kcns,pi,mi,dif,dx,amult,prod,rat common /const/kcns,pi,mi,dif,dx double precision + c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc, + gamma,delta,sumer, + uin1,uin2,ud1,ud2,ud3,ud4,qud1,qud2,qud3,qud4 parameter (c=1.6d-8,cs=2.0d-12,cp=1d-8,r=25d3,rp=50d0, + lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4, + rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2, + gamma=40.67286402d-9,delta=17.7493332d0) uin1=0.5d0*sin(2.0d3*pi*tn) uin2=2.0d0*sin(2.0d4*pi*tn) ud1=+y(3)-y(5)-y(7)-uin2 ud2=-y(4)+y(6)-y(7)-uin2 ud3=+y(4)+y(5)+y(7)+uin2 ud4=-y(3)-y(6)+y(7)+uin2 qud1=gamma*(exp(delta*ud1)-1.0d0) qud2=gamma*(exp(delta*ud2)-1.0d0) qud3=gamma*(exp(delta*ud3)-1.0d0) qud4=gamma*(exp(delta*ud4)-1.0d0) dy(1)=(y(8)-0.5d0*y(10)+0.5d0*y(11)+y(14)-y(1)/r)/c dy(2)=(y(9)-0.5d0*y(12)+0.5d0*y(13)+y(15)-y(2)/r)/c dy(3)=(y(10)-qud1+qud4)/cs dy(4)=(-y(11)+qud2-qud3)/cs dy(5)=(y(12)+qud1-qud3)/cs dy(6)=(-y(13)-qud2+qud4)/cs dy(7)=(-y(7)/rp+qud1+qud2-qud3-qud4)/cp dy(8)=-y(1)/lh dy(9)=-y(2)/lh dy(10)=(0.5d0*y(1)-y(3)-rg2*y(10))/ls2 dy(11)=(-0.5d0*y(1)+y(4)-rg3*y(11))/ls3 dy(12)=(0.5d0*y(2)-y(5)-rg2*y(12))/ls2 dy(13)=(-0.5d0*y(2)+y(6)-rg3*y(13))/ls3 dy(14)=(-y(1)+uin1-(ri+rg1)*y(14))/ls1 dy(15)=(-y(2)-(rc+rg1)*y(15))/ls1 return end subroutine pderv(tn,y,pd,n,meband) integer n, meband double precision tn,y(*),pd(meband,n) c c implicit double precision (a-h,o-z) double precision kcns,pi,mi,dif,dx,amult,prod,rat,a1,a2,a3,u common /const/kcns,pi,mi,dif,dx double precision + c,cs,cp,r,rp,lh,ls1,ls2,ls3,rg1,rg2,rg3,ri,rc, + gamma,delta,sumer, + uin2,ud1,ud2,ud3,ud4,qpud1,qpud2,qpud3,qpud4 parameter (c=1.6d-8,cs=2.0d-12,cp=1d-8,r=25d3,rp=50d0, + lh=4.45d0,ls1=2d-3,ls2=5d-4,ls3=5d-4, + rg1=36.3d0,rg2=17.3d0,rg3=17.3d0,ri=5d1,rc=6d2, + gamma=40.67286402d-9,delta=17.7493332d0) do 177 i=1,n do 177 j=1,n pd(i,j)=0.0D+0 177 continue uin2=2.0d0*sin(2.0d4*pi*tn) ud1=+y(3)-y(5)-y(7)-uin2 ud2=-y(4)+y(6)-y(7)-uin2 ud3=+y(4)+y(5)+y(7)+uin2 ud4=-y(3)-y(6)+y(7)+uin2 qpud1=gamma*delta*exp(delta*ud1) qpud2=gamma*delta*exp(delta*ud2) qpud3=gamma*delta*exp(delta*ud3) qpud4=gamma*delta*exp(delta*ud4) pd(1,1)=-1.0d0/(c*r) pd(1,8)=-1.0d0/c pd(1,10)=-0.5d0/c pd(1,11)=-pd(1,10) pd(1,14)=pd(1,8) pd(2,2)=pd(1,1) pd(2,9)=pd(1,8) pd(2,12)=pd(1,10) pd(2,13)=pd(1,11) pd(2,15)=pd(1,14) pd(3,3)=(-qpud1-qpud4)/cs pd(3,5)=qpud1/cs pd(3,6)=-qpud4/cs pd(3,7)=(qpud1+qpud4)/cs pd(3,10)=1.0d0/cs pd(4,4)=(-qpud2-qpud3)/cs pd(4,5)=-qpud3/cs pd(4,6)=qpud2/cs pd(4,7)=(-qpud2-qpud3)/cs pd(4,11)=-1.0d0/cs pd(5,3)=qpud1/cs pd(5,4)=-qpud3/cs pd(5,5)=(-qpud1-qpud3)/cs pd(5,7)=(-qpud1-qpud3)/cs pd(5,12)=1.0d0/cs pd(6,3)=-qpud4/cs pd(6,4)=qpud2/cs pd(6,6)=(-qpud2-qpud4)/cs pd(6,7)=(qpud2+qpud4)/cs pd(6,13)=-1.0d0/cs pd(7,3)=(qpud1+qpud4)/cp pd(7,4)=(-qpud2-qpud3)/cp pd(7,5)=(-qpud1-qpud3)/cp pd(7,6)=(qpud2+qpud4)/cp pd(7,7)=(-qpud1-qpud2-qpud3-qpud4-1.0d0/rp)/cp pd(8,1)=-1.0d0/lh pd(9,2)=pd(8,1) pd(10,1)=0.5d0/ls2 pd(10,3)=-1.0d0/ls2 pd(10,10)=-rg2/ls2 pd(11,1)=-0.5d0/ls3 pd(11,4)=1.0d0/ls3 pd(11,11)=-rg3/ls3 pd(12,2)=pd(10,1) pd(12,5)=pd(10,3) pd(12,12)=pd(10,10) pd(13,2)=pd(11,1) pd(13,6)=pd(11,4) pd(13,13)=pd(11,11) pd(14,1)=-1.0d0/ls1 pd(14,14)=-(ri+rg1)/ls1 pd(15,2)=pd(14,1) pd(15,15)=-(rc+rg1)/ls1 return end PROGRAM DRIMBD DOUBLE PRECISION pi, h, hstart, t, tend, tol DOUBLE PRECISION dx,hused,yend(20),ATOL,RTOL DOUBLE PRECISION mi,kcns,dif,err,errm,tout, SUMER parameter(n=1999) PARAMETER (LWORK=(31+12)*N+1,LIWORK=N) INTEGER iprob,iwork(liwork),lout,mbnd(4) double precision y(n),work(lwork) INTEGER NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /ipr/iprob,ntop,atol,rtol COMMON /const/kcns,pi,mi,dif,dx common/type/linear logical linear c c Workspace required by the mebdf solver. The arrays c iwork(), work() should have dimension greater than or equal c to liwork, lwork respectively. c open(6,file='mebdfheat.easynores50') maxstep=500000 ntop=50 c pi=4.0d0*datan(1.0d0) c c... Burgers' equation parameter. c write(6,*) ntop c c... Diffusion coefficient. c dif=2.0d-05 c dif=4.0d+0/500000.0d+0 c c... Stiffness parameter. c kcns=5.0d1 c c... iprob=1 Heat flow, iprob=2 Stiff parabolic, c iprob=3 Advection-Diffusion, iprob=4 Burgers' equation. c c write(*,*) 'Problem number (1-4):' c read(*,*) iprob c write(6,*) write(6,*) ' tol ',' fns ',' bsubs ',' jacs ',' fact', + ' stps ',' time ',' max err',' end err' write(6,*) c c... Solve for tolerances 10**(-i). c do 10 i=2,10 c c... Number of equations. c c n=1999 write(6,8255) n 8255 format(1x,i8) c c... Spatial Mesh (x-mesh) stepsize. c dx=1.0d0/(n+1.0d0) c mf = 24 ml = 1 mu = 1 c c...mu, ml upper, lower bandwidth respectively. c mbnd(1)=ml mbnd(2)=mu mbnd(3)=mu+ml+1 mbnd(4)=2*ml+mu+1 c c now set initial conditions. c 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 t=0.0d0 tend=2.0d+0 linear=.false. c tout=tend c c... Starting stepsize. c hstart=1.0d-5 h=hstart tol = 10.d0**(-i) maxder=7 atol=tol rtol=tol itol=2 index=1 lout=6 it1 = mclock() call driver(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,maxder,itol,rtol,atol) ERRM=0.0D+0 IF(INDEX.LT.0) GOTO 107 2 IF(t.GE.tEND) GOTO 107 IF(t+H.GT.tEND) THEN INDEX=2 TOUT=tEND ELSE INDEX=0 TOUT=tEND ENDIF CALL DRIVER(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,MAXDER,itol,rtol,atol) if(index.lt.0) goto 107 do 120 j=1,n sumer=0.0D+0 do 121 jj=1,ntop sumer=sumer+dexp(-pi*pi*t*jj*jj)*dsin(pi*jj*j*dx) 121 continue err=err+((y(j)-sumer)/(rtol*dabs(y(j))+atol))**2 120 continue err=dsqrt(err/n)*atol if(err.gt.errm) errm=err goto 2 107 continue it2 = mclock() time = .01*float(it2-it1) err=0.0d0 t=tout c c now compute error at end point. c do 1200 j=1,n sumer=0.0D+0 do 1201 jj=1,ntop sumer=sumer+dexp(-pi*pi*t*jj*jj)*dsin(pi*jj*j*dx) 1201 continue err=err+((y(j)-sumer)/(rtol*dabs(y(j))+atol))**2 1200 continue err=dsqrt(err/n)*atol write(6,600)-i,NFE,NBSOL,NJE,NDEC,NSTEP,time,errm,err us=log10(err) ut=log10(errm) write(6,*) us,ut 10 continue 600 format (I5,I5,I6,I6,I6,I5,F8.2,E9.3,E9.3) 500 format (5(F20.16)) stop end c c c... Right Hand side of y'=f(t,y) c c subroutine f(n,tn,y,dy) integer n double precision tn,y(*),dy(*) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision y0, yN double precision kcns,pi,mi,dif,dx,amult,prod,rat common /const/kcns,pi,mi,dif,dx 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) c c... Heat Equation (non insulated rod) c 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 dy(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 dy(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 dy(n)=(y(n-1)-2.0d0*y(n))/(dx*dx) + +y(n)**2-0.5D+0*y(n)*sumer-0.5D+0*sumer**2 return end c c c... Jacobian matrix df/dy c c subroutine pderv(tn,y,pd,n,meband) integer n, meband double precision tn,y(*),pd(meband,n) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision kcns,pi,mi,dif,dx,amult,prod,rat,a1,a2,a3,u common /const/kcns,pi,mi,dif,dx 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) prod=dx*dx amult=2.0d0*dx a1=1.0d0/prod a2=-2.0d0/prod do 141 j=2,n pd(1,j)=a1 141 continue u=a2 do 142 j=1,n sumer=0.0D+0 do 143 k=1,ntop sumer=sumer+dexp(-pi*pi*tn*k*k)*dsin(pi*j*dx*k) 143 continue pd(2,j)=u+2.0D+0*y(j)-0.5D+0*sumer 142 continue do 144 j=1,n-1 pd(3,j)=a1 144 continue return end c PROGRAM DRIMBD DOUBLE PRECISION pi, h, hstart, t, tend, tol DOUBLE PRECISION dx,hused,yend(20),ATOL,RTOL DOUBLE PRECISION mi,kcns,dif,err,errm,tout, SUMER parameter(n=49999) PARAMETER (LWORK=(31+12)*N+1,LIWORK=N) INTEGER iprob,iwork(liwork),lout,mbnd(4) double precision y(n),work(lwork) INTEGER NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /ipr/iprob,ntop,atol,rtol COMMON /const/kcns,pi,mi,dif,dx common/type/linear logical linear c c Workspace required by the mebdf solver. The arrays c iwork(), work() should have dimension greater than or equal c to liwork, lwork respectively. c open(6,file='mebdfheat.hardnores50') maxstep=500000 ntop=50 c pi=4.0d0*datan(1.0d0) c c... Burgers' equation parameter. c write(6,*) ntop c c... Diffusion coefficient. c c dif=2.0d-05 dif=4.0d+0/500000.0d+0 c c... Stiffness parameter. c kcns=5.0d1 c c... iprob=1 Heat flow, iprob=2 Stiff parabolic, c iprob=3 Advection-Diffusion, iprob=4 Burgers' equation. c c write(*,*) 'Problem number (1-4):' c read(*,*) iprob c write(6,*) write(6,*) ' tol ',' fns ',' bsubs ',' jacs ',' fact', + ' stps ',' time ',' max err',' end err' write(6,*) c c... Solve for tolerances 10**(-i). c do 10 i=2,10 c c... Number of equations. c c n=49999 write(6,8255) n 8255 format(1x,i8) c c... Spatial Mesh (x-mesh) stepsize. c dx=1.0d0/(n+1.0d0) c mf = 24 ml = 1 mu = 1 c c...mu, ml upper, lower bandwidth respectively. c mbnd(1)=ml mbnd(2)=mu mbnd(3)=mu+ml+1 mbnd(4)=2*ml+mu+1 c c now set initial conditions. c 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 t=0.0d0 tend=2.0d+0 linear=.false. c tout=tend c c... Starting stepsize. c hstart=1.0d-5 h=hstart tol = 10.d0**(-i) maxder=7 atol=tol rtol=tol itol=2 index=1 lout=6 it1 = mclock() call driver(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,maxder,itol,rtol,atol) ERRM=0.0D+0 IF(INDEX.LT.0) GOTO 107 2 IF(t.GE.tEND) GOTO 107 IF(t+H.GT.tEND) THEN INDEX=2 TOUT=tEND ELSE INDEX=0 TOUT=tEND ENDIF CALL DRIVER(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,MAXDER,itol,rtol,atol) if(index.lt.0) goto 107 do 120 j=1,n sumer=0.0D+0 do 121 jj=1,ntop sumer=sumer+dexp(-pi*pi*t*jj*jj)*dsin(pi*jj*j*dx) 121 continue err=err+((y(j)-sumer)/(rtol*dabs(y(j))+atol))**2 120 continue err=dsqrt(err/n)*atol if(err.gt.errm) errm=err goto 2 107 continue it2 = mclock() time = .01*float(it2-it1) err=0.0d0 t=tout c c now compute error at end point. c do 1200 j=1,n sumer=0.0D+0 do 1201 jj=1,ntop sumer=sumer+dexp(-pi*pi*t*jj*jj)*dsin(pi*jj*j*dx) 1201 continue err=err+((y(j)-sumer)/(rtol*dabs(y(j))+atol))**2 1200 continue err=dsqrt(err/n)*atol write(6,600)-i,NFE,NBSOL,NJE,NDEC,NSTEP,time,errm,err us=log10(err) ut=log10(errm) write(6,*) us,ut 10 continue 600 format (I5,I5,I6,I6,I6,I5,F8.2,E9.3,E9.3) 500 format (5(F20.16)) stop end c c c... Right Hand side of y'=f(t,y) c c subroutine f(n,tn,y,dy) integer n double precision tn,y(*),dy(*) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision y0, yN double precision kcns,pi,mi,dif,dx,amult,prod,rat common /const/kcns,pi,mi,dif,dx 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) c c... Heat Equation (non insulated rod) c 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 dy(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 dy(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 dy(n)=(y(n-1)-2.0d0*y(n))/(dx*dx) + +y(n)**2-0.5D+0*y(n)*sumer-0.5D+0*sumer**2 return end c c c... Jacobian matrix df/dy c c subroutine pderv(tn,y,pd,n,meband) integer n, meband double precision tn,y(*),pd(meband,n) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision kcns,pi,mi,dif,dx,amult,prod,rat,a1,a2,a3,u common /const/kcns,pi,mi,dif,dx 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) prod=dx*dx amult=2.0d0*dx a1=1.0d0/prod a2=-2.0d0/prod do 141 j=2,n pd(1,j)=a1 141 continue u=a2 do 142 j=1,n sumer=0.0D+0 do 143 k=1,ntop sumer=sumer+dexp(-pi*pi*tn*k*k)*dsin(pi*j*dx*k) 143 continue pd(2,j)=u+2.0D+0*y(j)-0.5D+0*sumer 142 continue do 144 j=1,n-1 pd(3,j)=a1 144 continue return end c PROGRAM DRIMBD DOUBLE PRECISION pi, h, hstart, t, tend, tol DOUBLE PRECISION dx,hused,yend(20),ATOL,RTOL DOUBLE PRECISION mi,kcns,dif,err,errm,tout,SUMER parameter(n=1999) parameter(lwork=(31+12)*n+1,liwork=n) INTEGER iprob,iwork(liwork),lout,mbnd(4) double precision y(n),work(lwork) INTEGER NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /ipr/iprob,ntop,rtol,atol COMMON /const/kcns,pi,mi,dif,dx common/type/linear logical linear c c Workspace required by the mebdf solver. The arrays c iwork(), work() should have dimension greater than or equal c to liwork, lwork respectively. c open(6,file='mebdfflood.easynores7') maxstep=500000 c c... Diffusion coefficient. c dif=2.0d-05 c dif=4.0d+0/500000.0d+0 c write(6,*) write(6,*) ' tol ',' fns ',' bsubs ',' jacs ',' fact', + ' stps ',' time ',' max err',' end err' write(6,*) ntolmn=2 ntolmx=9 ntoldf=2 nrloop=(ntolmx-ntolmn)*ntoldf+1 tolst=0.1d+0**ntolmn tolfc=0.1d+0**(1.0D+0/ntoldf) do 30 ntol=1,nrloop c... Number of equations. c dx=4.0d0/(n+1.0d0) write(6,*) n,dif c mf = 24 ml = 1 mu = 1 c c...mu, ml upper, lower bandwidth respectively. c mbnd(1)=ml mbnd(2)=mu mbnd(3)=mu+ml+1 mbnd(4)=2*ml+mu+1 c c now set initial conditions. c do 40 j=1,n y(j)=dexp(-(-2.0d0+j*dx)**2) 40 continue c t=0.0d0 tend=4.d0 linear=.true. c tout=tend c c... Starting stepsize. c hstart=1.0d-5 h=hstart maxder=7 atol=tolst rtol=atol itol=2 write(6,469)itol,maxder 469 format(1x,'the value of itol is',i8,'maxder is',i8) index=1 lout=6 it1 = mclock() call driver(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,maxder,itol,rtol,atol) ERRM=0.0D+0 IF(INDEX.LT.0) GOTO 107 2 IF(t.GE.tEND) GOTO 107 IF(t+H.GT.tEND) THEN INDEX=2 TOUT=tEND ELSE INDEX=3 TOUT=tEND ENDIF CALL DRIVER(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,MAXDER,itol,rtol,atol) if(index.lt.0) goto 107 err=0.0d0 c tn=t do 3330 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) + /((atol+rtol*dabs(y(j)))**2) 3330 continue c err=dsqrt(err/n)*atol errm=dmax1(errm,err) goto 2 107 continue it2 = mclock() time = .01*float(it2-it1) err=0.0d0 t=tout c c now compute error at end point. c do 130 j=1,n err=err+((y(j)-(1.d0/dsqrt(1.d0+4.d0*dif*t))* & dexp(-((-2.d0+j*dx-t)**2)/(1.0d0+4.0d0*dif*t)))**2) + /((atol+rtol*dabs(y(j)))**2) 130 continue err=dsqrt(err/n)*atol write(6,*) err,errm us=log10(err) ut=log10(errm) write(6,*) us,ut write(6,600)atol,NFE,NBSOL,NJE,NDEC,NSTEP,time,errm,err 10 continue 600 format (e9.3,I5,I6,I6,I6,I5,F8.2,E9.3,E9.3) 500 format (5(F20.16)) 601 continue tolst=tolst*tolfc 30 continue stop end c c c... Right Hand side of y'=f(t,y) c c subroutine f(n,tn,y,dy) integer n double precision tn,y(*),dy(*) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision y0, yN double precision kcns,pi,mi,dif,dx,amult,prod,rat common /const/kcns,pi,mi,dif,dx 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) c prod=dx*dx amult=2.0d0*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)) dy(1)=dif*(y0-2.d0*y(1)+y(2))/(dx*dx)-(y(2)-y0)/(2.d0*dx) rat=dif/prod do 30 i=2,n-1 dy(i)=rat*(y(i-1)-2.0d0*y(i)+y(i+1))- & (y(i+1)-y(i-1))/(amult) 30 continue dy(n)=dif*(y(n-1)-2.d0*y(n)+yN)/(dx*dx)-(yN-y(n-1))/(2.d0*dx) return end c c c... Jacobian matrix df/dy c c subroutine pderv(tn,y,pd,n,meband) integer n, meband double precision tn,y(*),pd(meband,n) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision kcns,pi,mi,dif,dx,amult,prod,rat,a1,a2,a3,u common /const/kcns,pi,mi,dif,dx 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) prod=dx*dx amult=2.0d0*dx a1=dif/prod a2=-1.0d0/amult a3=-2.0d0*a1 u=a1+a2 do 70 j=2,n pd(1,j)=u 70 continue do 80 j=1,n pd(2,j)=a3 80 continue u=a1-a2 do 90 j=1,n-1 pd(3,j)=u 90 continue return end PROGRAM DRIMBD DOUBLE PRECISION pi, h, hstart, t, tend, tol DOUBLE PRECISION dx,hused,yend(20),ATOL,RTOL DOUBLE PRECISION mi,kcns,dif,err,errm,tout,SUMER parameter(n=49999) parameter(lwork=(31+12)*n+1,liwork=n) INTEGER iprob,iwork(liwork),lout,mbnd(4) double precision y(n),work(lwork) INTEGER NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD COMMON /MEBDF2/HUSED,NQUSED,NSTEP,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET, + MAXORD,maxstep COMMON /ipr/iprob,ntop,rtol,atol COMMON /const/kcns,pi,mi,dif,dx common/type/linear logical linear c c Workspace required by the mebdf solver. The arrays c iwork(), work() should have dimension greater than or equal c to liwork, lwork respectively. c open(6,file='mebdfflood.nohardres7') maxstep=500000 c c... Diffusion coefficient. c c dif=2.0d-05 dif=4.0d+0/500000.0d+0 c write(6,*) write(6,*) ' tol ',' fns ',' bsubs ',' jacs ',' fact', + ' stps ',' time ',' max err',' end err' write(6,*) ntolmn=2 ntolmx=9 ntoldf=2 nrloop=(ntolmx-ntolmn)*ntoldf+1 tolst=0.1d+0**ntolmn tolfc=0.1d+0**(1.0D+0/ntoldf) do 30 ntol=1,nrloop c... Number of equations. c dx=4.0d0/(n+1.0d0) write(6,*) n,dif c mf = 24 ml = 1 mu = 1 c c...mu, ml upper, lower bandwidth respectively. c mbnd(1)=ml mbnd(2)=mu mbnd(3)=mu+ml+1 mbnd(4)=2*ml+mu+1 c c now set initial conditions. c do 40 j=1,n y(j)=dexp(-(-2.0d0+j*dx)**2) 40 continue c t=0.0d0 tend=4.d0 linear=.true. c tout=tend c c... Starting stepsize. c hstart=1.0d-5 h=hstart maxder=7 atol=tolst rtol=atol itol=2 write(6,469)itol,maxder 469 format(1x,'the value of itol is',i8,'maxder is',i8) index=1 lout=6 it1 = mclock() call driver(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,maxder,itol,rtol,atol) ERRM=0.0D+0 IF(INDEX.LT.0) GOTO 107 2 IF(t.GE.tEND) GOTO 107 IF(t+H.GT.tEND) THEN INDEX=2 TOUT=tEND ELSE INDEX=0 TOUT=tEND ENDIF CALL DRIVER(n,t,h,y,tout,tend,mf,index,lout,lwork, + work,liwork,iwork,mbnd,MAXDER,itol,rtol,atol) if(index.lt.0) goto 107 err=0.0d0 c tn=t do 3330 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) + /((atol+rtol*dabs(y(j)))**2) 3330 continue c err=dsqrt(err/n)*atol errm=dmax1(errm,err) goto 2 107 continue it2 = mclock() time = .01*float(it2-it1) err=0.0d0 t=tout c c now compute error at end point. c do 130 j=1,n err=err+((y(j)-(1.d0/dsqrt(1.d0+4.d0*dif*t))* & dexp(-((-2.d0+j*dx-t)**2)/(1.0d0+4.0d0*dif*t)))**2) + /((atol+rtol*dabs(y(j)))**2) 130 continue err=dsqrt(err/n)*atol write(6,*) err,errm us=log10(err) ut=log10(errm) write(6,*) us,ut write(6,600)atol,NFE,NBSOL,NJE,NDEC,NSTEP,time,errm,err 10 continue 600 format (e9.3,I5,I6,I6,I6,I5,F8.2,E9.3,E9.3) 500 format (5(F20.16)) 601 continue tolst=tolst*tolfc 30 continue stop end c c c... Right Hand side of y'=f(t,y) c c subroutine f(n,tn,y,dy) integer n double precision tn,y(*),dy(*) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision y0, yN double precision kcns,pi,mi,dif,dx,amult,prod,rat common /const/kcns,pi,mi,dif,dx 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) c prod=dx*dx amult=2.0d0*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)) dy(1)=dif*(y0-2.d0*y(1)+y(2))/(dx*dx)-(y(2)-y0)/(2.d0*dx) rat=dif/prod do 30 i=2,n-1 dy(i)=rat*(y(i-1)-2.0d0*y(i)+y(i+1))- & (y(i+1)-y(i-1))/(amult) 30 continue dy(n)=dif*(y(n-1)-2.d0*y(n)+yN)/(dx*dx)-(yN-y(n-1))/(2.d0*dx) return end c c c... Jacobian matrix df/dy c c subroutine pderv(tn,y,pd,n,meband) integer n, meband double precision tn,y(*),pd(meband,n) c c implicit double precision (a-h,o-z) integer iprob common /ipr/iprob,ntop,atol,rtol double precision kcns,pi,mi,dif,dx,amult,prod,rat,a1,a2,a3,u common /const/kcns,pi,mi,dif,dx 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) prod=dx*dx amult=2.0d0*dx a1=dif/prod a2=-1.0d0/amult a3=-2.0d0*a1 u=a1+a2 do 70 j=2,n pd(1,j)=u 70 continue do 80 j=1,n pd(2,j)=a3 80 continue u=a1-a2 do 90 j=1,n-1 pd(3,j)=u 90 continue return end