! D02NBF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d02nbfe_mod ! D02NBF 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 :: alpha = 0.04_nag_wp Real (Kind=nag_wp), Parameter :: beta = 1.0E4_nag_wp Real (Kind=nag_wp), Parameter :: gamma = 3.0E7_nag_wp Real (Kind=nag_wp), Parameter :: one = 1.0_nag_wp Real (Kind=nag_wp), Parameter :: two = 2.0_nag_wp 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) ! .. 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) = -alpha*y(1) + beta*y(2)*y(3) f(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2) f(3) = gamma*y(2)*y(2) Return End Subroutine fcn Subroutine jac(neq,t,y,h,d,p) ! .. 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) = one - hxd*(-alpha) p(1,2) = -hxd*(beta*y(3)) p(1,3) = -hxd*(beta*y(2)) p(2,1) = -hxd*(alpha) p(2,2) = one - hxd*(-beta*y(3)-two*gamma*y(2)) p(2,3) = -hxd*(-beta*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*(two*gamma*y(2)) p(3,3) = one - hxd*(0.0_nag_wp) Return End Subroutine jac End Module d02nbfe_mod Program d02nbfe ! D02NBF Example Main Program ! .. Use Statements .. Use nag_library, Only: d02nbf, d02nby, d02nbz, d02nsf, d02nvf, d02nyf, & nag_wp, x04abf Use d02nbfe_mod, Only: fcn, iset, itrace, jac, ldysav, neq, nin, nout, & nrw, nwkjac ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: h, h0, hmax, hmin, hu, t, tcrit, & tcur, tinit, tolsf, tout Integer :: i, icase, 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(:), & yinit(:), ysav(:,:) Real (Kind=nag_wp) :: con(6) Integer :: inform(23) Logical, Allocatable :: algequ(:) ! .. Executable Statements .. Write (nout,*) 'D02NBF Example Program Results' ! Skip heading in data file Read (nin,*) Read (nin,*) maxord, maxstp, mxhnil sdysav = maxord + 1 Allocate (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq), & yinit(neq),ydot(neq),ysav(ldysav,sdysav),algequ(neq)) outchn = nout Call x04abf(iset,outchn) Read (nin,*) petzld Read (nin,*) hmin, hmax, h0 Read (nin,*) tinit, tout Read (nin,*) itol Read (nin,*) yinit(1:neq) Read (nin,*) rtol(1), atol(1) ! Two cases. In both cases: ! integrate to tout without passing tout; ! use B.D.F formulae with a Newton method; ! use default values for the array con; ! use scalar tolerances; ! use NAG dummy routine D02NBY in place of MONITR subroutine. con(1:6) = 0.0_nag_wp tcrit = tout itask = 4 cases: Do icase = 1, 2 t = tinit y(1:neq) = yinit(1:neq) ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call d02nvf(neq,sdysav,maxord,'Newton',petzld,con,tcrit,hmin,hmax,h0, & maxstp,mxhnil,'Average-L2',rwork,ifail) Write (nout,*) ifail = 0 Select Case (icase) Case (1) ! First case. The Jacobian is evaluated internally. Call d02nsf(neq,neq,'Numerical',nwkjac,rwork,ifail) Write (nout,*) ' Numerical Jacobian' Case (2) ! Second case. The Jacobian is evaluated by jac. Call d02nsf(neq,neq,'Analytical',nwkjac,rwork,ifail) Write (nout,*) ' Analytic Jacobian' End Select Write (nout,99993)(i,i=1,neq) Write (nout,99999) t, y(1:neq) ! Soft fail and error messages only ifail = -1 Select Case (icase) Case (1) Call d02nbf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, & fcn,ysav,sdysav,d02nbz,wkjac,nwkjac,d02nby,itask,itrace,ifail) Case (2) Call d02nbf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform, & fcn,ysav,sdysav,jac,wkjac,nwkjac,d02nby,itask,itrace,ifail) End Select If (ifail==0) Then Write (nout,99999) t, y(1:neq) ifail = 0 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 End Do cases 99999 Format (1X,F8.3,3(F13.5,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) 99993 Format (/1X,' X ',3(' Y(',I1,') ')) End Program d02nbfe