C************************************************************************* C 02/04/2000 C************************************************************************* C DRIVER FOR MEBDFDAE ON OREGONATOR PROBLEM C************************************************************************* IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (ND=3,LWORK=(31+3*ND)*ND+3,LIWORK=ND+14) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK),MBND(4),TRUE(ND), + ERROR(ND),MASBND(4) C CHARACTER*(*) SOLVER, PROBLEM PARAMETER (SOLVER = ' MEBDFDAE',PROBLEM = 'OREGONATOR') C EXTERNAL F,PDERV,MAS C open(6,file='oreg.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... 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=3 C... ENDPOINT OF INTEGRATION X=0.0D0 XEND=360.0D+0 C... INITIAL VALUES CALL INIT(N,X,Y) C...REQUIRED TOLERANCE RTOL=TOLST ATOL=RTOL*1.0D-6 ITOL=2 WRITE(6,*) 'RESULT WITH THE FOLLOWING TOL :' WRITE(6,*) 'RTOL =' ,RTOL WRITE(6,*) 'ATOL =' ,ATOL WRITE(6,*) 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=30.0D+0 masbnd(1)=0 C it1=mclock() C... C... CALL OF THE SUBROUTINE DO 20 I=1,12 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+30.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*atol sumh=sumh*atol/12.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,DY,ipar,rpar,ierr) C... RIGHT-HAND SIDE OF OREGON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DY(N),ipar(*),rpar(*) C 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)) C RETURN END C------------------------------------------------------------------------- SUBROUTINE PDERV(X,Y,DFY,N,MEBND,ipar,rpar,ierr) C... JACOBIAN OF OREGON EQUATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(N),DFY(MEBND,N),ipar(*),rpar(*) C 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 C RETURN END C-------------------------------------------------------------------------- subroutine mas(m,am,n,ipar,rpar,ierr) dimension am(2,2),ipar(*),rpar(*) return end C----------------------------------------------------------------------- SUBROUTINE INIT(N,T,Y) DOUBLE PRECISION Y(N),T INTEGER N C Y(1)=1.0D0 Y(2)=2.0D0 Y(3)=3.0D+0 C RETURN END C----------------------------------------------------------------------- SUBROUTINE SOLN(N,T,TRUE,I) DOUBLE PRECISION TRUE(N),T INTEGER N,I C 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 C RETURN END C------------------------------------------------------------------------