! D02ZAF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d02zafe_mod ! D02ZAF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. USE nag_library, ONLY : nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: iset = 1, itrace = 0, neq = 3, & nin = 5, nout = 6 INTEGER, PARAMETER :: nrw = 50 + 4*neq INTEGER, PARAMETER :: nwkjac = neq*(neq+1) INTEGER, PARAMETER :: ldysav = neq CONTAINS SUBROUTINE fcn(neq,t,y,f,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) :: f(neq) REAL (KIND=nag_wp), INTENT (IN) :: y(neq) ! .. Executable Statements .. f(1) = -0.04E0_nag_wp*y(1) + 1.0E4_nag_wp*y(2)*y(3) f(2) = 0.04E0_nag_wp*y(1) - 1.0E4_nag_wp*y(2)*y(3) - & 3.0E7_nag_wp*y(2)*y(2) f(3) = 3.0E7_nag_wp*y(2)*y(2) RETURN END SUBROUTINE fcn SUBROUTINE jac(neq,t,y,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) ! .. Local Scalars .. REAL (KIND=nag_wp) :: hxd ! .. Executable Statements .. hxd = h*d p(1,1) = 1.0E0_nag_wp - hxd*(-0.04E0_nag_wp) p(1,2) = -hxd*(1.0E4_nag_wp*y(3)) p(1,3) = -hxd*(1.0E4_nag_wp*y(2)) p(2,1) = -hxd*(0.04E0_nag_wp) p(2,2) = 1.0E0_nag_wp - hxd*(-1.0E4_nag_wp*y(3)-6.0E7_nag_wp*y(2)) p(2,3) = -hxd*(-1.0E4_nag_wp*y(2)) ! Do not need to set P(3,1) since Jacobian preset to zero ! P(3,1) = - HXD*(0.0E0) p(3,2) = -hxd*(6.0E7_nag_wp*y(2)) p(3,3) = 1.0E0_nag_wp - hxd*(0.0E0_nag_wp) RETURN END SUBROUTINE jac SUBROUTINE monitr(neq,ldysav,t,hlast,hnext,y,ydot,ysav,r,acor,imon, & inln,hmin,hmax,nqu) ! .. Use Statements .. USE nag_library, ONLY : d02zaf ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: hlast, t REAL (KIND=nag_wp), INTENT (INOUT) :: hmax, hmin, hnext INTEGER, INTENT (INOUT) :: imon INTEGER, INTENT (OUT) :: inln INTEGER, INTENT (IN) :: ldysav, neq, nqu ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: acor(neq,2), r(neq), & ydot(neq), ysav(ldysav,*) REAL (KIND=nag_wp), INTENT (INOUT) :: y(neq) ! .. Local Scalars .. REAL (KIND=nag_wp) :: errloc INTEGER :: i, ifail ! .. Executable Statements .. inln = 3 IF (imon==1) THEN ifail = -1 errloc = d02zaf(neq,acor(1,2),acor(1,1),ifail) IF (ifail/=0) THEN imon = -2 ELSE IF (errloc>5.0E0_nag_wp) THEN WRITE (nout,99999) t, (y(i),i=1,neq), errloc ELSE WRITE (nout,99998) t, (y(i),i=1,neq) END IF END IF RETURN 99999 FORMAT (1X,F10.6,3(F13.7,2X)/1X,' ** WARNING scaled local error = ', & F13.5) 99998 FORMAT (1X,F10.6,3(F13.7,2X)) END SUBROUTINE monitr END MODULE d02zafe_mod PROGRAM d02zafe ! D02ZAF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : d02nbf, d02nsf, d02nvf, d02nyf, nag_wp, x04abf USE d02zafe_mod, ONLY : fcn, iset, itrace, jac, ldysav, monitr, neq, & nin, nout, nrw, nwkjac ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: h, h0, hmax, hmin, hu, t, tcrit, & tcur, tolsf, tout INTEGER :: i, ifail, imxer, itask, itol, & maxord, maxstp, mxhnil, niter, & nje, nq, nqu, nre, nst, outchn, & sdysav LOGICAL :: petzld ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: atol(:), rtol(:), rwork(:), & wkjac(:), y(:), ydot(:), ysav(:,:) REAL (KIND=nag_wp) :: con(6) INTEGER :: inform(23) LOGICAL, ALLOCATABLE :: algequ(:) ! .. Executable Statements .. WRITE (nout,*) 'D02ZAF Example Program Results' ! Skip heading in data file READ (nin,*) ! neq: number of differential equations READ (nin,*) maxord, maxstp, mxhnil sdysav = maxord + 1 ALLOCATE (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq), & ydot(neq),ysav(ldysav,sdysav),algequ(neq)) outchn = nout CALL x04abf(iset,outchn) ! Set algorithmic and problem parameters READ (nin,*) hmin, hmax, h0, t, tout READ (nin,*) petzld ! Initialisation ! Integrate to tout without passing tout. tcrit = tout itask = 4 ! Use default values for the array con. con(1:6) = 0.0_nag_wp ! Use BDF formulae with modified Newton method. ! Use averaged L2 norm for local error control. ifail = 0 CALL d02nvf(neq,sdysav,maxord,'Newton',petzld,con,tcrit,hmin,hmax,h0, & maxstp,mxhnil,'Average-L2',rwork,ifail) ! Setup for using analytic Jacobian ifail = 0 CALL d02nsf(neq,neq,'Analytical',nwkjac,rwork,ifail) WRITE (nout,*) WRITE (nout,*) ' Analytic Jacobian' WRITE (nout,*) ! Set tolerances. READ (nin,*) itol READ (nin,*) rtol(1), atol(1) ! Initial values for Y. READ (nin,*) y(1:neq) WRITE (nout,*) ' X Y(1) Y(2) Y(3)' WRITE (nout,99999) t, (y(i),i=1,neq) ! Solve the problem using MONITR to output intermediate results. ifail = -1 CALL d02nbf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,fcn, & ysav,sdysav,jac,wkjac,nwkjac,monitr,itask,itrace,ifail) IF (ifail==0) THEN ! Get integration statistics. CALL d02nyf(neq,neq,hu,h,tcur,tolsf,rwork,nst,nre,nje,nqu,nq,niter, & imxer,algequ,inform,ifail) WRITE (nout,*) WRITE (nout,99997) hu, h, tcur WRITE (nout,99996) nst, nre, nje WRITE (nout,99995) nqu, nq, niter WRITE (nout,99994) ' Max Err Comp = ', imxer WRITE (nout,*) ELSE WRITE (nout,*) WRITE (nout,99998) 'Exit D02NBF with IFAIL = ', ifail, ' and T = ', & t END IF 99999 FORMAT (1X,F10.6,3(F13.7,2X)) 99998 FORMAT (1X,A,I2,A,E12.5) 99997 FORMAT (1X,' HUSED = ',E12.5,' HNEXT = ',E12.5,' TCUR = ',E12.5) 99996 FORMAT (1X,' NST = ',I6,' NRE = ',I6,' NJE = ',I6) 99995 FORMAT (1X,' NQU = ',I6,' NQ = ',I6,' NITER = ',I6) 99994 FORMAT (1X,A,I4) END PROGRAM d02zafe