! D02MVF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d02mvfe_mod ! D02MVF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. USE nag_library, ONLY : nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. REAL (KIND=nag_wp), PARAMETER :: one = 1.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: two = 2.0_nag_wp INTEGER, PARAMETER :: iset = 1, neq = 6, nin = 5, & nout = 6 INTEGER, PARAMETER :: nrw = 50 + 4*neq INTEGER, PARAMETER :: nwkjac = neq*(neq+1) INTEGER, PARAMETER :: sdysav = 8 INTEGER, PARAMETER :: ldysav = neq CONTAINS SUBROUTINE resid(neq,t,y,ydot,r,ires) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t INTEGER, INTENT (INOUT) :: ires INTEGER, INTENT (IN) :: neq ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: r(neq) REAL (KIND=nag_wp), INTENT (IN) :: y(neq), ydot(neq) ! .. Executable Statements .. IF (ires==-1) THEN r(1) = -ydot(1) r(2) = -ydot(2) r(3) = -ydot(3) r(4) = -ydot(4) r(5:6) = 0.0E0_nag_wp 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) - one - ydot(4) r(5) = y(1)*y(3) + y(2)*y(4) r(6) = y(1)**2 + y(2)**2 - one END IF RETURN END SUBROUTINE resid SUBROUTINE jac(neq,t,y,ydot,h,d,p) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: d, h, t INTEGER, INTENT (IN) :: neq ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (INOUT) :: p(neq,neq) REAL (KIND=nag_wp), INTENT (IN) :: y(neq), ydot(neq) ! .. Local Scalars .. REAL (KIND=nag_wp) :: 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 SUBROUTINE jac END MODULE d02mvfe_mod PROGRAM d02mvfe ! D02MVF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : d02mvf, d02nby, d02ngf, d02nsf, nag_wp, x04abf USE d02mvfe_mod, ONLY : iset, jac, ldysav, neq, nin, nout, nrw, nwkjac, & one, resid, sdysav ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: h0, hmax, hmin, t, tcrit, tout INTEGER :: i, ifail, itask, itol, itrace, & maxord, maxstp, mxhnil, outchn ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: atol(:), rtol(:), rwork(:), & wkjac(:), y(:), ydot(:), ysav(:,:) REAL (KIND=nag_wp) :: con(3) INTEGER :: inform(23) LOGICAL :: lderiv(2) ! .. Executable Statements .. WRITE (nout,*) 'D02MVF Example Program Results' ! Skip heading in data file READ (nin,*) ! Allocations based on number of differential equations (neq) ALLOCATE (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq), & ydot(neq),ysav(ldysav,sdysav)) ! Read algorithmic parameters READ (nin,*) maxord, maxstp, mxhnil READ (nin,*) hmin, hmax, h0 READ (nin,*) rtol(1), atol(1) READ (nin,*) itrace, itol READ (nin,*) t, tout READ (nin,*) y(1:neq) READ (nin,*) itask ! Set initial derivatives and other values outchn = nout CALL x04abf(iset,outchn) con(1:3) = 0.0_nag_wp 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) - one ydot(5) = -3.0_nag_wp*y(4) ydot(6) = 0.0_nag_wp tcrit = tout ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 CALL d02mvf(neq,sdysav,maxord,con,tcrit,hmin,hmax,h0,maxstp,mxhnil, & 'AVERAGE-L2',rwork,ifail) WRITE (nout,*) WRITE (nout,99999) 'Pendulum problem with relative tolerance', rtol(1) WRITE (nout,99999) ' and absolute tolerance', atol(1) WRITE (nout,99998) (i,i=1,neq) WRITE (nout,99997) t, y(1:neq) ifail = 0 CALL d02nsf(neq,neq,'ANALYTIC',nwkjac,rwork,ifail) lderiv(1) = .TRUE. lderiv(2) = .TRUE. ifail = 0 CALL d02ngf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,resid, & ysav,sdysav,jac,wkjac,nwkjac,d02nby,lderiv,itask,itrace,ifail) WRITE (nout,99997) t, y(1:neq) 99999 FORMAT (1X,A40,' =',1P,E8.1) 99998 FORMAT (/1X,' t ',3X,6(' y',I1)) 99997 FORMAT (1X,F7.4,2X,6(F8.4)) END PROGRAM d02mvfe