* D02MVF Example Program Text * Mark 14 Release. NAG Copyright 1989. * .. Parameters .. INTEGER NOUT PARAMETER (NOUT=6) INTEGER NEQ, NEQMAX, NRW, NINF, NWKJAC, MAXORD, SDYSAV, + MAXSTP, MXHNIL, LDYSAV PARAMETER (NEQ=6,NEQMAX=NEQ,NRW=50+4*NEQMAX,NINF=23, + NWKJAC=NEQMAX*(NEQMAX+1),MAXORD=5, + SDYSAV=MAXORD+3,MAXSTP=5000,MXHNIL=5,LDYSAV=NEQ) DOUBLE PRECISION H0, HMAX, HMIN PARAMETER (H0=0.0D0,HMAX=0.0D0,HMIN=1.0D-10) * .. Local Scalars .. DOUBLE PRECISION DUM, PI, T, TCRIT, TOUT INTEGER I, IFAIL, ITASK, ITOL, ITRACE, OUTCHN * .. Local Arrays .. DOUBLE PRECISION ATOL(NEQ), CON(3), RTOL(NEQ), RWORK(NRW), + WKJAC(NWKJAC), Y(NEQ), YDOT(NEQ), + YSAV(LDYSAV,SDYSAV) INTEGER INFORM(NINF) LOGICAL LDERIV(2) * .. External Functions .. DOUBLE PRECISION X01AAF EXTERNAL X01AAF * .. External Subroutines .. EXTERNAL D02MVF, D02NBY, D02NGF, D02NSF, DAEJAC, DAERES, + X04ABF * .. Executable Statements .. WRITE (NOUT,*) 'D02MVF Example Program Results' OUTCHN = NOUT CALL X04ABF(1,OUTCHN) DO 20 I = 1, 3 CON(I) = 0.0D0 20 CONTINUE ITRACE = 0 PI = X01AAF(DUM) RTOL(1) = 1.0D-3 ATOL(1) = 1.0D-6 ITOL = 1 T = 0.0D0 TOUT = PI * Set initial values and derivatives Y(1) = 1.D0 Y(2) = 0.0D0 Y(3) = 0.D0 Y(4) = 0.D0 Y(5) = 0.0D0 Y(6) = 0.0D0 YDOT(1) = Y(3) - Y(6)*Y(1) YDOT(2) = Y(4) - Y(6)*Y(2) YDOT(3) = -Y(5)*Y(1) YDOT(4) = -Y(5)*Y(2) - 1.D0 YDOT(5) = -3.D0*Y(4) YDOT(6) = 0.0D0 ITASK = 4 TCRIT = TOUT IFAIL = 1 * CALL D02MVF(NEQMAX,SDYSAV,MAXORD,CON,TCRIT,HMIN,HMAX,H0,MAXSTP, + MXHNIL,'AVERAGE-L2',RWORK,IFAIL) * IF (IFAIL.EQ.0) THEN WRITE (NOUT,*) WRITE (NOUT,99999) + 'Pendulum problem with relative tolerance = ', RTOL(1) WRITE (NOUT,99999) + ' and absolute tolerance = ', ATOL(1) WRITE (NOUT,*) WRITE (NOUT,*) + ' t y1 y2 y3 y4 y5 y6' WRITE (NOUT,99998) T, (Y(I),I=1,6) IFAIL = 0 CALL D02NSF(NEQ,NEQMAX,'ANALYTIC',NWKJAC,RWORK,IFAIL) * IFAIL = 1 LDERIV(1) = .TRUE. LDERIV(2) = .TRUE. * CALL D02NGF(NEQ,LDYSAV,T,TOUT,Y,YDOT,RWORK,RTOL,ATOL,ITOL, + INFORM,DAERES,YSAV,SDYSAV,DAEJAC,WKJAC,NWKJAC, + D02NBY,LDERIV,ITASK,ITRACE,IFAIL) * IF (IFAIL.NE.0) THEN WRITE (NOUT,99997) IFAIL, ' at T = ', T ELSE WRITE (NOUT,99998) T, (Y(I),I=1,NEQ) END IF ELSE WRITE (NOUT,99996) IFAIL END IF * 99999 FORMAT (1X,A,1P,E7.1) 99998 FORMAT (1X,F6.4,3X,6(F8.4)) 99997 FORMAT (1X,'Integration terminated with IFAIL = ',I2,A,F6.4) 99996 FORMAT (1X,/1X,' ** D02MVF returned with IFAIL = ',I5) END * SUBROUTINE DAERES(NEQ,T,Y,YDOT,R,IRES) * .. Scalar Arguments .. DOUBLE PRECISION T INTEGER IRES, NEQ * .. Array Arguments .. DOUBLE PRECISION R(NEQ), Y(NEQ), YDOT(NEQ) * .. Executable Statements .. IF (IRES.EQ.-1) THEN R(1) = -YDOT(1) R(2) = -YDOT(2) R(3) = -YDOT(3) R(4) = -YDOT(4) R(5) = 0.0D0 R(6) = 0.0D0 ELSE R(1) = Y(3) - Y(6)*Y(1) - YDOT(1) R(2) = Y(4) - Y(6)*Y(2) - YDOT(2) R(3) = -Y(5)*Y(1) - YDOT(3) R(4) = -Y(5)*Y(2) - 1.D0 - YDOT(4) R(5) = Y(1)*Y(3) + Y(2)*Y(4) R(6) = Y(1)**2 + Y(2)**2 - 1.0D0 END IF RETURN END SUBROUTINE DAEJAC(NEQ,T,Y,YDOT,H,D,P) * .. Parameters .. DOUBLE PRECISION ONE, TWO PARAMETER (ONE=1.0D0,TWO=2.0D0) * .. Scalar Arguments .. DOUBLE PRECISION D, H, T INTEGER NEQ * .. Array Arguments .. DOUBLE PRECISION P(NEQ,NEQ), Y(NEQ), YDOT(NEQ) * .. Local Scalars .. DOUBLE PRECISION HXD * .. Executable Statements .. HXD = H*D P(1,1) = (ONE+HXD*Y(6)) P(1,3) = -HXD P(1,6) = HXD*Y(1) P(2,2) = (ONE+HXD*Y(6)) P(2,4) = -HXD P(2,6) = HXD*Y(2) P(3,1) = HXD*Y(5) P(3,3) = ONE P(3,5) = HXD*Y(1) P(4,2) = HXD*Y(5) P(4,4) = ONE P(4,5) = HXD*Y(2) P(5,1) = -HXD*Y(3) P(5,2) = -HXD*Y(4) P(5,3) = -HXD*Y(1) P(5,4) = -HXD*Y(2) P(6,1) = -TWO*HXD*Y(1) P(6,2) = -TWO*HXD*Y(2) RETURN END