* D02NMF Example Program Text * Mark 14 Revised. NAG Copyright 1989. * .. Parameters .. INTEGER NOUT PARAMETER (NOUT=6) INTEGER NEQ, NEQMAX, NRW, NINF, NWKJAC, NJCPVT, MAXORD, + NY2DIM, MAXSTP, MXHNIL PARAMETER (NEQ=3,NEQMAX=NEQ,NRW=50+4*NEQMAX,NINF=23, + NWKJAC=NEQMAX*(NEQMAX+1),NJCPVT=1,MAXORD=5, + NY2DIM=MAXORD+1,MAXSTP=200,MXHNIL=5) INTEGER LACORB, LSAVRB PARAMETER (LACORB=50+NEQMAX,LSAVRB=LACORB+NEQMAX) DOUBLE PRECISION H0, HMAX, HMIN, TCRIT PARAMETER (H0=0.0D0,HMAX=10.0D0,HMIN=1.0D-10,TCRIT=0.0D0) LOGICAL PETZLD PARAMETER (PETZLD=.FALSE.) * .. Local Scalars .. DOUBLE PRECISION H, HLAST, HNEXT, HU, T, TC, TCUR, TOLSF, TOUT, + XOUT INTEGER I, IFAIL, IFLAG, IMON, IMXER, INLN, IOUT, IRES, + IREVCM, ITASK, ITOL, ITRACE, LACOR1, LACOR2, + LACOR3, LSAVR1, LSAVR2, LSAVR3, NITER, NJE, NQ, + NQU, NRE, NST * .. Local Arrays .. DOUBLE PRECISION ATOL(NEQMAX), CONST(6), RTOL(NEQMAX), RWORK(NRW), + WKJAC(NWKJAC), Y(NEQMAX), YDOT(NEQMAX), + YSAVE(NEQMAX,NY2DIM) INTEGER INFORM(NINF), JACPVT(NJCPVT) LOGICAL ALGEQU(NEQMAX) * .. External Subroutines .. EXTERNAL D02NMF, D02NSF, D02NVF, D02NYF, D02XKF, X04ABF * .. Intrinsic Functions .. INTRINSIC DBLE, INT * .. Executable Statements .. WRITE (NOUT,*) 'D02NMF Example Program Results' WRITE (NOUT,*) CALL X04ABF(1,NOUT) * * Integrate to TOUT by overshooting TOUT (ITASK=1) using B.D.F. * formulae with a Newton method. Default values for the array CONST * are used. Employ scalar tolerances and the Jacobian is evaluated * internally. On the reverse communication call equivalent to the * MONITR call in forward communication routines carry out * interpolation using D02XKF. * T = 0.0D0 TOUT = 10.0D0 ITASK = 1 IOUT = 1 XOUT = 2.0D0 Y(1) = 1.0D0 Y(2) = 0.0D0 Y(3) = 0.0D0 ITOL = 1 RTOL(1) = 1.0D-4 ATOL(1) = 1.0D-7 DO 20 I = 1, 6 CONST(I) = 0.0D0 20 CONTINUE IFAIL = 0 * CALL D02NVF(NEQMAX,NY2DIM,MAXORD,'Newton',PETZLD,CONST,TCRIT,HMIN, + HMAX,H0,MAXSTP,MXHNIL,'Average-L2',RWORK,IFAIL) CALL D02NSF(NEQ,NEQMAX,'Numerical',NWKJAC,RWORK,IFAIL) * LACOR1 = LACORB + 1 LACOR2 = LACORB + 2 LACOR3 = LACORB + 3 LSAVR1 = LSAVRB + 1 LSAVR2 = LSAVRB + 2 LSAVR3 = LSAVRB + 3 WRITE (NOUT,*) ' X Y(1) Y(2) Y(3)' WRITE (NOUT,99999) T, (Y(I),I=1,NEQ) * * Soft fail and error messages only IREVCM = 0 ITRACE = 0 40 IFAIL = 1 * CALL D02NMF(NEQ,NEQMAX,T,TOUT,Y,YDOT,RWORK,RTOL,ATOL,ITOL,INFORM, + YSAVE,NY2DIM,WKJAC,NWKJAC,JACPVT,NJCPVT,IMON,INLN, + IRES,IREVCM,ITASK,ITRACE,IFAIL) * IF (IREVCM.NE.0) THEN IF (IREVCM.EQ.1 .OR. IREVCM.EQ.3) THEN * Equivalent to FCN evaluation in forward communication * routines RWORK(LSAVR1) = -0.04D0*Y(1) + 1.0D4*Y(2)*Y(3) RWORK(LSAVR2) = 0.04D0*Y(1) - 1.0D4*Y(2)*Y(3) - 3.0D7*Y(2) + *Y(2) RWORK(LSAVR3) = 3.0D7*Y(2)*Y(2) ELSE IF (IREVCM.EQ.4) THEN * Equivalent to FCN evaluation in forward communication * routines RWORK(LACOR1) = -0.04D0*Y(1) + 1.0D4*Y(2)*Y(3) RWORK(LACOR2) = 0.04D0*Y(1) - 1.0D4*Y(2)*Y(3) - 3.0D7*Y(2) + *Y(2) RWORK(LACOR3) = 3.0D7*Y(2)*Y(2) ELSE IF (IREVCM.EQ.5) THEN * Equivalent to FCN evaluation in forward communication * routines YDOT(1) = -0.04D0*Y(1) + 1.0D4*Y(2)*Y(3) YDOT(2) = 0.04D0*Y(1) - 1.0D4*Y(2)*Y(3) - 3.0D7*Y(2)*Y(2) YDOT(3) = 3.0D7*Y(2)*Y(2) ELSE IF (IREVCM.EQ.9) THEN * Equivalent to MONITR call in forward communication routines IF (IMON.EQ.1) THEN TC = RWORK(19) HLAST = RWORK(15) HNEXT = RWORK(16) NQU = INT(RWORK(10)) 60 CONTINUE IF (TC-HLAST.LT.XOUT .AND. XOUT.LE.TC) THEN IFLAG = 1 * CALL D02XKF(XOUT,RWORK(LSAVR1),NEQ,YSAVE,NEQMAX, + NY2DIM,RWORK(LACOR1),NEQ,TC,NQU,HLAST, + HNEXT,IFLAG) * IF (IFLAG.NE.0) THEN IMON = -2 ELSE WRITE (NOUT,99999) XOUT, (RWORK(LSAVRB+I),I=1,NEQ) IOUT = IOUT + 1 XOUT = DBLE(IOUT)*2.0D0 IF (IOUT.LT.6) GO TO 60 END IF END IF END IF ELSE IF (IREVCM.EQ.2 .OR. IREVCM.EQ.6 .OR. IREVCM.EQ.7 .OR. + IREVCM.EQ.8) THEN WRITE (NOUT,*) WRITE (NOUT,99995) 'Illegal value of IREVCM = ', IREVCM STOP END IF GO TO 40 ELSE IF (IFAIL.EQ.0) THEN * CALL D02NYF(NEQ,NEQMAX,HU,H,TCUR,TOLSF,RWORK,NST,NRE,NJE, + NQU,NQ,NITER,IMXER,ALGEQU,INFORM,IFAIL) * WRITE (NOUT,*) WRITE (NOUT,99997) ' HUSED = ', HU, ' HNEXT = ', H, + ' TCUR = ', TCUR WRITE (NOUT,99996) ' NST = ', NST, ' NRE = ', NRE, + ' NJE = ', NJE WRITE (NOUT,99996) ' NQU = ', NQU, ' NQ = ', NQ, + ' NITER = ', NITER WRITE (NOUT,99995) ' Max err comp = ', IMXER WRITE (NOUT,*) ELSE WRITE (NOUT,*) WRITE (NOUT,99998) 'Exit D02NMF with IFAIL = ', IFAIL, + ' and T = ', T END IF END IF STOP * 99999 FORMAT (1X,F8.3,3(F13.5,2X)) 99998 FORMAT (1X,A,I2,A,D12.5) 99997 FORMAT (1X,A,D12.5,A,D12.5,A,D12.5) 99996 FORMAT (1X,A,I6,A,I6,A,I6) 99995 FORMAT (1X,A,I4) END