C************************************************************************* C 02/04/2000 C************************************************************************* C DRIVER FOR MEBDFDAE ON THE E5 PROBLEM C************************************************************************* IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (ND=4,LWORK=(34+3*ND)*ND+3,LIWORK=ND+14) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),TRUE(ND), + ERROR(ND),MASBND(4) CHARACTER*(*) SOLVER, PROBLEM PARAMETER (SOLVER = ' MEBDFDAE',PROBLEM = 'E5') C EXTERNAL F,PDERV,MAS C open(6,file='e5.res') rewind(6) C write(6,2050)PROBLEM,SOLVER 2050 FORMAT(1X,'COMPUTATIONAL STATISTICS OF THE ',A, & ' PROBLEM USING ', A/) WRITE(6,2051)ND 2051 FORMAT(1X,'NUMBER OF EQUATIONS :',I2) WRITE(6,*)' ' C... LOOP FOR DIFFERENT TOLERANCES NTOLMN=2 NTOLMX=10 NTOLDF=2 NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1 TOLST=0.1D0**NTOLMN TOLFC=0.1D0**(1.D0/NTOLDF) DO 30 NTOL=1,NRLOOP SUMM=0.0D+0 SUMH=0.0D+0 C... DIMENSION OF THE SYSTEM N=4 C... ENDPOINT OF INTEGRATION X=0.0D0 XEND=10000000000000.0D+0 C... REQUIRED TOLERANCE RTOL=TOLST ATOL=1.7D-24 ITOL=2 WRITE(6,*) 'RESULT WITH THE FOLLOWING TOL :' WRITE(6,*) 'RTOL =' ,RTOL WRITE(6,*) 'ATOL =' ,ATOL WRITE(6,*) C... INITIAL VALUES CALL INIT(N,X,Y) C... SET DEFAULT VALUES MF=21 INDEX=1 LOUT=6 MAXDER=6 C... MAXIMAL NUMBER OF STEPS IWORK(14)=100000 WORK(1)=0.0D+0 H=1.0D-6 XOUT=10.0D+0 MASBND(1)=0 C it1=mclock() C... C... CALL OF THE SUBROUTINE C... DO 20 I=1,7 220 CONTINUE CALL MEBDF(N,X,H,Y,XOUT,XEND,MF,INDEX,LOUT,LWORK, + WORK,LIWORK,IWORK,MBND,MASBND,MAXDER,ITOL,RTOL,ATOL, + rpar,ipar,f,pderv,mas,ierr) IF (INDEX.EQ.1) THEN INDEX=0 GOTO 220 ELSE XOUT=XOUT*100.0D+0 END IF C... C... CALCULATE AND PRINT THE ERROR AT THE END POINT C... CALL SOLN(N,X,TRUE,I) sum=0.0D+0 do 270 k=1,n error(k)=dabs(true(k)-y(k)) sum=sum+((error(k)/(rtol*dabs(y(k))+atol))**2) 270 continue sum=dsqrt(sum/dble(n)) sumh=sumh+sum summ=max(summ,sum) 20 CONTINUE it2=mclock() TIME=(it2-it1)/100.0D+0 summ=summ*rtol sumh=sumh*rtol/7.0D+0 us=-log10(summ) ut=-log10(sumh) C... C... STATISTICS OF THE SIMULATION C... HUSED = WORK(2) NQUSED = IWORK(4) NSTEP = IWORK(5) NFAIL = IWORK(6) NFE = IWORK(7) NJE = IWORK(8) NDEC = IWORK(9) NBSOL = IWORK(10) MAXORD = IWORK(13) C... C... THE CURRENT RUN IS COMPLETE, SO PRINT THE COMPUTATIONAL STAT- C... ISTICS FOR MEBDF AND GO ON TO THE NEXT RUN WRITE(LOUT,8)HUSED,NQUSED,MAXORD,NSTEP,NFAIL,NFE,NJE,NDEC, + NBSOL,US,TIME 8 FORMAT(1H / 1 ' LAST STEP SIZE ', D13.6,/, 2 ' LAST ORDER OF THE METHOD ', I10,/, 3 ' MAXIMUM ORDER USED SO FAR ', I10,/, 4 ' TOTAL NUMBER OF STEPS TAKEN ', I10,/, 5 ' TOTAL NUMBER OF FAILED STEPS ', I10,/, 6 ' NUMBER OF FUNCTION EVALUATIONS ', I10,/, 7 ' NUMBER OF JACOBIAN EVALUATIONS ', I10,/, 8 ' NUMBER OF FACTORIZATION ', I10,/, 9 ' NUMBER OF BACKSOLVES ', I10,/, + ' NUMBER OF CORRECT DIGITS ', F10.3,/, + ' CPU TIME ', F10.4) WRITE(LOUT,*) '---------------------------------------------' C... C... NEW TOLERANCE C... 25 TOLST=TOLST*TOLFC 30 CONTINUE close(6) STOP END C------------------------------------------------------------------------ SUBROUTINE F(N,X,Y,DF,IPAR,RPAR,IER) C... RIGHT-HAND SIDE OF E5 EQUATION IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),DF(N),IPAR(*),RPAR(*) C 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) C RETURN END C----------------------------------------------------------------------- SUBROUTINE PDERV(X,Y,DFY,N,MEBND,ipar,rpar,ierr) C --- JACOBIAN OF E5 EQUATION IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),DFY(MEBND,N),IPAR(*),RPAR(*) C 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 C RETURN END C-------------------------------------------------------------------- subroutine mas(n,am,m,ipar,rpar,ierr) dimension am(2,2),ipar(*),rpar(*) return end C---------------------------------------------------------------------- SUBROUTINE INIT(N,T,Y) DOUBLE PRECISION Y(N),T C Y(1)=1.76D-3 Y(2)=0.0D0 y(3)=0.0D+0 y(4)=0.0D+0 C RETURN END C---------------------------------------------------------------------- SUBROUTINE SOLN(N,T,TRUE,I) DOUBLE PRECISION TRUE(N),T INTEGER N C 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 ENDIF C RETURN END C--------------------------------------------------------------------------